Scene.crop fails for ll_bbox near edge of geostationary source area

[djhoese]

Describe the bug A Geo2Grid user is cropping a GOES-17 ABI CONUS image with a ll_bbox using the Scene.crop method. When he uses the bounds [-130, 30, -113, 51] everything is fine and he gets a good image. When he uses [-130, 30, -113, 52] then he ends up with a strip of data (<150 rows when it should be ~3000).

To Reproduce

Download some data: https://noaa-goes17.s3.amazonaws.com/ABI-L1b-RadC/2020/258/18/OR_ABI-L1b-RadC-M6C02_G17_s20202581801176_e20202581803549_c20202581803580.nc

 from satpy import Scene
scn = Scene(reader='abi_l1b', filenames='OR_ABI-L1b-RadC-M6C02_G17_s20202581801176_e20202581803549_c20202581803580.nc')
scn.load(['C02'])
crop_scn = scn.crop(ll_bbox=[-130, 30, -113, 52])
print(crop_scn['C02'].attrs['area'])

Result:

Area ID: GOES-West
Description: 0.5km at nadir
Projection ID: abi_fixed_grid
Projection: {'ellps': 'GRS80', 'h': '35786023', 'lon_0': '-137', 'no_defs': 'None', 'proj': 'geos', 'sweep': 'x', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 2967
Number of rows: 73
Area extent: (661826.7094, 3059132.3901, 2148306.5327, 3095705.7056)

Expected behavior A better estimation of the crop region so that I have about the same size region regardless of 1 degree of latitude.

Actual results See above.

Screenshots If applicable, add screenshots to help explain your problem.

Environment Info:

• OS: all 3
• Satpy Version: master
• PyResample Version: master

Additional context

So I've tracked this down to two parts: the size of the lon/lat AreaDefinition we create to perform ll_bbox slicing and the number of points used to estimate the boundary around this area.

We create a lon/lat area of size 100 by 100:

https://github.com/pytroll/satpy/blob/5b45c32fe0600bd658ce74f75a67898e17fb202d/satpy/scene.py#L477-L479

We then have pyresample estimate the boundary of this by taking a coordinate for every 100 points:

https://github.com/pytroll/pyresample/blob/46de898ec00a6bae243af18ba34dd8f3e53d06f2/pyresample/geometry.py#L2036

area_boundary = AreaDefBoundary(area_to_cover, 100)

Increasing the number of rows/columns in the AreaDefinition or decreasing the frequency of coordinates we take in the boundary improves the estimation here. However, that also decreases performance.

Options I see:

1. Increase the static size of the area definition.
2. Decrease the frequency.
3. Make the area definition the same number of rows/columns as the source area (tried this, terrible for performance, gets worse with larger source areas).
4. Make the frequency a percentage (10%?) of the source area, not a number of pixels.

[djhoese]

Given the small change in latitude bounds I wonder if this is a bug in the intersection code instead?

[djhoese]

The user let me know that this particular box is outside the CONUS region (I haven't mapped it out). I'd still say this is pretty unexpected.

[mraspaud]

I'm wondering it that's not an expected behaviour. The user asks to crop the scene to match the ll box given. If the ll box is partly outside the scene, I expect crop to just return the intersection. From what you say, the "faulty" ll box is like that. Or am I missing something?

[djhoese]

I thought it was outside the projection space, but from what I can tell that 52 latitude is in the space at both of those longitudes. So the box, at these corner points, should be valid. Maybe the boundary is choosing different points (in the middle of the "box") that aren't in the space?