From 915e34f2ac851bee93eda66c19c0b73e00562981 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Mon, 27 Nov 2023 22:56:07 +0100 Subject: [PATCH] Fix boundary ordering logic for swath and projections ! --- docs/source/howtos/spherical_geometry.rst | 4 +- pyresample/_config.py | 1 + pyresample/boundary/area_boundary.py | 3 + pyresample/boundary/base_boundary.py | 19 ++- pyresample/boundary/geographic_boundary.py | 150 ++++++++++++------ pyresample/boundary/projection_boundary.py | 35 ++-- pyresample/boundary/utils.py | 84 ++++++++++ pyresample/future/geometry/_subset.py | 6 +- pyresample/future/geometry/area.py | 1 + pyresample/geometry.py | 148 +++++------------ pyresample/slicer.py | 30 ++-- .../test_boundary/test_geographic_boundary.py | 73 ++++----- pyresample/test/test_boundary/test_order.py | 133 ++++++++++++++++ pyresample/test/test_boundary/test_utils.py | 145 +++++++++++++++++ pyresample/test/test_geometry/test_area.py | 65 ++++---- pyresample/test/test_geometry/test_swath.py | 2 +- pyresample/test/test_slicer.py | 11 +- 17 files changed, 636 insertions(+), 274 deletions(-) create mode 100644 pyresample/boundary/utils.py create mode 100644 pyresample/test/test_boundary/test_order.py create mode 100644 pyresample/test/test_boundary/test_utils.py diff --git a/docs/source/howtos/spherical_geometry.rst b/docs/source/howtos/spherical_geometry.rst index 0bf413470..6ca535671 100644 --- a/docs/source/howtos/spherical_geometry.rst +++ b/docs/source/howtos/spherical_geometry.rst @@ -74,11 +74,11 @@ satellite passes. See trollschedule_ how to generate a list of satellite overpas >>> from pyresample.spherical_utils import GetNonOverlapUnions >>> area_boundary = area_def.geographic_boundary(vertices_per_side=100) # doctest: +SKIP - >>> area_boundary = area_boundary.contour_poly # doctest: +SKIP + >>> area_boundary = area_boundary.polygon # doctest: +SKIP >>> list_of_polygons = [] >>> for mypass in passes: # doctest: +SKIP - >>> list_of_polygons.append(mypass.boundary.contour_poly) # doctest: +SKIP + >>> list_of_polygons.append(mypass.geographic_boundary().polygon) # doctest: +SKIP >>> non_overlaps = GetNonOverlapUnions(list_of_polygons) # doctest: +SKIP >>> non_overlaps.merge() # doctest: +SKIP diff --git a/pyresample/_config.py b/pyresample/_config.py index 9a4a2b1ae..4d92c927b 100644 --- a/pyresample/_config.py +++ b/pyresample/_config.py @@ -40,6 +40,7 @@ "features": { "future_geometries": False, }, + "force_boundary_computations": False, }], paths=_CONFIG_PATHS, ) diff --git a/pyresample/boundary/area_boundary.py b/pyresample/boundary/area_boundary.py index 5c9a8f9f6..b3087b829 100644 --- a/pyresample/boundary/area_boundary.py +++ b/pyresample/boundary/area_boundary.py @@ -43,6 +43,9 @@ def contour(self): @property def contour_poly(self): """Get the Spherical polygon corresponding to the Boundary.""" + warnings.warn("'contour_poly' is deprecated." + + "Use the 'boundary().polygon' property instead!.", + PendingDeprecationWarning, stacklevel=2) if self._contour_poly is None: self._contour_poly = SphPolygon( np.deg2rad(np.vstack(self.contour()).T)) diff --git a/pyresample/boundary/base_boundary.py b/pyresample/boundary/base_boundary.py index ffeb8f61e..7306237fc 100644 --- a/pyresample/boundary/base_boundary.py +++ b/pyresample/boundary/base_boundary.py @@ -30,18 +30,25 @@ class BaseBoundary: """Base class for boundary objects.""" __slots__ = ["_sides_x", "_sides_y"] - def __init__(self, sides_x, sides_y, order=None): + def __init__(self, area, vertices_per_side=None): + + sides_x, sides_y = self._compute_boundary_sides(area, vertices_per_side) self._sides_x = BoundarySides(sides_x) self._sides_y = BoundarySides(sides_y) + self._area = area - self.is_clockwise = self._check_is_boundary_clockwise(sides_x, sides_y) + self.is_clockwise = self._check_is_boundary_clockwise(sides_x, sides_y, area) self.is_counterclockwise = not self.is_clockwise - self._set_order(order) + self._set_order(order=None) # FIX ! - def _check_is_boundary_clockwise(self, sides_x, sides_y): + def _check_is_boundary_clockwise(self, sides_x, sides_y, area): """Check if the boundary is clockwise or counterclockwise.""" raise NotImplementedError() + def _compute_boundary_sides(self, area, vertices_per_side): + """Compute boundary sides.""" + raise NotImplementedError() + def _set_order(self, order): """Set the order of the boundary vertices.""" if self.is_clockwise: @@ -112,9 +119,9 @@ def contour(self, closed=False): y = np.hstack((y, y[0])) return x, y - def _to_shapely_polygon(self): + def to_shapely_polygon(self): """Define a Shapely Polygon.""" from shapely.geometry import Polygon - self = self.set_counterclockwise() # FIXME: add exception for pole wrapping polygons + self = self.set_counterclockwise() x, y = self.contour(closed=True) return Polygon(zip(x, y)) diff --git a/pyresample/boundary/geographic_boundary.py b/pyresample/boundary/geographic_boundary.py index abe6ebb88..3adc13510 100644 --- a/pyresample/boundary/geographic_boundary.py +++ b/pyresample/boundary/geographic_boundary.py @@ -24,44 +24,86 @@ from pyresample.boundary.area_boundary import Boundary as OldBoundary from pyresample.boundary.base_boundary import BaseBoundary +from pyresample.spherical import SphPolygon logger = logging.getLogger(__name__) -def _is_corner_is_clockwise(lon1, lat1, corner_lon, corner_lat, lon2, lat2): - """Determine if coordinates follow a clockwise path. +def _get_swath_local_inside_point_from_sides(sides_x, sides_y, start_idx=0): + """Retrieve point inside boundary close to the point used to determine order. + + This is required for global areas (spanning more than 180 degrees) and swaths. + The start_idx indicates the opposites sides corners (start_idx, start_idx+2) from + which to try identify a point inside the polygon. + """ + # NOTE: the name of the object here refers to start_idx=0 + from pyresample.spherical import Arc, SCoordinate + idx1 = start_idx + idx2 = (start_idx + 2) % 4 + + top_corner = sides_x[idx1][0], sides_y[idx1][0] + top_right_point = sides_x[idx1][1], sides_y[idx1][1] + bottom_corner = sides_x[idx2][-1], sides_y[idx2][-1] + bottom_right_point = sides_x[idx2][-2], sides_y[idx2][-2] + point_top_corner = SCoordinate(*np.deg2rad(top_corner)) + point_top_right_point = SCoordinate(*np.deg2rad(top_right_point)) + point_bottom_corner = SCoordinate(*np.deg2rad(bottom_corner)) + point_bottom_right_point = SCoordinate(*np.deg2rad(bottom_right_point)) + arc1 = Arc(point_top_corner, point_bottom_right_point) + arc2 = Arc(point_top_right_point, point_bottom_corner) + point_inside = arc1.intersection(arc2) + if point_inside is not None: + point_inside = point_inside.vertices_in_degrees[0] + return point_inside + + +def _try_get_local_inside_point(sides_x, sides_y): + """Try to get a local inside point from one of the 4 boundary sides corners.""" + for start_idx in range(0, 4): + point_inside = _get_swath_local_inside_point_from_sides(sides_x, sides_y, start_idx=start_idx) + if point_inside is not None: + return point_inside, start_idx + else: + return None, None + + +def _is_clockwise_order(first_point, second_point, point_inside): + """Determine if polygon coordinates follow a clockwise path. This uses :class:`pyresample.spherical.Arc` to determine the angle - between the first line segment (Arc) from (lon1, lat1) to - (corner_lon, corner_lat) and the second line segment from - (corner_lon, corner_lat) to (lon2, lat2). A straight line would - produce an angle of 0, a clockwise path would have a negative angle, - and a counter-clockwise path would have a positive angle. + between a polygon arc segment and a point known to be inside the polygon. + Note: pyresample.spherical assumes angles are positive if counterclockwise. + Note: if the longitude distance between the first_point/second_point and point_inside is + larger than 180°, the function likely return a wrong unexpected result ! """ - import math - from pyresample.spherical import Arc, SCoordinate - point1 = SCoordinate(math.radians(lon1), math.radians(lat1)) - point2 = SCoordinate(math.radians(corner_lon), math.radians(corner_lat)) - point3 = SCoordinate(math.radians(lon2), math.radians(lat2)) - arc1 = Arc(point1, point2) - arc2 = Arc(point2, point3) - angle = arc1.angle(arc2) + point1 = SCoordinate(*np.deg2rad(first_point)) + point2 = SCoordinate(*np.deg2rad(second_point)) + point3 = SCoordinate(*np.deg2rad(point_inside)) + arc12 = Arc(point1, point2) + arc23 = Arc(point2, point3) + angle = arc12.angle(arc23) is_clockwise = -np.pi < angle < 0 return is_clockwise -def _is_boundary_clockwise(sides_lons, sides_lats): - """Determine if the boundary sides are clockwise.""" - is_clockwise = _is_corner_is_clockwise( - lon1=sides_lons[0][-2], - lat1=sides_lats[0][-2], - corner_lon=sides_lons[0][-1], - corner_lat=sides_lats[0][-1], - lon2=sides_lons[1][1], - lat2=sides_lats[1][1]) - return is_clockwise +def _check_is_clockwise(area, sides_x, sides_y): + from pyresample import SwathDefinition + + if isinstance(area, SwathDefinition): + point_inside, start_idx = _try_get_local_inside_point(sides_x, sides_y) + first_point = sides_x[start_idx][0], sides_y[start_idx][0] + second_point = sides_x[start_idx][1], sides_y[start_idx][1] + return _is_clockwise_order(first_point, second_point, point_inside) + else: + if area.is_geostationary: + point_inside = area.get_lonlat(row=int(area.shape[0] / 2), col=int(area.shape[1] / 2)) + first_point = sides_x[0][0], sides_y[0][0] + second_point = sides_x[0][1], sides_y[0][1] + return _is_clockwise_order(first_point, second_point, point_inside) + else: + return True class GeographicBoundary(BaseBoundary, OldBoundary): @@ -74,17 +116,35 @@ class GeographicBoundary(BaseBoundary, OldBoundary): # from the old interface for compatibility to AreaBoundary @classmethod - def _check_is_boundary_clockwise(cls, sides_x, sides_y): + def _check_is_boundary_clockwise(cls, sides_x, sides_y, area): """GeographicBoundary specific implementation.""" - return _is_boundary_clockwise(sides_lons=sides_x, sides_lats=sides_y) + return _check_is_clockwise(area, sides_x, sides_y) + + @classmethod + def _compute_boundary_sides(cls, area, vertices_per_side): + sides_lons, sides_lats = area._get_geographic_sides(vertices_per_side=vertices_per_side) + return sides_lons, sides_lats - def __init__(self, sides_lons, sides_lats, order=None, crs=None): - super().__init__(sides_x=sides_lons, sides_y=sides_lats, order=order) + def __init__(self, area, vertices_per_side=None): + super().__init__(area=area, vertices_per_side=vertices_per_side) self.sides_lons = self._sides_x self.sides_lats = self._sides_y - self.crs = crs or CRS(proj="longlat", ellps="WGS84") - self._contour_poly = None # Backcompatibility with old AreaBoundary + + # Define CRS + if self.is_swath: + crs = self._area.crs + else: + crs = CRS(proj="longlat", ellps="WGS84") # FIXME: AreaDefinition.get_lonlat for geographic projections? + self.crs = crs + + # Backcompatibility with old AreaBoundary + self._contour_poly = None + + @property + def is_swath(self): + """Determine if is the boundary of a swath.""" + return self._area.__class__.__name__ == "SwathDefinition" @property def lons(self): @@ -96,28 +156,22 @@ def lats(self): """Retrieve boundary latitude vertices.""" return self._y - def _to_spherical_polygon(self): - self = self.set_clockwise() # TODO: add exception for pole wrapping polygons - raise NotImplementedError("This will return a SPolygon in pyresample 2.0") - - def polygon(self, shapely=False): - """Return the boundary polygon.""" - # ALTERNATIVE: - # - shapely: to_shapely_polygon(), to_shapely_line(), - # - pyresample spherical: to_polygon(), to_line(), polygon, line - if shapely: - return self._to_shapely_polygon() - else: - return self._to_spherical_polygon() - - def plot(self, ax=None, subplot_kw=None, **kwargs): + @property + def polygon(self): + """Return the boundary spherical polygon.""" + self = self.set_clockwise() + if self._contour_poly is None: + self._contour_poly = SphPolygon(np.deg2rad(self.vertices)) + return self._contour_poly + + def plot(self, ax=None, subplot_kw=None, alpha=0.6, **kwargs): """Plot the the boundary.""" import cartopy.crs as ccrs from pyresample.visualization.geometries import plot_geometries - geom = self.polygon(shapely=True) + geom = self.to_shapely_polygon() crs = ccrs.Geodetic() p = plot_geometries(geometries=[geom], crs=crs, - ax=ax, subplot_kw=subplot_kw, **kwargs) + ax=ax, subplot_kw=subplot_kw, alpha=alpha, **kwargs) return p diff --git a/pyresample/boundary/projection_boundary.py b/pyresample/boundary/projection_boundary.py index 7620f50ff..c1fd07055 100644 --- a/pyresample/boundary/projection_boundary.py +++ b/pyresample/boundary/projection_boundary.py @@ -45,17 +45,27 @@ class ProjectionBoundary(BaseBoundary): """ @classmethod - def _check_is_boundary_clockwise(cls, sides_x, sides_y): + def _check_is_boundary_clockwise(cls, sides_x, sides_y, area=None): """GeographicBoundary specific implementation.""" return _is_projection_boundary_clockwise(sides_x=sides_x, sides_y=sides_y) - def __init__(self, sides_x, sides_y, crs, order=None, cartopy_crs=None): - super().__init__(sides_x=sides_x, sides_y=sides_y, order=order) + @classmethod + def _compute_boundary_sides(cls, area, vertices_per_side): + sides_x, sides_y = area._get_projection_sides(vertices_per_side=vertices_per_side) + return sides_x, sides_y + + def __init__(self, area, vertices_per_side=None): + super().__init__(area=area, vertices_per_side=vertices_per_side) self.sides_x = self._sides_x self.sides_y = self._sides_y - self.crs = crs - self.cartopy_crs = cartopy_crs + self.crs = self._area.crs + self.cartopy_crs = self._area.to_cartopy_crs() + + @property + def _point_inside(self): + x, y = self._area.get_proj_coords(data_slice=(int(self.shape[0] / 2), int(self.shape[1] / 2))) + return x.squeeze(), y.squeeze() @property def x(self): @@ -67,14 +77,7 @@ def y(self): """Retrieve boundary y vertices.""" return self._y - def polygon(self, shapely=True): - """Return the boundary polygon.""" - if shapely: - return self._to_shapely_polygon() - else: - raise NotImplementedError("Only shapely polygon available.") - - def plot(self, ax=None, subplot_kw=None, crs=None, **kwargs): + def plot(self, ax=None, subplot_kw=None, crs=None, alpha=0.6, **kwargs): """Plot the the boundary. crs must be a Cartopy CRS !""" from pyresample.visualization.geometries import plot_geometries @@ -82,8 +85,10 @@ def plot(self, ax=None, subplot_kw=None, crs=None, **kwargs): raise ValueError("Projection Cartopy 'crs' is required to display projection boundary.") if crs is None: crs = self.cartopy_crs + if subplot_kw is None: + subplot_kw = {"projection": crs} - geom = self.polygon(shapely=True) + geom = self.to_shapely_polygon() p = plot_geometries(geometries=[geom], crs=crs, - ax=ax, subplot_kw=subplot_kw, **kwargs) + ax=ax, subplot_kw=subplot_kw, alpha=alpha, **kwargs) return p diff --git a/pyresample/boundary/utils.py b/pyresample/boundary/utils.py new file mode 100644 index 000000000..f6e425e2f --- /dev/null +++ b/pyresample/boundary/utils.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (c) 2014-2023 Pyresample developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +"""Utility to extract boundary mask and indices.""" +import numpy as np + + +def _find_boundary_mask(mask): + """Find boundary of a binary mask.""" + mask = mask.astype(int) + # Pad with zeros (enable to detect Trues at mask boundaries) + padded_mask = np.pad(mask, ((1, 1), (1, 1)), mode='constant', constant_values=0) + + # Shift the image in four directions and compare with the original + shift_up = np.roll(padded_mask, -1, axis=0) + shift_down = np.roll(padded_mask, 1, axis=0) + shift_left = np.roll(padded_mask, -1, axis=1) + shift_right = np.roll(padded_mask, 1, axis=1) + + # Find the boundary points + padded_boundary_mask = ((padded_mask != shift_up) | (padded_mask != shift_down) | + (padded_mask != shift_left) | (padded_mask != shift_right)) & padded_mask + + boundary_mask = padded_boundary_mask[1:-1,1:-1] + return boundary_mask + + +def find_boundary_mask(lons, lats): + """Find the boundary mask.""" + valid_mask = np.isfinite(lons) & np.isfinite(lats) + return _find_boundary_mask(valid_mask) + + +def get_ordered_contour(contour_mask): + """Return the ordered indices of a contour mask.""" + # Count number of rows and columns + rows, cols = contour_mask.shape + # Function to find the next contour point + def next_point(current, last, visited): + for dx, dy in [(-1, 0), (0, 1), (1, 0), (0, -1), (1, -1), (-1, -1), (1, 1), (-1, 1)]: + next_pt = (current[0] + dx, current[1] + dy) + if next_pt != last and next_pt not in visited and 0 <= next_pt[0] < rows and 0 <= next_pt[1] < cols and contour_mask[next_pt]: + return next_pt + return None + # Initialize + contour = [] + visited = set() # Keep track of visited points + # Find the starting point + start = tuple(np.argwhere(contour_mask)[0]) + contour.append(start) + visited.add(start) + # Initialize last and current points + last = start + current = next_point(last, None, visited) + while current and current != start: + contour.append(current) + visited.add(current) + next_pt = next_point(current, last, visited) + if not next_pt: # Break if no next point found + break + last, current = current, next_pt + return np.array(contour) + + +def find_boundary_contour_indices(lons, lats): + """Find the boundary contour ordered indices.""" + boundary_mask = find_boundary_mask(lons, lats) + boundary_contour_idx = get_ordered_contour(boundary_mask) + return boundary_contour_idx + \ No newline at end of file diff --git a/pyresample/future/geometry/_subset.py b/pyresample/future/geometry/_subset.py index c401f54e5..c46374489 100644 --- a/pyresample/future/geometry/_subset.py +++ b/pyresample/future/geometry/_subset.py @@ -49,8 +49,8 @@ def get_area_slices( data_boundary = _get_area_boundary(src_area) area_boundary = _get_area_boundary(area_to_cover) - intersection = data_boundary.contour_poly.intersection( - area_boundary.contour_poly) + intersection = data_boundary.polygon.intersection( + area_boundary.polygon) if intersection is None: logger.debug('Cannot determine appropriate slicing. ' "Data and projection area do not overlap.") @@ -103,7 +103,7 @@ def _get_area_boundary(area_to_cover: AreaDefinition) -> GeographicBoundary: vertices_per_side = None else: vertices_per_side = max(max(*area_to_cover.shape) // 100 + 1, 3) - return area_to_cover.geographic_boundary(vertices_per_side=vertices_per_side, order="clockwise") + return area_to_cover.geographic_boundary(vertices_per_side=vertices_per_side) except ValueError as err: raise NotImplementedError("Can't determine boundary of area to cover") from err diff --git a/pyresample/future/geometry/area.py b/pyresample/future/geometry/area.py index 6042debae..caab271b1 100644 --- a/pyresample/future/geometry/area.py +++ b/pyresample/future/geometry/area.py @@ -26,6 +26,7 @@ from pyresample.geometry import AreaDefinition as LegacyAreaDefinition # noqa from pyresample.geometry import ( # noqa DynamicAreaDefinition, + _get_geostationary_bounding_box, _get_geostationary_bounding_box_in_lonlats, get_full_geostationary_bounding_box_in_proj_coords, get_geostationary_angle_extent, diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 1ada96580..bf9489f1c 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -34,7 +34,7 @@ from pyproj import Geod, Proj from pyproj.aoi import AreaOfUse -from pyresample import CHUNK_SIZE +from pyresample import CHUNK_SIZE, config from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP from pyresample.area_config import create_area_def from pyresample.boundary import SimpleBoundary @@ -108,7 +108,7 @@ def __init__(self, lons=None, lats=None, nprocs=1): self.ndim = None self.cartesian_coords = None self.hash = None - self._boundary_mask = None + self._boundary_contour_idx = None def __getitem__(self, key): """Slice a 2D geographic definition.""" @@ -251,7 +251,8 @@ def get_lonlats(self, data_slice=None, chunks=None, **kwargs): # lons/lats are xarray DataArray objects, use numpy/dask array underneath lons = lons.data lats = lats.data - + # TODO + # --> if data_slice and chunks provided, why here first rechunk all array and then subset? if chunks is not None: import dask.array as da if isinstance(lons, da.Array): @@ -411,10 +412,11 @@ def _get_geostationary_boundary_sides(self, vertices_per_side=None, coordinates= def _compute_boundary_mask(self): """Compute valid boundary mask for AreaDefinition(s) with sides out of the Earth disk.""" - if self._boundary_mask is None: + from pyresample.boundary.utils import find_boundary_contour_indices + if self._boundary_contour_idx is None: lons, lats = self.get_lonlats() # all in memory ! - self._boundary_mask = find_boundary_mask(lons, lats) - return self._boundary_mask + self._boundary_contour_idx = find_boundary_contour_indices(lons, lats) + return self._boundary_contour_idx def _get_geographic_sides(self, vertices_per_side: Optional[int] = None) -> tuple: """Return the geographic boundary sides of the current area. @@ -434,6 +436,8 @@ def _get_geographic_sides(self, vertices_per_side: Optional[int] = None) -> tupl Each list element is a numpy array representing a specific side of the geometry. The order of the sides are [top", "right", "bottom", "left"] """ + if len(self.lons.shape) == 1: + raise ValueError("The area must have 2 dimensions to retrieve the boundary sides.") is_swath = self.__class__.__name__ == "SwathDefinition" if not is_swath and _is_any_corner_out_of_earth_disk(self): # Geostationary @@ -442,14 +446,16 @@ def _get_geographic_sides(self, vertices_per_side: Optional[int] = None) -> tupl coordinates="geographic") # Polar Projections, Global Planar Projections (Mollweide, Robinson) # - Retrieve dummy right and left sides - boundary_mask = self._compute_boundary_mask() - lons, lats = self.get_lonlats() - lons = lons[boundary_mask] - lats = lats[boundary_mask] - sides_lons = self._get_dummy_sides(lons, vertices_per_side=vertices_per_side) - sides_lats = self._get_dummy_sides(lats, vertices_per_side=vertices_per_side) - else: - sides_lons, sides_lats = self._get_sides(coord_fun=self.get_lonlats, vertices_per_side=vertices_per_side) + if config.get("force_boundary_computations", False): + boundary_contour_idx = self._compute_boundary_mask() + lons, lats = self.get_lonlats() + lons = lons[boundary_contour_idx[:, 0], boundary_contour_idx[:, 1]] + lats = lats[boundary_contour_idx[:, 0], boundary_contour_idx[:, 1]] + sides_lons = self._get_dummy_sides(lons, vertices_per_side=vertices_per_side) + sides_lats = self._get_dummy_sides(lats, vertices_per_side=vertices_per_side) + return sides_lons, sides_lats + + sides_lons, sides_lats = self._get_sides(coord_fun=self.get_lonlats, vertices_per_side=vertices_per_side) return sides_lons, sides_lats def _get_sides(self, coord_fun, vertices_per_side): @@ -567,7 +573,8 @@ def boundary(self, *, vertices_per_side=None, force_clockwise=False, frequency=N (i.e. full disc geostationary area, Robinson projection, polar projections, ...) by default only 50 points are selected. force_clockwise: - Perform minimal checks and reordering of coordinates to ensure + DEPRECATED. + Performed minimal checks and reordering of coordinates to ensure that the returned coordinates follow a clockwise direction. This is important for compatibility with :class:`pyresample.spherical.SphPolygon` where operations depend @@ -581,13 +588,9 @@ def boundary(self, *, vertices_per_side=None, force_clockwise=False, frequency=N warnings.warn("The `boundary` method is pending deprecation. Use `geographic_boundary` instead", PendingDeprecationWarning, stacklevel=2) vertices_per_side = vertices_per_side or frequency - if force_clockwise: - order = "clockwise" - else: - order = None - return self.geographic_boundary(vertices_per_side=vertices_per_side, order=order) + return self.geographic_boundary(vertices_per_side=vertices_per_side) - def geographic_boundary(self, vertices_per_side=None, order=None): + def geographic_boundary(self, vertices_per_side=None): """Retrieve the GeographicBoundary object. Parameters @@ -598,25 +601,9 @@ def geographic_boundary(self, vertices_per_side=None, order=None): If the area object is an AreaDefinition with any corner out of the Earth disk (i.e. full disc geostationary area, Robinson projection, polar projections, ...) by default only 50 points are selected. - order: - Specify the desired order of the boundary polygon vertices in GeographicBoundary. - If order=None, the sides order is used. - If order="clockwise", the boundary polygon vertices are returned by - GeographicBoundary in clockwise order. - If order="counterclockwise", the boundary polygon vertices are returned by - GeographicBoundary in counterclockwise order. """ from pyresample.boundary import GeographicBoundary - - sides_lons, sides_lats = self._get_geographic_sides(vertices_per_side=vertices_per_side) - if self.__class__.__name__ == "SwathDefinition": - crs = self.crs - else: - crs = None # default to WGS84 for AreaDefinition - return GeographicBoundary(sides_lons=sides_lons, - sides_lats=sides_lats, - order=order, - crs=crs) + return GeographicBoundary(area=self, vertices_per_side=vertices_per_side) def get_cartesian_coords(self, nprocs=None, data_slice=None, cache=False): """Retrieve cartesian coordinates of geometry definition. @@ -1761,18 +1748,20 @@ def _get_projection_sides(self, vertices_per_side: Optional[int] = None) -> tupl coordinates="projection") # Polar Projections, Global Planar Projections (Mollweide, Robinson) # - Retrieve dummy right and left sides - boundary_mask = self._compute_boundary_mask() - x, y = self.get_proj_coords() - x = x[boundary_mask] - y = y[boundary_mask] - sides_x = self._get_dummy_sides(x, vertices_per_side=vertices_per_side) - sides_y = self._get_dummy_sides(y, vertices_per_side=vertices_per_side) - else: - sides_x, sides_y = self._get_sides(coord_fun=self.get_proj_coords, - vertices_per_side=vertices_per_side) + if config.get("force_boundary_computations", False): + boundary_contour_idx = self._compute_boundary_mask() + x, y = self.get_proj_coords() + x = x[boundary_contour_idx[:, 0], boundary_contour_idx[:, 1]] + y = y[boundary_contour_idx[:, 0], boundary_contour_idx[:, 1]] + sides_x = self._get_dummy_sides(x, vertices_per_side=vertices_per_side) + sides_y = self._get_dummy_sides(y, vertices_per_side=vertices_per_side) + return sides_x, sides_y + + sides_x, sides_y = self._get_sides(coord_fun=self.get_proj_coords, + vertices_per_side=vertices_per_side) return sides_x, sides_y - def projection_boundary(self, vertices_per_side=None, order=None): + def projection_boundary(self, vertices_per_side=None): """Retrieve the ProjectionBoundary object. Parameters @@ -1782,25 +1771,13 @@ def projection_boundary(self, vertices_per_side=None, order=None): By default (None) the full width and height will be provided. If the area object is an AreaDefinition with any corner out of the Earth disk (i.e. full disc geostationary area, Robinson projection, polar projections, ...) - by default only 50 points are selected. - order: - Specify the desired order of the boundary polygon vertices in GeographicBoundary. - If order=None, the sides order is used. - If order="clockwise", the boundary polygon vertices are returned by - GeographicBoundary in clockwise order. - If order="counterclockwise", the boundary polygon vertices are returned by - GeographicBoundary in counterclockwise order. + by default only 50 points are selected.. """ from pyresample.boundary import ProjectionBoundary if self.crs.is_geographic: - return self.geographic_boundary(vertices_per_side=vertices_per_side, order=order) - sides_x, sides_y = self._get_projection_sides(vertices_per_side=vertices_per_side) - return ProjectionBoundary(sides_x=sides_x, - sides_y=sides_y, - crs=self.crs, - order=order, - cartopy_crs=self.to_cartopy_crs() - ) + return self.geographic_boundary(vertices_per_side=vertices_per_side) + return ProjectionBoundary(area=self, + vertices_per_side=vertices_per_side) def get_edge_bbox_in_projection_coordinates(self, vertices_per_side: Optional[int] = None, frequency: Optional[int] = None): @@ -2905,6 +2882,8 @@ def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50): try: x, y = intersection.boundary.xy except NotImplementedError: + # geos_area is fully out of Earth disk + # --> FIXME: Why we do not raise an error here ? return np.array([]), np.array([]) return np.asanyarray(x[:-1]), np.asanyarray(y[:-1]) @@ -3015,47 +2994,6 @@ def _is_any_corner_out_of_earth_disk(area_def): return True -def _find_boundary_indices(mask): - """Assume mask does not have any row/column with only False values.""" - # For rows - first_valid_row = np.argmax(mask, axis=0) - last_valid_row = mask.shape[0] - 1 - np.argmax(mask[::-1], axis=0) - # For columns - first_valid_col = np.argmax(mask, axis=1) - last_valid_col = mask.shape[1] - 1 - np.argmax(mask[:, ::-1], axis=1) - - return first_valid_row, last_valid_row, first_valid_col, last_valid_col - - -def find_boundary_mask(lons, lats): - """Find the boundary mask.""" - valid_mask = np.isfinite(lons) & np.isfinite(lats) - - # Filter out rows and columns without valid values - valid_rows = np.any(valid_mask, axis=1) - valid_cols = np.any(valid_mask, axis=0) - filtered_mask = valid_mask[valid_rows][:, valid_cols] - - # Find boundary indices - fvr, lvr, fvc, lvc = _find_boundary_indices(filtered_mask) - - # Prepare row and column indices for gathering coordinates - col_indices = np.arange(filtered_mask.shape[1]) - row_indices = np.arange(filtered_mask.shape[0]) - - # Gather coordinates using advanced indexing - filtered_boundary_mask = np.zeros_like(filtered_mask, dtype=bool) - filtered_boundary_mask[fvr, col_indices] = True - filtered_boundary_mask[lvr, col_indices] = True - filtered_boundary_mask[row_indices, fvc] = True - filtered_boundary_mask[row_indices, lvc] = True - - # Reinsert the boundary mask into the original mask's shape - boundary_mask = np.zeros_like(valid_mask, dtype=bool) - boundary_mask[np.ix_(valid_rows, valid_cols)] = filtered_boundary_mask - return boundary_mask - - def combine_area_extents_vertical(area1, area2): """Combine the area extents of areas 1 and 2.""" if (area1.area_extent[0] == area2.area_extent[0] and area1.area_extent[2] == area2.area_extent[2]): diff --git a/pyresample/slicer.py b/pyresample/slicer.py index fcc6f8099..09f7c7ab4 100644 --- a/pyresample/slicer.py +++ b/pyresample/slicer.py @@ -128,7 +128,7 @@ def _get_chunk_polygons_for_swath_to_crop(swath_to_crop): line_slice = expand_slice(line_slice) col_slice = expand_slice(col_slice) smaller_swath = swath_to_crop[line_slice, col_slice] - smaller_poly = smaller_swath.geographic_boundary(vertices_per_side=10).polygon(shapely=True) + smaller_poly = smaller_swath.geographic_boundary(vertices_per_side=10).to_shapely_polygon() res.append((smaller_poly, (line_slice, col_slice))) return res @@ -148,8 +148,9 @@ def get_polygon_to_contain(self): if self.area_to_crop.is_geostationary: geo_boundary = self.area_to_crop.projection_boundary(vertices_per_side=360) x_geos, y_geos = geo_boundary.contour(closed=True) - # POSSIBLE BUG: Here I expect that some coordinates could be NaN ! + # POSSIBLE BUG: some coordinates could be NaN ! # - if points of the geostationary disk are out of the CRS bounds of the area_to_contain + # - the order could not be counterclockwise (as expected by shapely) x_geos, y_geos = self._transformer.transform(x_geos, y_geos, direction=TransformDirection.INVERSE) geos_poly = Polygon(zip(x_geos, y_geos)) poly = Polygon(zip(x, y)) @@ -158,10 +159,14 @@ def get_polygon_to_contain(self): raise IncompatibleAreas('No slice on area.') x, y = zip(*poly.exterior.coords) + # Return poly_to_contain (in area_to_crop CRS) return Polygon(zip(*self._transformer.transform(x, y))) def get_slices_from_polygon(self, poly_to_contain): - """Get the slices based on the polygon.""" + """Get the slices based on the polygon boundary of area_to_contain. + + poly_to_contain is expected to have been projected in the CRS of area_to_crop + """ if not poly_to_contain.is_valid: raise IncompatibleAreas("Area outside of domain.") try: @@ -170,22 +175,22 @@ def get_slices_from_polygon(self, poly_to_contain): buffer_size = np.max(self.area_to_contain.resolution) else: buffer_size = 0 - buffered_poly = poly_to_contain.buffer(buffer_size) - bounds = buffered_poly.bounds + buffered_poly_to_contain = poly_to_contain.buffer(buffer_size) + bounds_poly_to_contain = buffered_poly_to_contain.bounds except ValueError as err: raise InvalidArea("Invalid area") from err - poly_to_crop = self.area_to_crop.projection_boundary(vertices_per_side=10).polygon(shapely=True) - if not poly_to_crop.intersects(buffered_poly): + poly_to_crop = self.area_to_crop.projection_boundary(vertices_per_side=10).to_shapely_polygon() + if not poly_to_crop.intersects(buffered_poly_to_contain): raise IncompatibleAreas("Areas not overlapping.") - bounds = self._sanitize_polygon_bounds(bounds) - slice_x, slice_y = self._create_slices_from_bounds(bounds) + x_bounds, y_bounds = self._sanitize_polygon_bounds(bounds_poly_to_contain) + slice_x, slice_y = self._create_slices_from_bounds(x_bounds, y_bounds) return slice_x, slice_y - def _sanitize_polygon_bounds(self, bounds): + def _sanitize_polygon_bounds(self, bounds_poly_to_contain): """Reset the bounds within the shape of the area.""" try: - (minx, miny, maxx, maxy) = bounds + (minx, miny, maxx, maxy) = bounds_poly_to_contain except ValueError as err: raise IncompatibleAreas('No slice on area.') from err x_bounds, y_bounds = self.area_to_crop.get_array_coordinates_from_projection_coordinates(np.array([minx, maxx]), @@ -196,9 +201,8 @@ def _sanitize_polygon_bounds(self, bounds): return x_bounds, y_bounds @staticmethod - def _create_slices_from_bounds(bounds): + def _create_slices_from_bounds(x_bounds, y_bounds): """Create slices from bounds.""" - x_bounds, y_bounds = bounds try: slice_x = slice(int(np.floor(max(np.min(x_bounds), 0))), int(np.ceil(np.max(x_bounds)))) diff --git a/pyresample/test/test_boundary/test_geographic_boundary.py b/pyresample/test/test_boundary/test_geographic_boundary.py index 0af54eac7..b61be7a22 100644 --- a/pyresample/test/test_boundary/test_geographic_boundary.py +++ b/pyresample/test/test_boundary/test_geographic_boundary.py @@ -17,19 +17,18 @@ # You should have received a copy of the GNU Lesser General Public License # along with this program. If not, see . """Test the GeographicBoundary objects.""" -import unittest import numpy as np import pytest +from pyresample import SwathDefinition from pyresample.boundary import GeographicBoundary -class TestGeographicBoundary(unittest.TestCase): - """Test 'GeographicBoundary' class.""" +class TestSwathGeographicBoundary(): + """Test 'GeographicBoundary' class for SwathDefinition.""" - def test_creation(self): - """Test GeographicBoundary creation.""" + def setup_method(self): sides_lons = [np.array([1.0, 1.5, 2.0]), np.array([2.0, 3.0]), np.array([3.0, 3.5, 4.0]), @@ -38,41 +37,40 @@ def test_creation(self): np.array([7.0, 8.0]), np.array([8.0, 8.5, 9.0]), np.array([9.0, 6.0])] + lons = np.array([[1.0, 1.5, 2.0], + [4.0, 3.5, 3.0]]) + lats = np.array([[6.0, 6.5, 7.0], + [9.0, 8.5, 8.0]]) + + self.lons = lons + self.lats = lats + self.area = SwathDefinition(lons, lats) + self.sides_lons = sides_lons + self.sides_lats = sides_lats + self.point_inside = (2.5, 7.5) + def test_creation(self): + """Test GeographicBoundary creation.""" # Define GeographicBoundary - boundary = GeographicBoundary(sides_lons, sides_lats) + boundary = GeographicBoundary(self.area) # Assert sides coincides - for b_lon, src_lon in zip(boundary.sides_lons, sides_lons): + for b_lon, src_lon in zip(boundary.sides_lons, self.sides_lons): assert np.allclose(b_lon, src_lon) - for b_lat, src_lat in zip(boundary.sides_lats, sides_lats): + for b_lat, src_lat in zip(boundary.sides_lats, self.sides_lats): assert np.allclose(b_lat, src_lat) - def test_number_sides_required(self): - """Test GeographicBoundary requires 4 sides .""" - sides_lons = [np.array([1.0, 1.5, 2.0]), - np.array([2.0, 3.0]), - np.array([4.0, 1.0])] - sides_lats = [np.array([6.0, 6.5, 7.0]), - np.array([7.0, 8.0]), - np.array([9.0, 6.0])] + def test_minimum_swath_size(self): + """Test swath has minimum 2 dimensions.""" + area = SwathDefinition(self.lons[0], self.lats[1]) with pytest.raises(ValueError): - GeographicBoundary(sides_lons, sides_lats) + GeographicBoundary(area) def test_vertices_property(self): """Test GeographicBoundary vertices property.""" - sides_lons = [np.array([1.0, 1.5, 2.0]), - np.array([2.0, 3.0]), - np.array([3.0, 3.5, 4.0]), - np.array([4.0, 1.0])] - sides_lats = [np.array([6.0, 6.5, 7.0]), - np.array([7.0, 8.0]), - np.array([8.0, 8.5, 9.0]), - np.array([9.0, 6.0])] # Define GeographicBoundary - boundary = GeographicBoundary(sides_lons, sides_lats) - + boundary = GeographicBoundary(self.area) # Assert vertices expected_vertices = np.array([[1., 6.], [1.5, 6.5], @@ -84,32 +82,17 @@ def test_vertices_property(self): def test_contour(self): """Test that GeographicBoundary.contour(closed=False) returns the correct (lon,lat) tuple.""" - sides_lons = [np.array([1.0, 1.5, 2.0]), - np.array([2.0, 3.0]), - np.array([3.0, 3.5, 4.0]), - np.array([4.0, 1.0])] - sides_lats = [np.array([6.0, 6.5, 7.0]), - np.array([7.0, 8.0]), - np.array([8.0, 8.5, 9.0]), - np.array([9.0, 6.0])] # Define GeographicBoundary - boundary = GeographicBoundary(sides_lons, sides_lats) + boundary = GeographicBoundary(self.area) + # Assert contour lons, lats = boundary.contour() assert np.allclose(lons, np.array([1., 1.5, 2., 3., 3.5, 4.])) assert np.allclose(lats, np.array([6., 6.5, 7., 8., 8.5, 9.])) def test_contour_closed(self): """Test that GeographicBoundary.contour(closed=True) returns the correct (lon,lat) tuple.""" - sides_lons = [np.array([1.0, 1.5, 2.0]), - np.array([2.0, 3.0]), - np.array([3.0, 3.5, 4.0]), - np.array([4.0, 1.0])] - sides_lats = [np.array([6.0, 6.5, 7.0]), - np.array([7.0, 8.0]), - np.array([8.0, 8.5, 9.0]), - np.array([9.0, 6.0])] # Define GeographicBoundary - boundary = GeographicBoundary(sides_lons, sides_lats) + boundary = GeographicBoundary(self.area) lons, lats = boundary.contour(closed=True) assert np.allclose(lons, np.array([1., 1.5, 2., 3., 3.5, 4., 1.])) assert np.allclose(lats, np.array([6., 6.5, 7., 8., 8.5, 9., 6.])) diff --git a/pyresample/test/test_boundary/test_order.py b/pyresample/test/test_boundary/test_order.py new file mode 100644 index 000000000..bc7815758 --- /dev/null +++ b/pyresample/test/test_boundary/test_order.py @@ -0,0 +1,133 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Nov 27 15:48:41 2023 + +@author: ghiggi +""" +import pytest +from pyresample.boundary.geographic_boundary import _is_clockwise_order + +class Test_Is_Clockwise_Order: + + def test_vertical_edges(self): + """Test with vertical polygon edges.""" + # North Hemisphere + first_point = (0, 10) + second_point = (0, 20) + point_inside = (1, 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + # South Hemisphere + first_point = (0, -20) + second_point = (0, -10) + point_inside = (1, - 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + @pytest.mark.parametrize("lon",[-180, -90, 0, 90, 180]) + def test_horizontal_edges(self, lon): + """Test with horizontal polygon edges.""" + # Point in northern hemisphere + first_point = (lon, 0) + second_point = (lon-10, 0) + point_inside = (1, 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + # Point in northern hemisphere + first_point = (lon, 0) + second_point = (lon+10, 0) + point_inside = (1, -15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + def test_diagonal_edges(self): + """Test with diagonal polygon edges.""" + point_inside = (20, 15) + + # Edge toward right (above point) --> clockwise + first_point = (0, 0) + second_point = (20, 20) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + # Edge toward right (below point) --> not clockwise + first_point = (0, 0) + second_point = (20, 10) + assert not _is_clockwise_order(first_point, second_point, point_inside) + + def test_polygon_edges_on_antimeridian(self): + """Test polygon edges touching the antimeridian edges.""" + ## Right side of antimeridian + # North Hemisphere + first_point = (-180, 10) + second_point = (-180, 20) + point_inside = (-179, 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + # South Hemisphere + first_point = (-180, -20) + second_point = (-180, -10) + point_inside = (-179, - 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + ## Left side of antimeridian + # North Hemisphere + first_point = (-180, 20) + second_point = (-180, 10) + point_inside = (179, 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + # South Hemisphere + first_point = (-180, -10) + second_point = (-180,-20) + point_inside = (179, - 15) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + @pytest.mark.parametrize("lon",[179, 180, -180, -179]) + def test_polygon_around_antimeridian(self, lon): + """Test polygon edges crossing antimeridian.""" + # North Hemisphere + first_point = (170, 10) + second_point = (-170, 10) + point_inside = (lon, 5) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + # South Hemisphere + first_point = (-170, -10) + second_point = (170, -10) + point_inside = (lon, - 5) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + @pytest.mark.parametrize("lon_pole",[-180, 90, 45, 0, 45, 90, 180]) + @pytest.mark.parametrize("lat",[85, 0, -85]) + def test_polygon_around_north_pole(self, lon_pole, lat): + """Test polygon edges around north pole (right to left).""" + point_inside = (lon_pole, 90) + first_point = (0, lat) + second_point = (-10, lat) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + @pytest.mark.parametrize("lon_pole",[-180, 90, 45, 0, 45, 90, 180]) + @pytest.mark.parametrize("lat",[85, 0, -85]) + def test_polygon_around_south_pole(self, lon_pole, lat): + """Test polygon edges around south pole (left to right).""" + point_inside = (lon_pole, -90) + first_point = (0, lat) + second_point = (10, lat) + assert _is_clockwise_order(first_point, second_point, point_inside) + assert not _is_clockwise_order(second_point, first_point , point_inside) + + + + + \ No newline at end of file diff --git a/pyresample/test/test_boundary/test_utils.py b/pyresample/test/test_boundary/test_utils.py new file mode 100644 index 000000000..024516e79 --- /dev/null +++ b/pyresample/test/test_boundary/test_utils.py @@ -0,0 +1,145 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# pyresample, Resampling of remote sensing image data in python +# +# Copyright (C) 2010-2022 Pyresample developers +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . +"""Test the boundary utility.""" +import numpy as np +import pytest +from pyresample.boundary.utils import ( + find_boundary_mask, + find_boundary_contour_indices, + get_ordered_contour, +) + + +@pytest.mark.parametrize("lonlat, expected", [ + # Case: All True values + ((np.array([[1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4]]), + np.array([[1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4]])), + np.array([[True, True, True, True], + [True, False, False, True], + [True, False, False, True], + [True, True, True, True]])), + + # Case: Multiple True values in the center + ((np.array([[np.inf, np.inf, np.inf, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, np.inf, np.inf, np.inf]]), + np.array([[np.inf, np.inf, np.inf, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, np.inf, np.inf, np.inf]])), + np.array([[False, False, False, False], + [False, True, True, False], + [False, True, True, False], + [False, False, False, False]])), +]) +def test_find_boundary_mask(lonlat, expected): + """Test boundary mask for lon lat array with non finite values.""" + lons, lats = lonlat + result = find_boundary_mask(lons, lats) + np.testing.assert_array_equal(result, expected, err_msg=f"Expected {expected}, but got {result}") + + +@pytest.mark.parametrize("boundary_mask, expected", [ + # Case: All True values + (np.array([[True, True, True, True], + [True, False, False, True], + [True, False, False, True], + [True, True, True, True]]), + np.array([[0, 0], [0, 1], [0, 2], [0, 3], + [1, 3], [2, 3], [3, 3], + [3, 2], [3, 1], [3, 0], + [2, 0], [1, 0]]) + ), + # Case: Multiple True values in the center + (np.array([[False, False, False, False], + [False, True, True, False], + [False, True, True, False], + [False, False, False, False]]), + np.array([[1, 1], [1, 2], [2, 2], [2, 1]]) + ), + # Case: Square on angle + (np.array([[True, True, True, False], + [True, False, True, False], + [True, True, True, False], + [False, False, False, False]]), + np.array([[0, 0], [0, 1], [0, 2], [1, 2], [2, 2], [2, 1], [2, 0], [1, 0]]) + ), + # Case: Cross Pattern + (np.array([[False, True, True, False], + [True, False, False, True], + [True, False, False, True], + [False, True, True, False]]), + np.array([[0, 1], [0, 2], [1, 3], [2, 3], [3, 2], [3, 1], [2, 0], [1, 0]]) + ), + # Case: Possibile infinit loop if not checking visited + (np.array([[1, 1, 1, 1, 0], + [1, 0, 1, 0, 0], + [1, 0, 1, 0, 0], + [1, 1, 0, 0, 0], + [1, 0, 0, 0, 0], + [1, 0, 0, 0, 0], + [0, 0, 0, 0, 0]]), + np.array([[0, 0], [0, 1], [0, 2], [0, 3], [1, 2], [2, 2], [3, 1], [3, 0], [2, 0], [1, 0]]) + ), +]) +def test_get_ordered_contour(boundary_mask, expected): + """Test order of the boundary contour indices (clockwise).""" + result = get_ordered_contour(boundary_mask) + np.testing.assert_array_equal(result, expected, err_msg=f"Expected {expected}, but got {result}") + + +@pytest.mark.parametrize("lonlat, expected", [ + # Case: All True values + ((np.array([[1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4]]), + np.array([[1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4], + [1, 2, 3, 4]])), + np.array([[0, 0], [0, 1], [0, 2], [0, 3], + [1, 3], [2, 3], [3, 3], + [3, 2], [3, 1], [3, 0], + [2, 0], [1, 0]]) + ), + # Case: Multiple True values in the center + ((np.array([[np.inf, np.inf, np.inf, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, np.inf, np.inf, np.inf]]), + np.array([[np.inf, np.inf, np.inf, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, 2, 3, np.inf], + [np.inf, np.inf, np.inf, np.inf]])), + np.array([[1, 1], [1, 2], [2, 2], [2, 1]]) + ), +]) +def test_find_boundary_contour_indices(lonlat, expected): + """Test order of the boundary contour indices (clockwise).""" + lons, lats = lonlat + result = find_boundary_contour_indices(lons, lats) + np.testing.assert_array_equal(result, expected, err_msg=f"Expected {expected}, but got {result}") \ No newline at end of file diff --git a/pyresample/test/test_geometry/test_area.py b/pyresample/test/test_geometry/test_area.py index a6cb41f74..3e2ae11ac 100644 --- a/pyresample/test/test_geometry/test_area.py +++ b/pyresample/test/test_geometry/test_area.py @@ -29,10 +29,9 @@ from pyresample import geo_filter, parse_area_file from pyresample.future.geometry import AreaDefinition, SwathDefinition from pyresample.future.geometry.area import ( - _get_geostationary_bounding_box_in_lonlats, + _get_geostationary_bounding_box, get_full_geostationary_bounding_box_in_proj_coords, get_geostationary_angle_extent, - get_geostationary_bounding_box_in_proj_coords, ignore_pyproj_proj_warnings, ) from pyresample.future.geometry.base import get_array_hashable @@ -262,7 +261,7 @@ def geos_mesoscale_area(): @pytest.fixture def truncated_geos_area(): - """Create a truncated geostationary area.""" + """Create a truncated geostationary area (SEVIRI above 30° lat).""" proj_dict = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'} area_extent = (5567248.0742, 5570248.4773, -5570248.4773, 1393687.2705) @@ -276,7 +275,7 @@ def truncated_geos_area(): @pytest.fixture def truncated_geos_area_in_space(): - """Create a truncated geostationary area.""" + """Create a geostationary area entirely out of the Earth disk !.""" proj_dict = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'} area_extent = (5575000, 5575000, 5570000, 5570000) @@ -1669,43 +1668,41 @@ def test_get_full_geostationary_bbox(self, truncated_geos_area): def test_get_geostationary_bbox_works_with_truncated_area(self, truncated_geos_area): """Ensure the geostationary bbox works when truncated.""" - lon, lat = _get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20) - expected_lon = np.array( - [-64.24072434653284, -68.69662326361153, -65.92516214783112, -60.726360278290336, - -47.39851775032484, 9.500000000000018, 66.39851775032487, 79.72636027829033, - 84.92516214783113, 87.69662326361151, 83.24072434653286]) - expected_lat = np.array( - [14.554922655532085, 17.768795771961937, 35.34328897185421, 52.597860701318254, 69.00533141646078, - 79.1481121862375, 69.00533141646076, 52.597860701318254, 35.34328897185421, 17.768795771961933, - 14.554922655532085]) + # NOTE: the results change if nb_points is changed + lon, lat = _get_geostationary_bounding_box(truncated_geos_area, coordinates="geographic", nb_points=5) + expected_lon = np.array([48.23903923, 27.09551769, 9.48612352, -8.12549639, + -29.28016639, -40.25837094, -47.39851775, 84.92516215, + 58.87273432]) + expected_lat = np.array([13.3415869, 12.89359038, 12.76978964, 12.89400717, 13.34272893, + 13.67772269, 69.00533142, 35.34328897, 13.66502884]) np.testing.assert_allclose(lon, expected_lon) np.testing.assert_allclose(lat, expected_lat) def test_get_geostationary_bbox_works_with_truncated_area_proj_coords(self, truncated_geos_area): """Ensure the geostationary bbox works when truncated.""" - x, y = get_geostationary_bounding_box_in_proj_coords(truncated_geos_area, 20) + # NOTE: the results change if nb_points is changed + x, y = _get_geostationary_bounding_box(truncated_geos_area, coordinates="projection", nb_points=5) - expected_x = np.array( - [-5209128.302753595, -5164828.965702432, -4393465.934674804, -3192039.8468840676, -1678154.6586309497, - 3.325297262895822e-10, 1678154.6586309501, 3192039.846884068, 4393465.934674805, 5164828.965702432, - 5209128.302753594]) - expected_y = np.array( - [1393687.2705, 1672427.7900638399, 3181146.6955466354, 4378472.798117005, 5147203.47659387, - 5412090.016106332, 5147203.476593869, 4378472.798117005, 3181146.695546635, 1672427.7900638392, - 1393687.2705]) + expected_x = np.array([3.71099865e+06, 1.85474922e+06, -1.50020155e+03, -1.85774963e+06, + -3.71399905e+06, -4.41458214e+06, -1.67815466e+06, 4.39346593e+06, + 4.39346593e+06]) + expected_y = np.array([1393687.2705, 1393687.2705, 1393687.2705, + 1393687.2705, 1393687.2705, 1393687.2705, 5147203.47659387, + 3181146.69554663, 1393687.2705]) np.testing.assert_allclose(x, expected_x) np.testing.assert_allclose(y, expected_y) def test_get_geostationary_bbox_does_not_contain_inf(self, truncated_geos_area): """Ensure the geostationary bbox does not contain np.inf.""" - lon, lat = _get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20) + lon, lat = _get_geostationary_bounding_box(truncated_geos_area, coordinates="geographic", nb_points=20) assert not any(np.isinf(lon)) assert not any(np.isinf(lat)) def test_get_geostationary_bbox_returns_empty_lonlats_in_space(self, truncated_geos_area_in_space): """Ensure the geostationary bbox is empty when in space.""" - lon, lat = _get_geostationary_bounding_box_in_lonlats(truncated_geos_area_in_space, 20) + lon, lat = _get_geostationary_bounding_box(truncated_geos_area_in_space, coordinates="geographic", + nb_points=20) assert len(lon) == 0 assert len(lat) == 0 @@ -1722,7 +1719,7 @@ def test_get_geostationary_bbox(self): geos_area.crs = CRS(proj_dict) geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.] - lon, lat = _get_geostationary_bounding_box_in_lonlats(geos_area, 20) + lon, lat = _get_geostationary_bounding_box(geos_area, coordinates="geographic", nb_points=20) expected_lon = np.array([-78.19662326, -75.42516215, -70.22636028, -56.89851775, 0., 56.89851775, 70.22636028, 75.42516215, 78.19662326, 79.23372832, 78.19662326, @@ -1747,7 +1744,7 @@ def test_get_geostationary_bbox(self): geos_area.crs = CRS(proj_dict) geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.] - lon, lat = _get_geostationary_bounding_box_in_lonlats(geos_area, 20) + lon, lat = _get_geostationary_bounding_box(geos_area, coordinates="geographic", nb_points=20) np.testing.assert_allclose(lon, expected_lon + lon_0) def test_get_geostationary_angle_extent(self): @@ -2113,7 +2110,7 @@ def test_get_geographic_sides_call_geostationary_utility(self, request, area_def def test_polar_south_pole_projection(self, south_pole_area): """Test boundary for polar projection around the South Pole.""" areadef = south_pole_area - boundary = areadef.geographic_boundary(order=None) + boundary = areadef.geographic_boundary() # Check boundary shape height, width = areadef.shape @@ -2131,7 +2128,7 @@ def test_north_pole_projection(self, north_pole_area): """Test boundary for polar projection around the North Pole.""" areadef = north_pole_area - boundary = areadef.geographic_boundary(order=None) + boundary = areadef.geographic_boundary() # Check boundary shape height, width = areadef.shape @@ -2151,13 +2148,13 @@ def test_full_disc_geostationary_projection(self, geos_fd_area): # Check default boundary shape default_n_vertices = 50 - boundary = areadef.geographic_boundary(vertices_per_side=None, order=None) + boundary = areadef.geographic_boundary(vertices_per_side=None, ) assert boundary.vertices.shape == (default_n_vertices, 2) # Check minimum boundary vertices n_vertices = 3 minimum_n_vertices = 4 - boundary = areadef.geographic_boundary(vertices_per_side=n_vertices, order=None) + boundary = areadef.geographic_boundary(vertices_per_side=n_vertices, ) assert boundary.vertices.shape == (minimum_n_vertices, 2) # Check odd number of vertices per side @@ -2168,7 +2165,7 @@ def test_full_disc_geostationary_projection(self, geos_fd_area): # Check boundary vertices n_vertices = 10 - boundary = areadef.geographic_boundary(vertices_per_side=n_vertices, order=None) + boundary = areadef.geographic_boundary(vertices_per_side=n_vertices, ) # Check boundary vertices is in correct order expected_vertices = np.array([[-7.54251621e+01, 3.53432890e+01], @@ -2186,7 +2183,7 @@ def test_full_disc_geostationary_projection(self, geos_fd_area): def test_global_platee_caree_projection(self, global_platee_caree_area): """Test boundary for global platee caree projection.""" areadef = global_platee_caree_area - boundary = areadef.geographic_boundary(order=None) + boundary = areadef.geographic_boundary() # Check boundary shape height, width = areadef.shape @@ -2211,7 +2208,7 @@ def test_global_platee_caree_projection(self, global_platee_caree_area): def test_minimal_global_platee_caree_projection(self, global_platee_caree_minimum_area): """Test boundary for global platee caree projection.""" areadef = global_platee_caree_minimum_area - boundary = areadef.geographic_boundary(order=None) + boundary = areadef.geographic_boundary() # Check boundary shape height, width = areadef.shape @@ -2228,7 +2225,7 @@ def test_minimal_global_platee_caree_projection(self, global_platee_caree_minimu def test_local_area_projection(self, local_meter_area): """Test local area projection in meter.""" areadef = local_meter_area - boundary = areadef.geographic_boundary(order=None) + boundary = areadef.geographic_boundary() # Check boundary shape height, width = areadef.shape diff --git a/pyresample/test/test_geometry/test_swath.py b/pyresample/test/test_geometry/test_swath.py index a16d3905f..36f09578c 100644 --- a/pyresample/test/test_geometry/test_swath.py +++ b/pyresample/test/test_geometry/test_swath.py @@ -604,7 +604,7 @@ def test_swath_definition(self, create_test_swath): # Define SwathDefinition and retrieve GeographicBoundary swath_def = create_test_swath(lons, lats) - boundary = swath_def.geographic_boundary(order=None) + boundary = swath_def.geographic_boundary() # Check boundary shape height, width = swath_def.shape diff --git a/pyresample/test/test_slicer.py b/pyresample/test/test_slicer.py index af082f896..db198296b 100644 --- a/pyresample/test/test_slicer.py +++ b/pyresample/test/test_slicer.py @@ -146,7 +146,7 @@ def test_barely_touching_chunks_intersection(self): assert x_slice.start > 0 and x_slice.stop < 100 assert y_slice.start > 0 and y_slice.stop >= 100 - def test_slicing_an_area_with_infinite_bounds(self): + def test_slicing_with_dst_area_with_infinite_edges(self): """Test slicing an area with infinite bounds.""" src_area = AreaDefinition('dst', 'dst area', None, {'ellps': 'WGS84', 'proj': 'merc'}, @@ -166,7 +166,14 @@ def test_slicing_an_area_with_infinite_bounds(self): slicer = create_slicer(src_area, dst_area) with pytest.raises(IncompatibleAreas): - slicer.get_slices() + slice_x, slice_y = slicer.get_slices() + + # Unreasonable slices if force_boundary_computations=True + # area_to_crop = src_area + # area_to_contain = dst_area + # slice_x, slice_y = slicer.get_slices() + # assert slice_x.start > 0 and slice_x.stop < 100 + # assert slice_y.start > 0 and slice_y.stop >= 100 def test_slicing_works_with_extents_of_different_units(self): """Test a problematic case."""