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.

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
Show file tree
Hide file tree
Changes from all commits
Commits
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
33 changes: 15 additions & 18 deletions pyresample/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
from pyresample import CHUNK_SIZE
from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP
from pyresample.area_config import create_area_def
from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary
from pyresample.boundary import Boundary, SimpleBoundary
from pyresample.utils import check_slice_orientation, load_cf_area
from pyresample.utils.proj4 import (
get_geostationary_height,
Expand Down Expand Up @@ -2595,14 +2595,10 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None):
y_slice = _ensure_integer_slice(y_slice)
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))
area_boundary = self._get_area_to_cover_boundary(area_to_cover)
intersection = data_boundary.contour_poly.intersection(area_boundary.contour_poly)
data_boundary = _get_area_boundary(self)
area_boundary = _get_area_boundary(area_to_cover)
intersection = data_boundary.contour_poly.intersection(
area_boundary.contour_poly)
if intersection is None:
logger.debug('Cannot determine appropriate slicing. '
"Data and projection area do not overlap.")
Expand All @@ -2622,15 +2618,6 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None):
return (check_slice_orientation(x_slice),
check_slice_orientation(y_slice))

@staticmethod
def _get_area_to_cover_boundary(area_to_cover: AreaDefinition) -> Boundary:
try:
if area_to_cover.is_geostationary:
return Boundary(*get_geostationary_bounding_box_in_lonlats(area_to_cover))
return AreaDefBoundary(area_to_cover, 100)
except ValueError:
raise NotImplementedError("Can't determine boundary of area to cover")

def crop_around(self, other_area):
"""Crop this area around `other_area`."""
xslice, yslice = self.get_area_slices(other_area)
Expand Down Expand Up @@ -2730,6 +2717,16 @@ def geocentric_resolution(self, ellps='WGS84', radius=None):
return res


def _get_area_boundary(area_to_cover: AreaDefinition) -> Boundary:
try:
if area_to_cover.is_geostationary:
return Boundary(*get_geostationary_bounding_box_in_lonlats(area_to_cover))
boundary_shape = max(max(*area_to_cover.shape) // 100 + 1, 3)
return area_to_cover.boundary(frequency=boundary_shape, force_clockwise=True)
except ValueError:
raise NotImplementedError("Can't determine boundary of area to cover")


def _make_slice_divisible(sli, max_size, factor=2):
"""Make the given slice even in size."""
rem = (sli.stop - sli.start) % factor
Expand Down
12 changes: 12 additions & 0 deletions pyresample/test/test_geometry/test_area.py
Original file line number Diff line number Diff line change
Expand Up @@ -1821,6 +1821,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)

def test_area_to_cover_all_nan_bounds(self, geos_src_area, create_test_area):
"""Check area slicing when the target doesn't have a valid boundary."""
area_def = geos_src_area
Expand Down