pytroll / pyresample

Geospatial image resampling in Python
http://pyresample.readthedocs.org
GNU Lesser General Public License v3.0
344 stars 95 forks source link

Preserve get_area_slices behavior when area to cover has an invalid boundary #524

Closed djhoese closed 1 year ago

djhoese commented 1 year ago

As described by @lobsiger on the mailing list, Satpy now fails resampling (specifically reducing data before resampling) because of changes in #465. The main problem is that areas that used to "work" would result in a NotImplementedError inside AreaDefinition.get_area_slices because no intersection could be determined with the source area and nothing "smarter" is currently implemented. This no intersection determination was actually because the target area ("area to cover") had inf in its boundary coordinates and the math just worked out that no intersection was possible. The notImplementedError is a signal to Satpy to not attempt data reduction as no area subset slices could be generated...but that its OK that this happened.

After #465, the infinity values in the coordinates is now raising a ValueError. This PR catches this and produces a NotImplementedError to mimic the previous behavior.

An alternative solution would be to update Satpy to catch not only the NotImplementedError, but ValueError and maybe RuntimeError as possible reasons to ignore/not attempt data reduction.

@mraspaud @pnuu @ghiggi thoughts?

codecov[bot] commented 1 year ago

Codecov Report

Merging #524 (0fa73f0) into main (7a434c1) will increase coverage by 0.01%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main     #524      +/-   ##
==========================================
+ Coverage   94.32%   94.34%   +0.01%     
==========================================
  Files          79       79              
  Lines       12976    12987      +11     
==========================================
+ Hits        12240    12252      +12     
+ Misses        736      735       -1     
Flag Coverage Δ
unittests 94.34% <100.00%> (+0.01%) :arrow_up:

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files Coverage Δ
pyresample/geometry.py 87.56% <100.00%> (+0.05%) :arrow_up:
pyresample/test/test_geometry/test_area.py 99.20% <100.00%> (+<0.01%) :arrow_up:

... and 1 file with indirect coverage changes

coveralls commented 1 year ago

Coverage Status

coverage: 93.941% (+0.005%) from 93.936% when pulling 0fa73f03d6540fb173254c904271748431fe8b4a on djhoese:bugfix-noboundary-areaslices into 7a434c116a256b155c75db2ec99765abdd02035d on pytroll:main.

ghiggi commented 1 year ago

Hey @djhoese

This looks like a good temporary fix.

I thought that in the medium-term I think would be nice to do some reorg within pyresample related to the extraction of boundary sides coordinates and boundary object creation.

For the boundary sides coordinates, we currently use get_bbox_lonlats and get_geostationary_bounding_box_in_lonlats to deal with Inf in projection coordinates of geostationary areas.

I guess that we could refactor somewhat the code in get_bbox_lonlats to:

  1. call the get_geostationary_bounding_box_in_lonlats (when an area is geostationary) within such a method
  2. and if the area is a full globe projection that has Inf popping up in lon/lat coordinates (after conversion from projection coordinates resulting out of Earth disk), we search within the area lons/lats array the most exterior valid coordinates. This operation would take a bit more computations (I guess all-in-memory lon/lats arrays) but would enable all downstream computations.

As a result, we could just call get_bbox_lonlats throughout the codebase.

And as a second step, we could standardize the creation of Boundary objects. Currently, we use :

If we would use the area.boundary() method calling the refactored get_bbox_lonlats, we would remove all the if-else logics currently reappearing many times throughout the code.

djhoese commented 1 year ago

@ghiggi That makes a lot of sense. I like it. I'm not sure why we didn't move the geos-specific stuff into the get_bbox_lonlats in the first place. So I think we agree that this little hack/change of mine gives us a quick fix for backwards compatibility. The get_bbox_lonlats deserves a rewrite and should be smart enough to raise a ValueError or RuntimeError if a "decent" set of boundary lon/lats can be produced. Or should the Boundary classes produce the error if infinities appear in the bounding coordinates?

mraspaud commented 1 year ago

I think a more specific error would be nice, so that Satpy can catch it. Eg BooleanComputationError, so that satpy can handle this specific case without risk for catching something unrelated. But as a temporary fix this is acceptable.

mraspaud commented 1 year ago

Merging this. @djhoese do you want to add issue to act on the points mentioned above?

djhoese commented 1 year ago

I'll make an issue for the overall conversation @ghiggi is talking about. For the more specific error, I'm not sure Satpy should care in this case. We can't determine the slices so (satpy) stop worrying about it. Internally pyresample might care, but I think maybe this comes after some of the suggestions @ghiggi is talking about get implemented.