Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow cropping non-geos areas #516

Merged
merged 4 commits into from
Jun 28, 2023
Merged

Allow cropping non-geos areas #516

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4e245a9
Allow cropping non-geos areas
mraspaud May 23, 2023
f325262
Merge branch 'main' into feature-non-geos-crop
mraspaud Jun 28, 2023
814ad27
Refactor area boundary generation
mraspaud Jun 28, 2023
44af636
Make static method a function
mraspaud Jun 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 3 additions & 5 deletions pyresample/geometry.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2587,11 +2587,9 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None):
return x_slice, y_slice

if not self.is_geostationary:
raise NotImplementedError("Source projection must be 'geos' if "
"source/target projections are not "
"equal.")

data_boundary = Boundary(*get_geostationary_bounding_box_in_lonlats(self))
data_boundary = self.boundary(frequency=100)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

else:
data_boundary = Boundary(*get_geostationary_bounding_box_in_lonlats(self))
Copy link
Contributor

@ghiggi ghiggi May 25, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 ...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this simplification is better suited for another PR

if area_to_cover.is_geostationary:
area_boundary = Boundary(
*get_geostationary_bounding_box_in_lonlats(area_to_cover))
Expand Down
12 changes: 12 additions & 0 deletions pyresample/test/test_geometry/test_area.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1814,6 +1814,18 @@ def test_on_flipped_geos_area(self, create_test_area):
assert slice_lines == expected_slice_lines
assert slice_cols == expected_slice_cols

def test_non_geos_can_be_cropped(self, create_test_area):
"""Test that non-geos areas can be cropped also."""
src_area = create_test_area(dict(proj="utm", zone=33),
10980, 10980,
(499980.0, 6490200.0, 609780.0, 6600000.0))
crop_area = create_test_area({'proj': 'latlong'},
100, 100,
(15.9689, 58.5284, 16.4346, 58.6995))
slice_x, slice_y = src_area.get_area_slices(crop_area)
assert slice_x == slice(5630, 8339)
assert slice_y == slice(9261, 10980)


class TestBoundary:
"""Test 'boundary' method for AreaDefinition classes."""
Expand Down