Allow cropping non-geos areas #516

mraspaud · 2023-05-23T13:44:40Z

Tests added
Tests passed
Passes git diff origin/main **/*py | flake8 --diff
Fully documented

mraspaud · 2023-05-23T13:47:00Z

@ameraner you might want to have a look at this also

codecov · 2023-05-23T14:05:42Z

Codecov Report

Merging #516 (44af636) into main (2a6d769) will decrease coverage by 0.10%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main     #516      +/-   ##
==========================================
- Coverage   94.34%   94.25%   -0.10%     
==========================================
  Files          82       82              
  Lines       13026    13030       +4     
==========================================
- Hits        12290    12282       -8     
- Misses        736      748      +12

Flag	Coverage Δ
unittests	`94.25% <100.00%> (-0.10%)`	⬇️

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

Impacted Files	Coverage Δ
pyresample/geometry.py	`87.64% <100.00%> (+0.05%)`	⬆️
pyresample/test/test_geometry/test_area.py	`99.21% <100.00%> (+<0.01%)`	⬆️

... and 1 file with indirect coverage changes

djhoese · 2023-05-23T14:45:18Z

If I'm not mistaken this will mean that reduce_data in Satpy will now cause a performance hit for all LEO data right? For some reason I thought reduce_data depended on this failing to not have reduce_data be performed on LEO/SwathDefinition data.

djhoese · 2023-05-23T14:46:48Z

Yep, this NotImplementedError handling here:

https://github.com/pytroll/satpy/blob/4f6a7c35ce8a620428f51ff086e3276a73c72c97/satpy/scene.py#L903-L926

ameraner

LGTM, thank you for "unlocking" this useful feature! I left just a clarification question.

On a side note, the code for SwathDef data now falls into here

pyresample/pyresample/geometry.py

Line 2597 in 4e245a9

area_boundary = AreaDefBoundary(area_to_cover, 100)

which then reaches a PendingDeprecationWarning

pyresample/pyresample/boundary.py

Lines 119 to 131 in 3fa9529

    
           class AreaDefBoundary(AreaBoundary): 
        
               """Boundaries for area definitions (pyresample).""" 
        
               def __init__(self, area, frequency=1): 
        
                   lons, lats = area.get_bbox_lonlats() 
        
                   warnings.warn("'AreaDefBoundary' will be removed in the future. " + 
        
                                 "Use the Swath/AreaDefinition 'boundary' method instead!.", 
        
                                 PendingDeprecationWarning, stacklevel=2) 
        
                   AreaBoundary.__init__(self, 
        
                                         *zip(lons, lats)) 
        
                   if frequency != 1: 
        
                       self.decimate(frequency)

ameraner · 2023-05-24T07:15:34Z

pyresample/geometry.py

-                                      "equal.")
-
-        data_boundary = Boundary(*get_geostationary_bounding_box_in_lonlats(self))
+            data_boundary = self.boundary(frequency=100)


what is the rationale for choosing 100 here? (we also hardcode it a bit later

pyresample/pyresample/geometry.py

Line 2597 in 4e245a9

area_boundary = AreaDefBoundary(area_to_cover, 100)

). I see that the geos case has some more handling of the frequency here

pyresample/pyresample/geometry.py

Lines 1559 to 1569 in 4e245a9

def _get_geo_boundary_sides(self, frequency=None):

"""Retrieve the boundary sides list for geostationary projections."""

# Define default frequency

if frequency is None:

frequency = 50

# Ensure at least 4 points are used

if frequency < 4:

frequency = 4

# Ensure an even number of vertices for side creation

if (frequency % 2) != 0:

frequency = frequency + 1

.

We set that default just to restrict the number of vertices (several thousands)... otherwise the polygon intersection would takes too long

Yes, we want to keep the number of vertices low to speed up the process.

mraspaud · 2023-05-24T13:10:17Z

If I'm not mistaken this will mean that reduce_data in Satpy will now cause a performance hit for all LEO data right? For some reason I thought reduce_data depended on this failing to not have reduce_data be performed on LEO/SwathDefinition data.

So maybe we should check that the it is an areadef not swathdef here? or should that be done in satpy to keeps things cleaner here?

djhoese · 2023-05-25T02:26:40Z

Oh I'm wrong. There is an AreaDefinition check at the top of the method. This is probably good then.

ghiggi · 2023-05-25T09:15:18Z

@mraspaud: FYI there are several edge cases where this method currently fails:

When a swath/area intersects twice with another area, the polygon intersection routine currently runs infinitely if the intersections are more than 1
When an area is within (or touches) another, I think that the intersection now runs infinitely ...

If you implement get_array_indices_from_lonlat method also for SwathDefinition, then the method would work also when area_to_cover is SwathDefinition ;)

ghiggi · 2023-05-25T09:26:43Z

pyresample/geometry.py

-        data_boundary = Boundary(*get_geostationary_bounding_box_in_lonlats(self))
+            data_boundary = self.boundary(frequency=100)
+        else:
+            data_boundary = Boundary(*get_geostationary_bounding_box_in_lonlats(self))


If we add/call get_geostationary_bounding_box_in_lonlats in https://github.com/mraspaud/pyresample/blob/4e245a906bf3aeff072f6b61febdd1df298fb425/pyresample/geometry.py#L299 get_bbox_lonlats when is geostationary you can simplify all the if_else conditions here just do self.boundary() with a sensible frequency defaults ...

Here below:

AreaDefBoundary(area_to_cover, 100) fails if area_to_cover is SwathDefinition.
Maybe use area_to_cover.boundary() also here in all cases (with sensible frequency defaults)

I think this simplification is better suited for another PR

mraspaud · 2023-05-25T11:10:28Z

Oh I'm wrong. There is an AreaDefinition check at the top of the method. This is probably good then.

I can't see it, could you point me to it please?

djhoese · 2023-05-26T01:33:51Z

Here:

pyresample/pyresample/geometry.py

Lines 2570 to 2571 in 3fa9529

    
           if not isinstance(area_to_cover, AreaDefinition): 
        
               raise NotImplementedError('Only AreaDefinitions can be used')

djhoese · 2023-06-27T02:04:34Z

Hhhmmm maybe I should have merged this before my big config PR. Sorry, but I've caused some conflicts.

mraspaud · 2023-06-28T10:07:53Z

Fixed.

coveralls · 2023-06-28T10:43:41Z

coverage: 93.864% (-0.09%) from 93.953% when pulling 44af636 on mraspaud:feature-non-geos-crop into 2a6d769 on pytroll:main.

mraspaud · 2023-06-28T11:40:08Z

@djhoese please review my two last commits to see if you agree with the refactoring I made.

This refactoring took me some time because I found a bug/inconsistency in the usage of the frequency argument in AreaDefBoundary and the boundary method of area defintion: in the first case it means the step we use to determine the length of the boundary, in the second case it is the length of the boundary itself. I'm working on a fix now (which will be another PR).

ghiggi · 2023-06-28T12:06:15Z

I remember I noticed such inconsistencies already a while ago. For GEO areas is the total number of vertices, while for swath and other areas is the number of steps on each boundary side right? I personally would really like to have this as the number of vertices per boundary side, and allow to specify a tuple (y,x) or (x,y) instead of a single value for all dimensions.

I also remember somewhat that quite ago, by using the step approach, we were discarding boundary corners values, but I think this was fixed by one of mine or your PR @mraspaud.

mraspaud · 2023-06-28T12:10:23Z

@ghiggi I think both approaches are valid, the usage just depends on what the user needs at that point.

in #526 I have now just replaced frequency with bbox_shape as a more accurate name. I don't change any functionality or behaviour though, so allowing to pass a tuple for example isn't supported (yet?)

mraspaud · 2023-06-28T12:37:29Z

which then reaches a PendingDeprecationWarning

Indeed, this is now fixed.

djhoese · 2023-06-28T18:40:12Z

@mraspaud I think I understand the last 2 commits. I assume the main motivation was to:

Remove use of AreaDefBoundary
Use one method for getting the boundaries for both the source and cover(target) areas.

So if that is true then 👍 from me.

Allow cropping non-geos areas

4e245a9

mraspaud added the enhancement label May 23, 2023

mraspaud requested a review from djhoese May 23, 2023 13:44

mraspaud self-assigned this May 23, 2023

ameraner approved these changes May 24, 2023

View reviewed changes

ghiggi reviewed May 25, 2023

View reviewed changes

djhoese added this to the v1.28.0 milestone Jun 19, 2023

Merge branch 'main' into feature-non-geos-crop

f325262

mraspaud added 2 commits June 28, 2023 13:35

Refactor area boundary generation

814ad27

Make static method a function

44af636

mraspaud merged commit 2048b12 into pytroll:main Jun 28, 2023

mraspaud deleted the feature-non-geos-crop branch June 28, 2023 20:48

djhoese mentioned this pull request Nov 24, 2023

Question about cropping resampled scenes pytroll/tutorial-satpy-half-day#13

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Allow cropping non-geos areas #516

Allow cropping non-geos areas #516

mraspaud commented May 23, 2023

mraspaud commented May 23, 2023

codecov bot commented May 23, 2023 •

edited

Loading

djhoese commented May 23, 2023

djhoese commented May 23, 2023

ameraner left a comment

ameraner May 24, 2023

ghiggi May 25, 2023

mraspaud Jun 28, 2023

mraspaud commented May 24, 2023

djhoese commented May 25, 2023

ghiggi commented May 25, 2023 •

edited

Loading

ghiggi May 25, 2023 •

edited

Loading

ghiggi May 25, 2023

mraspaud May 25, 2023

mraspaud commented May 25, 2023

djhoese commented May 26, 2023

djhoese commented Jun 27, 2023

mraspaud commented Jun 28, 2023

coveralls commented Jun 28, 2023 •

edited

Loading

mraspaud commented Jun 28, 2023

ghiggi commented Jun 28, 2023

mraspaud commented Jun 28, 2023

mraspaud commented Jun 28, 2023

djhoese commented Jun 28, 2023

	class AreaDefBoundary(AreaBoundary):
	"""Boundaries for area definitions (pyresample)."""

	def __init__(self, area, frequency=1):
	lons, lats = area.get_bbox_lonlats()
	warnings.warn("'AreaDefBoundary' will be removed in the future. " +
	"Use the Swath/AreaDefinition 'boundary' method instead!.",
	PendingDeprecationWarning, stacklevel=2)
	AreaBoundary.__init__(self,
	*zip(lons, lats))

	if frequency != 1:
	self.decimate(frequency)

	def _get_geo_boundary_sides(self, frequency=None):
	"""Retrieve the boundary sides list for geostationary projections."""
	# Define default frequency
	if frequency is None:
	frequency = 50
	# Ensure at least 4 points are used
	if frequency < 4:
	frequency = 4
	# Ensure an even number of vertices for side creation
	if (frequency % 2) != 0:
	frequency = frequency + 1

Allow cropping non-geos areas #516

Allow cropping non-geos areas #516

Conversation

mraspaud commented May 23, 2023

mraspaud commented May 23, 2023

codecov bot commented May 23, 2023 • edited Loading

Codecov Report

djhoese commented May 23, 2023

djhoese commented May 23, 2023

ameraner left a comment

Choose a reason for hiding this comment

ameraner May 24, 2023

Choose a reason for hiding this comment

ghiggi May 25, 2023

Choose a reason for hiding this comment

mraspaud Jun 28, 2023

Choose a reason for hiding this comment

mraspaud commented May 24, 2023

djhoese commented May 25, 2023

ghiggi commented May 25, 2023 • edited Loading

ghiggi May 25, 2023 • edited Loading

Choose a reason for hiding this comment

ghiggi May 25, 2023

Choose a reason for hiding this comment

mraspaud May 25, 2023

Choose a reason for hiding this comment

mraspaud commented May 25, 2023

djhoese commented May 26, 2023

djhoese commented Jun 27, 2023

mraspaud commented Jun 28, 2023

coveralls commented Jun 28, 2023 • edited Loading

mraspaud commented Jun 28, 2023

ghiggi commented Jun 28, 2023

mraspaud commented Jun 28, 2023

mraspaud commented Jun 28, 2023

djhoese commented Jun 28, 2023

codecov bot commented May 23, 2023 •

edited

Loading

ghiggi commented May 25, 2023 •

edited

Loading

ghiggi May 25, 2023 •

edited

Loading

coveralls commented Jun 28, 2023 •

edited

Loading