diff --git a/docs/source/howtos/spherical_geometry.rst b/docs/source/howtos/spherical_geometry.rst
index d2c0056ce..73fbdacf2 100644
--- a/docs/source/howtos/spherical_geometry.rst
+++ b/docs/source/howtos/spherical_geometry.rst
@@ -72,14 +72,13 @@ satellite passes. See trollschedule_ how to generate a list of satellite overpas
`area_def` is an :class:`~pyresample.geometry.AreaDefinition` object.
>>> from pyresample.spherical_utils import GetNonOverlapUnions
- >>> from pyresample.boundary import AreaDefBoundary
- >>> area_boundary = AreaDefBoundary(area_def, frequency=100) # doctest: +SKIP
- >>> area_boundary = area_boundary.contour_poly # doctest: +SKIP
+ >>> area_boundary = area_def.boundary(vertices_per_side=100) # 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.boundary().polygon) # doctest: +SKIP
>>> non_overlaps = GetNonOverlapUnions(list_of_polygons) # doctest: +SKIP
>>> non_overlaps.merge() # doctest: +SKIP
diff --git a/pyresample/__init__.py b/pyresample/__init__.py
index b229e0951..20aa92a79 100644
--- a/pyresample/__init__.py
+++ b/pyresample/__init__.py
@@ -56,8 +56,13 @@
from .version import get_versions # noqa
-__all__ = ['grid', 'image', 'kd_tree', 'utils', 'plot', 'geo_filter', 'geometry', 'CHUNK_SIZE',
- 'load_area', 'create_area_def', 'get_area_def', 'parse_area_file', 'convert_def_to_yaml']
+_root_path = os.path.dirname(os.path.realpath(__file__))
+
+__all__ = [
+ 'grid', 'image', 'kd_tree', 'utils', 'plot', 'geo_filter', 'geometry', 'CHUNK_SIZE',
+ 'load_area', 'create_area_def', 'get_area_def', 'parse_area_file', 'convert_def_to_yaml',
+ "_root_path",
+]
__version__ = get_versions()['version']
del get_versions
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/bilinear/_base.py b/pyresample/bilinear/_base.py
index 447428b1b..08d0b4718 100644
--- a/pyresample/bilinear/_base.py
+++ b/pyresample/bilinear/_base.py
@@ -599,14 +599,14 @@ def get_valid_indices_from_lonlat_boundaries(
target_geo_def, source_lons, source_lats, radius_of_influence):
"""Get valid indices from lonlat boundaries."""
# Resampling from swath to grid or from grid to grid
- lonlat_boundary = target_geo_def.get_boundary_lonlats()
-
- # Combine reduced and legal values
- return data_reduce.get_valid_index_from_lonlat_boundaries(
- lonlat_boundary[0],
- lonlat_boundary[1],
- source_lons, source_lats,
- radius_of_influence)
+ try:
+ sides_lons, sides_lats = target_geo_def.boundary().sides
+ valid_indices = data_reduce.get_valid_index_from_lonlat_boundaries(sides_lons, sides_lats,
+ source_lons, source_lats,
+ radius_of_influence)
+ except Exception:
+ valid_indices = np.ones(source_lons.size, dtype=bool)
+ return valid_indices
def get_slicer(data):
diff --git a/pyresample/boundary/__init__.py b/pyresample/boundary/__init__.py
new file mode 100644
index 000000000..8dec4345b
--- /dev/null
+++ b/pyresample/boundary/__init__.py
@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2014-2021 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 .
+"""The Boundary classes."""
+
+from pyresample.boundary.area_boundary import AreaBoundary, AreaDefBoundary
+from pyresample.boundary.planar_boundary import PlanarBoundary
+from pyresample.boundary.simple_boundary import SimpleBoundary
+from pyresample.boundary.spherical_boundary import SphericalBoundary
+
+__all__ = [
+ "SphericalBoundary",
+ "PlanarBoundary",
+ # Deprecated
+ "SimpleBoundary",
+ "AreaBoundary",
+ "AreaDefBoundary",
+]
diff --git a/pyresample/boundary.py b/pyresample/boundary/area_boundary.py
similarity index 86%
rename from pyresample/boundary.py
rename to pyresample/boundary/area_boundary.py
index 87f9b2ae1..b3087b829 100644
--- a/pyresample/boundary.py
+++ b/pyresample/boundary/area_boundary.py
@@ -1,7 +1,7 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
-# Copyright (c) 2014-2021 Pyresample developers
+# 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
@@ -15,8 +15,7 @@
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
-"""The Boundary classes."""
-
+"""Deprecated Boundary, AreaBoundary, AreaDefBoundary class."""
import logging
import warnings
@@ -44,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))
@@ -62,6 +64,9 @@ class AreaBoundary(Boundary):
def __init__(self, *sides):
Boundary.__init__(self)
+ warnings.warn("'AreaBoundary' will be removed in the future. " +
+ "Use the Swath/AreaDefinition 'boundary' method instead!.",
+ PendingDeprecationWarning, stacklevel=2)
# Check 4 sides are provided
if len(sides) != 4:
raise ValueError("AreaBoundary expects 4 sides.")
@@ -117,28 +122,13 @@ def decimate(self, ratio):
class AreaDefBoundary(AreaBoundary):
- """Boundaries for area definitions (pyresample)."""
+ """Boundaries for a pyresample AreaDefinition."""
def __init__(self, area, frequency=1):
- lons, lats = area.get_bbox_lonlats()
+ sides_lons, sides_lats = area.boundary().sides
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))
-
+ AreaBoundary.__init__(self, *zip(sides_lons, sides_lats))
if frequency != 1:
self.decimate(frequency)
-
-
-class SimpleBoundary(object):
- """Container for geometry boundary.
-
- Labelling starts in upper left corner and proceeds clockwise
- """
-
- def __init__(self, side1, side2, side3, side4):
- self.side1 = side1
- self.side2 = side2
- self.side3 = side3
- self.side4 = side4
diff --git a/pyresample/boundary/base_boundary.py b/pyresample/boundary/base_boundary.py
new file mode 100644
index 000000000..7306237fc
--- /dev/null
+++ b/pyresample/boundary/base_boundary.py
@@ -0,0 +1,127 @@
+#!/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 .
+"""Define the BaseBoundary class."""
+
+import logging
+
+import numpy as np
+
+from pyresample.boundary.sides import BoundarySides
+
+logger = logging.getLogger(__name__)
+
+
+class BaseBoundary:
+ """Base class for boundary objects."""
+ __slots__ = ["_sides_x", "_sides_y"]
+
+ 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, area)
+ self.is_counterclockwise = not self.is_clockwise
+ self._set_order(order=None) # FIX !
+
+ 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:
+ self._actual_order = "clockwise"
+ else:
+ self._actual_order = "counterclockwise"
+
+ if order is None:
+ self._wished_order = self._actual_order
+ else:
+ if order not in ["clockwise", "counterclockwise"]:
+ raise ValueError("Valid 'order' is 'clockwise' or 'counterclockwise'")
+ self._wished_order = order
+
+ def set_clockwise(self):
+ """Set clockwise order for vertices retrieval."""
+ self._wished_order = "clockwise"
+ return self
+
+ def set_counterclockwise(self):
+ """Set counterclockwise order for vertices retrieval."""
+ self._wished_order = "counterclockwise"
+ return self
+
+ @property
+ def sides(self):
+ """Return the boundary sides as a tuple of (sides_x, sides_y) arrays."""
+ return self._sides_x, self._sides_y
+
+ @property
+ def _x(self):
+ """Retrieve boundary x vertices."""
+ xs = self._sides_x.vertices
+ if self._wished_order == self._actual_order:
+ return xs
+ else:
+ return xs[::-1]
+
+ @property
+ def _y(self):
+ """Retrieve boundary y vertices."""
+ ys = self._sides_y.vertices
+ if self._wished_order == self._actual_order:
+ return ys
+ else:
+ return ys[::-1]
+
+ @property
+ def vertices(self):
+ """Return boundary vertices 2D array [x, y]."""
+ vertices = np.vstack((self._x, self._y)).T
+ vertices = vertices.astype(np.float64, copy=False) # Important for spherical ops.
+ return vertices
+
+ def contour(self, closed=False):
+ """Return the (x, y) tuple of the boundary object.
+
+ If excludes the last element of each side because it's included in the next side.
+ If closed=False (the default), the last vertex is not equal to the first vertex
+ If closed=True, the last vertex is set to be equal to the first
+ closed=True is required for shapely Polygon creation.
+ closed=False is required for pyresample SPolygon creation.
+ """
+ x = self._x
+ y = self._y
+ if closed:
+ x = np.hstack((x, x[0]))
+ y = np.hstack((y, y[0]))
+ return x, y
+
+ def to_shapely_polygon(self):
+ """Define a Shapely Polygon."""
+ from shapely.geometry import Polygon
+ self = self.set_counterclockwise()
+ x, y = self.contour(closed=True)
+ return Polygon(zip(x, y))
diff --git a/pyresample/boundary/planar_boundary.py b/pyresample/boundary/planar_boundary.py
new file mode 100644
index 000000000..7b10074f5
--- /dev/null
+++ b/pyresample/boundary/planar_boundary.py
@@ -0,0 +1,94 @@
+#!/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 .
+"""Define the PlanarBoundary class."""
+
+import logging
+
+import numpy as np
+
+from pyresample.boundary.base_boundary import BaseBoundary
+
+logger = logging.getLogger(__name__)
+
+
+def _is_projection_boundary_clockwise(sides_x, sides_y):
+ """Determine if the boundary is clockwise-defined in planar coordinates."""
+ from shapely.geometry import Polygon
+ x = np.concatenate([xs[:-1] for xs in sides_x])
+ y = np.concatenate([ys[:-1] for ys in sides_y])
+ x = np.hstack((x, x[0]))
+ y = np.hstack((y, y[0]))
+ polygon = Polygon(zip(x, y))
+ return not polygon.exterior.is_ccw
+
+
+class PlanarBoundary(BaseBoundary):
+ """Projection Boundary object.
+
+ The inputs must be the x and y sides of the projection.
+ It expects the projection coordinates to be planar (i.e. metric, radians).
+ """
+
+ @classmethod
+ def _check_is_boundary_clockwise(cls, sides_x, sides_y, area=None):
+ """SphericalBoundary specific implementation."""
+ return _is_projection_boundary_clockwise(sides_x=sides_x, sides_y=sides_y)
+
+ @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 = 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):
+ """Retrieve boundary x vertices."""
+ return self._x
+
+ @property
+ def y(self):
+ """Retrieve boundary y vertices."""
+ return self._y
+
+ 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
+
+ if self.cartopy_crs is None and crs is None:
+ 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.to_shapely_polygon()
+ p = plot_geometries(geometries=[geom], crs=crs,
+ ax=ax, subplot_kw=subplot_kw, alpha=alpha, **kwargs)
+ return p
diff --git a/pyresample/boundary/sides.py b/pyresample/boundary/sides.py
new file mode 100644
index 000000000..e52499ea6
--- /dev/null
+++ b/pyresample/boundary/sides.py
@@ -0,0 +1,86 @@
+#!/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 .
+"""Define the BoundarySides class."""
+
+import logging
+
+import numpy as np
+
+logger = logging.getLogger(__name__)
+
+
+class BoundarySides:
+ """A class to represent the sides of an area boundary.
+
+ The sides are stored as a tuple of 4 numpy arrays, each representing the
+ coordinate (geographic or projected) of the vertices of the boundary side.
+ The sides must be stored in the order (top, right, left, bottom),
+ which refers to the side position with respect to the coordinate array.
+ The first row of the coordinate array correspond to the top side, the last row to the bottom side,
+ the first column to the left side and the last column to the right side.
+ Please note that the last vertex of each side must be equal to the first vertex of the next side.
+ """
+ __slots__ = ['_sides']
+
+ def __init__(self, sides):
+ """Initialize the BoundarySides object."""
+ if len(sides) != 4 or not all(isinstance(side, np.ndarray) and side.ndim == 1 for side in sides):
+ raise ValueError("Sides must be a list of four numpy arrays.")
+
+ if not all(np.array_equal(sides[i][-1], sides[(i + 1) % 4][0]) for i in range(4)):
+ raise ValueError("The last element of each side must be equal to the first element of the next side.")
+
+ self._sides = tuple(sides) # Store as a tuple
+
+ @property
+ def top(self):
+ """Return the vertices of the top side."""
+ return self._sides[0]
+
+ @property
+ def right(self):
+ """Return the vertices of the right side."""
+ return self._sides[1]
+
+ @property
+ def bottom(self):
+ """Return the vertices of the bottom side."""
+ return self._sides[2]
+
+ @property
+ def left(self):
+ """Return the vertices of the left side."""
+ return self._sides[3]
+
+ @property
+ def vertices(self):
+ """Return the vertices of the concatenated sides.
+
+ Note that the last element of each side is discarded to avoid duplicates.
+ """
+ return np.concatenate([side[:-1] for side in self._sides])
+
+ def __iter__(self):
+ """Return an iterator over the sides."""
+ return iter(self._sides)
+
+ def __getitem__(self, index):
+ """Return the side at the given index."""
+ if not isinstance(index, int) or not 0 <= index < 4:
+ raise IndexError("Index must be an integer from 0 to 3.")
+ return self._sides[index]
diff --git a/pyresample/boundary/simple_boundary.py b/pyresample/boundary/simple_boundary.py
new file mode 100644
index 000000000..ff66302c5
--- /dev/null
+++ b/pyresample/boundary/simple_boundary.py
@@ -0,0 +1,35 @@
+#!/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 .
+"""The deprecated SimpleBoundary class."""
+
+import logging
+
+logger = logging.getLogger(__name__)
+
+
+class SimpleBoundary(object):
+ """Container for geometry boundary.
+
+ Labelling starts in upper left corner and proceeds clockwise
+ """
+
+ def __init__(self, side1, side2, side3, side4):
+ self.side1 = side1
+ self.side2 = side2
+ self.side3 = side3
+ self.side4 = side4
diff --git a/pyresample/boundary/spherical_boundary.py b/pyresample/boundary/spherical_boundary.py
new file mode 100644
index 000000000..48a852930
--- /dev/null
+++ b/pyresample/boundary/spherical_boundary.py
@@ -0,0 +1,177 @@
+#!/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 .
+"""Define the SphericalBoundary class."""
+
+import logging
+
+import numpy as np
+from pyproj import CRS
+
+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 _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 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 !
+ """
+ from pyresample.spherical import Arc, SCoordinate
+ 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 _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 SphericalBoundary(BaseBoundary, OldBoundary):
+ """SphericalBoundary object.
+
+ The inputs must be the list of longitude and latitude boundary sides.
+ """
+ # NOTES:
+ # - Boundary provide the ancient method contour_poly and draw
+ # from the old interface for compatibility to AreaBoundary
+
+ @classmethod
+ def _check_is_boundary_clockwise(cls, sides_x, sides_y, area):
+ """SphericalBoundary specific implementation."""
+ 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, 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
+
+ # 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):
+ """Retrieve boundary longitude vertices."""
+ return self._x
+
+ @property
+ def lats(self):
+ """Retrieve boundary latitude vertices."""
+ return self._y
+
+ @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.to_shapely_polygon()
+ crs = ccrs.Geodetic()
+ p = plot_geometries(geometries=[geom], crs=crs,
+ 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..b884b2b18
--- /dev/null
+++ b/pyresample/boundary/utils.py
@@ -0,0 +1,85 @@
+#!/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
diff --git a/pyresample/data_reduce.py b/pyresample/data_reduce.py
index 5e4b0518d..154f0e2dc 100644
--- a/pyresample/data_reduce.py
+++ b/pyresample/data_reduce.py
@@ -63,7 +63,7 @@ def get_valid_index_from_cartesian_grid(cart_grid, lons, lats,
Parameters
----------
- chart_grid : numpy array
+ cart_grid : numpy array
Grid of area cartesian coordinates
lons : numpy array
Swath lons
@@ -103,7 +103,8 @@ def _get_lats(z):
return valid_index
-def swath_from_lonlat_grid(grid_lons, grid_lats, lons, lats, data,
+def swath_from_lonlat_grid(grid_lons, grid_lats,
+ lons, lats, data,
radius_of_influence):
"""Make coarse data reduction of swath data by comparison with lon lat grid.
@@ -137,20 +138,21 @@ def swath_from_lonlat_grid(grid_lons, grid_lats, lons, lats, data,
return lons, lats, data
-def swath_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, data,
+def swath_from_lonlat_boundaries(boundary_lons, boundary_lats,
+ lons, lats, data,
radius_of_influence):
"""Make coarse data reduction of swath data by comparison with lon lat boundary.
Parameters
----------
boundary_lons : numpy array
- Grid of area lons
+ Longitude BoundarySide object
boundary_lats : numpy array
- Grid of area lats
+ Latitude BoundarySide object
lons : numpy array
- Swath lons
+ Swath longitude
lats : numpy array
- Swath lats
+ Swath latitude
data : numpy array
Swath data
radius_of_influence : float
@@ -162,7 +164,9 @@ def swath_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, data,
Reduced swath data and coordinate set
"""
valid_index = get_valid_index_from_lonlat_boundaries(boundary_lons,
- boundary_lats, lons, lats, radius_of_influence)
+ boundary_lats,
+ lons, lats,
+ radius_of_influence)
lons = lons[valid_index]
lats = lats[valid_index]
@@ -212,10 +216,10 @@ def get_valid_index_from_lonlat_grid(grid_lons, grid_lats, lons, lats, radius_of
def get_valid_index_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, radius_of_influence):
"""Find relevant indices from grid boundaries using the winding number theorem."""
- valid_index = _get_valid_index(boundary_lons.side1, boundary_lons.side2,
- boundary_lons.side3, boundary_lons.side4,
- boundary_lats.side1, boundary_lats.side2,
- boundary_lats.side3, boundary_lats.side4,
+ valid_index = _get_valid_index(boundary_lons.top, boundary_lons.right,
+ boundary_lons.bottom, boundary_lons.left,
+ boundary_lats.top, boundary_lats.right,
+ boundary_lats.bottom, boundary_lats.left,
lons, lats, radius_of_influence)
return valid_index
diff --git a/pyresample/future/geometry/_subset.py b/pyresample/future/geometry/_subset.py
index bb04cf69b..4639f9849 100644
--- a/pyresample/future/geometry/_subset.py
+++ b/pyresample/future/geometry/_subset.py
@@ -1,6 +1,7 @@
"""Functions and tools for subsetting a geometry object."""
from __future__ import annotations
+import logging
import math
from typing import TYPE_CHECKING, Any
@@ -10,13 +11,14 @@
# must be imported inside functions in the geometry modules if needed
# to avoid circular dependencies
from pyresample._caching import cache_to_json_if
-from pyresample.boundary import Boundary
-from pyresample.geometry import get_geostationary_bounding_box_in_lonlats, logger
+from pyresample.boundary import SphericalBoundary
from pyresample.utils import check_slice_orientation
if TYPE_CHECKING:
from pyresample import AreaDefinition
+logger = logging.getLogger(__name__)
+
@cache_to_json_if("cache_geometry_slices")
def get_area_slices(
@@ -47,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.")
@@ -95,12 +97,13 @@ def _get_slice_starts_stops(src_area, area_to_cover):
return xstart, xstop, ystart, ystop
-def _get_area_boundary(area_to_cover: AreaDefinition) -> Boundary:
+def _get_area_boundary(area_to_cover: AreaDefinition) -> SphericalBoundary:
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)
+ vertices_per_side = None
+ else:
+ vertices_per_side = max(max(*area_to_cover.shape) // 100 + 1, 3)
+ return area_to_cover.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 d0ff7c85b..cc4e1e42b 100644
--- a/pyresample/future/geometry/area.py
+++ b/pyresample/future/geometry/area.py
@@ -26,9 +26,9 @@
from pyresample.geometry import AreaDefinition as LegacyAreaDefinition # noqa
from pyresample.geometry import ( # noqa
DynamicAreaDefinition,
+ _get_geostationary_bounding_box,
get_full_geostationary_bounding_box_in_proj_coords,
get_geostationary_angle_extent,
- get_geostationary_bounding_box_in_lonlats,
get_geostationary_bounding_box_in_proj_coords,
ignore_pyproj_proj_warnings,
)
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 2c2b52254..0ea3bfb10 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,6 +108,7 @@ def __init__(self, lons=None, lats=None, nprocs=1):
self.ndim = None
self.cartesian_coords = None
self.hash = None
+ self._boundary_contour_idx = None
def __getitem__(self, key):
"""Slice a 2D geographic definition."""
@@ -250,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):
@@ -273,8 +275,15 @@ def get_lonlats_dask(self, chunks=None):
chunks = CHUNK_SIZE # FUTURE: Use a global config object instead
return self.get_lonlats(chunks=chunks)
+ def get_proj_coords(self, data_slice=None, chunks=None, **kwargs):
+ """Return AreaDefinition projection coordinates."""
+ class_name = self.__class__.__name__
+ raise NotImplementedError(f"get_proj_coords is not available for geometry {class_name}.")
+
def get_boundary_lonlats(self):
"""Return Boundary objects."""
+ warnings.warn("'get_boundary_lonlats' is deprecated. Please use "
+ "'area.boundary().sides'.", DeprecationWarning, stacklevel=2)
s1_lon, s1_lat = self.get_lonlats(data_slice=(0, slice(None)))
s2_lon, s2_lat = self.get_lonlats(data_slice=(slice(None), -1))
s3_lon, s3_lat = self.get_lonlats(data_slice=(-1, slice(None, None, -1)))
@@ -282,14 +291,21 @@ def get_boundary_lonlats(self):
return (SimpleBoundary(s1_lon.squeeze(), s2_lon.squeeze(), s3_lon.squeeze(), s4_lon.squeeze()),
SimpleBoundary(s1_lat.squeeze(), s2_lat.squeeze(), s3_lat.squeeze(), s4_lat.squeeze()))
- def get_bbox_lonlats(self, vertices_per_side: Optional[int] = None, force_clockwise: bool = True,
+ @property
+ def is_geostationary(self):
+ """Whether this geometry is in a geostationary satellite projection or not."""
+ return False
+
+ def get_bbox_lonlats(self, vertices_per_side: Optional[int] = None,
+ force_clockwise: bool = True,
frequency: Optional[int] = None) -> tuple:
- """Return the bounding box lons and lats.
+ """Return the bounding box lons and lats sides.
Args:
vertices_per_side:
The number of points to provide for each side. By default (None)
- the full width and height will be provided.
+ the full width and height will be provided, except for geostationary
+ projections where by default only 50 points are selected.
frequency:
Deprecated, use vertices_per_side
force_clockwise:
@@ -320,45 +336,174 @@ def get_bbox_lonlats(self, vertices_per_side: Optional[int] = None, force_clockw
if frequency is not None:
warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
PendingDeprecationWarning, stacklevel=2)
+
vertices_per_side = vertices_per_side or frequency
- lons, lats = self._get_bbox_elements(self.get_lonlats, vertices_per_side)
+ sides_lons, sides_lats = self._get_geographic_sides(vertices_per_side=vertices_per_side)
+ warnings.warn("`get_bbox_lonlats` is pending deprecation. Use `area.boundary().sides` instead",
+ PendingDeprecationWarning, stacklevel=2)
+ # NOTE: This check is erroneous:
+ # - For SwathDefinition must be checked with respect to a local internal point.
+ # - For AreaDefinition, the order already defines the enclosing area !
+ # - For many AreaDefinition, this is failing (infinite edges ...)
if force_clockwise and not self._corner_is_clockwise(
- lons[0][-2], lats[0][-2], lons[0][-1], lats[0][-1], lons[1][1], lats[1][1]):
+ sides_lons[0][-2], sides_lats[0][-2],
+ sides_lons[0][-1], sides_lats[0][-1],
+ sides_lons[1][1], sides_lats[1][1]):
# going counter-clockwise
- lons, lats = self._reverse_boundaries(lons, lats)
- return lons, lats
+ sides_lons, sides_lats = self._reverse_boundaries(sides_lons, sides_lats)
+ return sides_lons, sides_lats
+
+ def _get_dummy_sides(self, arr, vertices_per_side):
+ """Retrieve a 'dummy' boundary side list for AreaDefinition with boundaries out of the Earth disk.
+
+ The second and fourth sides are always of length 2.
+ """
+ N = len(arr)
+ if vertices_per_side is None:
+ vertices_per_side = int(N / 2)
+ top_side = arr[slice(0, vertices_per_side)]
+ bottom_side = arr[slice(vertices_per_side, N)]
+ else:
+ # Sample vertices
+ # - Split the array into two halves
+ first_half = arr[:N // 2]
+ second_half = arr[N // 2:]
+ # - Adjust the number of vertices if necessary
+ vertices_per_side = min(vertices_per_side, len(first_half), len(second_half))
+ # - Sample points from each half
+ top_side_indices = np.round(np.linspace(0, len(first_half) - 1, vertices_per_side)).astype(int)
+ bottom_side_indices = np.round(np.linspace(0, len(second_half) - 1, vertices_per_side)).astype(int)
+ top_side = first_half[top_side_indices]
+ bottom_side = second_half[bottom_side_indices]
+ # - Create side object
+ sides = [top_side, np.array([top_side[-1], bottom_side[0]]),
+ bottom_side, np.array([bottom_side[-1], top_side[0]])]
+ return sides
+
+ def _get_geostationary_boundary_sides(self, vertices_per_side=None, coordinates="geographic"):
+ """Retrieve the boundary sides list for geostationary projections with out-of-Earth disk coordinates.
+
+ If vertices_per_side is too small, there is the risk to loose boundary side points
+ at the intersection corners between the CRS bounds polygon and the area
+ extent polygon (which could exclude relevant regions of the geos area).
+
+ The boundary sides right (1) and side left (3) are set to length 2.
+ """
+ # TODO:
+ # - Evaluate nb_points required for FULL DISC and CONUS area !
+ # - Maybe just use vertices_per_side for _get_geostationary_bounding_box and
+ # then return all the returned vertices --> _get_geostationary_dummy_sides(x, None)
+
+ # Define default frequency
+ if vertices_per_side is None:
+ vertices_per_side = 50
+ # Ensure at least 4 points are used
+ if vertices_per_side < 4:
+ vertices_per_side = 4
+ # Ensure an even number of vertices for side creation
+ if (vertices_per_side % 2) != 0:
+ vertices_per_side = vertices_per_side + 1
+ # Retrieve coordinates (x,y) or (lon, lat)
+ x, y = _get_geostationary_bounding_box(self, coordinates=coordinates, nb_points=vertices_per_side)
+ # Ensure that a portion of the area is within the Earth disk.
+ if x.shape[0] < 4:
+ raise ValueError("The geostationary projection area is entirely out of the Earth disk.")
+ # Retrieve dummy sides for GEO
+ # - _get_geostationary_bounding_box does not guarantee to return nb_points and even points!
+ 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
+
+ def _compute_boundary_mask(self):
+ """Compute valid boundary mask for AreaDefinition(s) with sides out of the Earth disk."""
+ 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_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.
+
+ Args:
+ vertices_per_side:
+ The number of points to provide for each side.
+ 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.
- def _get_bbox_elements(self, coord_fun, vertices_per_side: Optional[int] = None) -> tuple:
- s1_slice, s2_slice, s3_slice, s4_slice = self._get_bbox_slices(vertices_per_side)
- s1_dim1, s1_dim2 = coord_fun(data_slice=s1_slice)
- s2_dim1, s2_dim2 = coord_fun(data_slice=s2_slice)
- s3_dim1, s3_dim2 = coord_fun(data_slice=s3_slice)
- s4_dim1, s4_dim2 = coord_fun(data_slice=s4_slice)
- dim1, dim2 = zip(*[(s1_dim1.squeeze(), s1_dim2.squeeze()),
- (s2_dim1.squeeze(), s2_dim2.squeeze()),
- (s3_dim1.squeeze(), s3_dim2.squeeze()),
- (s4_dim1.squeeze(), s4_dim2.squeeze())])
- if hasattr(dim1[0], 'compute') and da is not None:
- dim1, dim2 = da.compute(dim1, dim2)
- clean_dim1, clean_dim2 = self._filter_bbox_nans(dim1, dim2)
- return clean_dim1, clean_dim2
-
- def _filter_bbox_nans(
+ Returns:
+ The output structure is a tuple of two lists of four elements each.
+ The first list contains the longitude coordinates.
+ The second list contains the latitude coordinates.
+ Each list element is a numpy array representing a specific side of the geometry.
+ The order of the sides are [top", "right", "bottom", "left"]
+ """
+ is_swath = self.__class__.__name__ == "SwathDefinition"
+ if not is_swath and _is_any_corner_out_of_earth_disk(self):
+ # Geostationary
+ if self.is_geostationary:
+ return self._get_geostationary_boundary_sides(vertices_per_side=vertices_per_side,
+ coordinates="geographic")
+ # Polar Projections, Global Planar Projections (Mollweide, Robinson)
+ # - Retrieve dummy right and left sides
+ 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):
+ """Return the boundary sides."""
+ top_slice, right_slice, bottom_slice, left_slice = self._get_bbox_slices(vertices_per_side)
+ top_dim1, top_dim2 = coord_fun(data_slice=top_slice)
+ right_dim1, right_dim2 = coord_fun(data_slice=right_slice)
+ bottom_dim1, bottom_dim2 = coord_fun(data_slice=bottom_slice)
+ left_dim1, left_dim2 = coord_fun(data_slice=left_slice)
+ sides_dim1, sides_dim2 = zip(*[(top_dim1.squeeze(), top_dim2.squeeze()),
+ (right_dim1.squeeze(), right_dim2.squeeze()),
+ (bottom_dim1.squeeze(), bottom_dim2.squeeze()),
+ (left_dim1.squeeze(), left_dim2.squeeze())])
+ if hasattr(sides_dim1[0], 'compute') and da is not None:
+ sides_dim1, sides_dim2 = da.compute(sides_dim1, sides_dim2)
+ sides_dim1, sides_dim2 = self._filter_sides_nans(sides_dim1, sides_dim2)
+ return sides_dim1, sides_dim2
+
+ def _filter_sides_nans(
self,
dim1_sides: list[np.ndarray],
dim2_sides: list[np.ndarray],
) -> tuple[list[np.ndarray], list[np.ndarray]]:
+ """Remove nan and inf values present in each side."""
new_dim1_sides = []
new_dim2_sides = []
for dim1_side, dim2_side in zip(dim1_sides, dim2_sides):
- is_valid_mask = ~(np.isnan(dim1_side) | np.isnan(dim2_side))
+ is_valid_mask = ~(~np.isfinite(dim1_side) | ~np.isfinite(dim1_side))
if not is_valid_mask.any():
- raise ValueError("Can't compute swath bounding coordinates. At least one side is completely invalid.")
+ raise ValueError("Can't compute boundary coordinates. At least one side is completely invalid.")
new_dim1_sides.append(dim1_side[is_valid_mask])
new_dim2_sides.append(dim2_side[is_valid_mask])
return new_dim1_sides, new_dim2_sides
def _get_bbox_slices(self, vertices_per_side):
+ """Return the slices for the bounding box of the current area.
+
+ Args:
+ vertices_per_side:
+ The number of points to provide for each side.
+ Results:
+ The output structure is a tuple of four slices, one for each side.
+ The order of the sides are [top", "right", "bottom", "left"]
+ """
+ if len(self.shape) == 1:
+ raise ValueError("The area must have 2 dimensions to retrieve the boundary sides.")
height, width = self.shape
if vertices_per_side is None:
row_num = height
@@ -366,11 +511,11 @@ def _get_bbox_slices(self, vertices_per_side):
else:
row_num = vertices_per_side
col_num = vertices_per_side
- s1_slice = (0, np.linspace(0, width - 1, col_num, dtype=int))
- s2_slice = (np.linspace(0, height - 1, row_num, dtype=int), -1)
- s3_slice = (-1, np.linspace(width - 1, 0, col_num, dtype=int))
- s4_slice = (np.linspace(height - 1, 0, row_num, dtype=int), 0)
- return s1_slice, s2_slice, s3_slice, s4_slice
+ top_slice = (0, np.linspace(0, width - 1, col_num, dtype=int))
+ right_slice = (np.linspace(0, height - 1, row_num, dtype=int), -1)
+ bottom_slice = (-1, np.linspace(width - 1, 0, col_num, dtype=int))
+ left_slice = (np.linspace(height - 1, 0, row_num, dtype=int), 0)
+ return top_slice, right_slice, bottom_slice, left_slice
@staticmethod
def _reverse_boundaries(sides_lons: list, sides_lats: list) -> tuple:
@@ -411,24 +556,29 @@ def _corner_is_clockwise(lon1, lat1, corner_lon, corner_lat, lon2, lat2):
def get_edge_lonlats(self, vertices_per_side=None, frequency=None):
"""Get the concatenated boundary of the current swath."""
if frequency is not None:
- warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
+ warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead.",
PendingDeprecationWarning, stacklevel=2)
+ msg = "`get_edge_lonlats` is pending deprecation"
+ msg += "Use `area.boundary(vertices_per_side=vertices_per_side).contour()` instead."
+ warnings.warn(msg, PendingDeprecationWarning, stacklevel=2)
vertices_per_side = vertices_per_side or frequency
- lons, lats = self.get_bbox_lonlats(vertices_per_side=vertices_per_side, force_clockwise=False)
- blons = np.ma.concatenate(lons)
- blats = np.ma.concatenate(lats)
- return blons, blats
+ lons, lats = self.boundary(vertices_per_side=vertices_per_side).contour()
+ return lons, lats
- def boundary(self, vertices_per_side=None, force_clockwise=False, frequency=None):
+ def boundary(self, *, vertices_per_side=None, force_clockwise=None, frequency=None):
"""Retrieve the AreaBoundary object.
Parameters
----------
vertices_per_side:
- (formerly `frequency`) The number of points to provide for each side. By default (None)
- the full width and height will be provided.
+ (formerly `frequency`) The number of points to provide for each side.
+ 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.
force_clockwise:
- Perform minimal checks and reordering of coordinates to ensure
+ DEPRECATED. IS NOT USED ANYMORE !
+ 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
@@ -436,14 +586,16 @@ def boundary(self, vertices_per_side=None, force_clockwise=False, frequency=None
operations assume that coordinates are clockwise.
Default is False.
"""
- from pyresample.boundary import AreaBoundary
+ from pyresample.boundary import SphericalBoundary
if frequency is not None:
warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
PendingDeprecationWarning, stacklevel=2)
+ if force_clockwise is not None:
+ warnings.warn("The `force_clockwise` argument is not used anymore. "
+ "Please remove the argument from the boundary() call !!!",
+ PendingDeprecationWarning, stacklevel=2)
vertices_per_side = vertices_per_side or frequency
- lon_sides, lat_sides = self.get_bbox_lonlats(vertices_per_side=vertices_per_side,
- force_clockwise=force_clockwise)
- return AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
+ return SphericalBoundary(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.
@@ -500,7 +652,7 @@ def get_cartesian_coords(self, nprocs=None, data_slice=None, cache=False):
@property
def corners(self):
- """Return the corners of the current area."""
+ """Return the corners (pixel centroids) of the current area."""
from pyresample.spherical_geometry import Coordinate
return [Coordinate(*self.get_lonlat(0, 0)),
Coordinate(*self.get_lonlat(0, -1)),
@@ -983,7 +1135,7 @@ def compute_optimal_bb_area(self, proj_dict=None, resolution=None):
proj_dict = self.compute_bb_proj_params(proj_dict)
area = DynamicAreaDefinition(area_id, description, proj_dict)
- lons, lats = self.get_edge_lonlats()
+ lons, lats = self.boundary(vertices_per_side=None).contour()
return area.freeze((lons, lats), shape=(height, width))
@@ -997,7 +1149,7 @@ class DynamicAreaDefinition(object):
Note that if the provided projection is geographic (lon/lat degrees) and
the provided longitude and latitude data crosses the anti-meridian
(-180/180), the resulting area will be the smallest possible in order to
- contain that data and avoid a large area spanning from -180 to 180
+ contain that data and boundary spanning from -180 to 180
longitude. This means the resulting AreaDefinition will have a right-most
X extent greater than 180 degrees. This does not apply to data crossing
the north or south pole as there is no "smallest" area in this case.
@@ -1563,64 +1715,63 @@ def is_geostationary(self):
return False
return 'geostationary' in coord_operation.method_name.lower()
- def _get_geo_boundary_sides(self, vertices_per_side=None):
- """Retrieve the boundary sides list for geostationary projections."""
- # Define default frequency
- if vertices_per_side is None:
- vertices_per_side = 50
- # Ensure at least 4 points are used
- if vertices_per_side < 4:
- vertices_per_side = 4
- # Ensure an even number of vertices for side creation
- if (vertices_per_side % 2) != 0:
- vertices_per_side = vertices_per_side + 1
- lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=vertices_per_side)
- # Retrieve dummy sides for GEO (side1 and side3 always of length 2)
- side02_step = int(vertices_per_side / 2) - 1
- lon_sides = [lons[slice(0, side02_step + 1)],
- lons[slice(side02_step, side02_step + 1 + 1)],
- lons[slice(side02_step + 1, side02_step * 2 + 1 + 1)],
- np.append(lons[side02_step * 2 + 1], lons[0])
- ]
- lat_sides = [lats[slice(0, side02_step + 1)],
- lats[slice(side02_step, side02_step + 1 + 1)],
- lats[slice(side02_step + 1, side02_step * 2 + 1 + 1)],
- np.append(lats[side02_step * 2 + 1], lats[0])
- ]
- return lon_sides, lat_sides
-
- def boundary(self, *, vertices_per_side=None, force_clockwise=False, frequency=None):
- """Retrieve the AreaBoundary object.
+ def _get_projection_sides(self, vertices_per_side: Optional[int] = None) -> tuple:
+ """Return the projection boundary sides of the current area.
+
+ Args:
+ vertices_per_side:
+ The number of points to provide for each side.
+ 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.
+
+ Returns:
+ The output structure is a tuple of two lists of four elements each.
+ The first list contains the projection x coordinates.
+ The second list contains the projection y coordinates.
+ 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 _is_any_corner_out_of_earth_disk(self):
+ # Geostationary
+ if self.is_geostationary:
+ return self._get_geostationary_boundary_sides(vertices_per_side=vertices_per_side,
+ coordinates="projection")
+ # Polar Projections, Global Planar Projections (Mollweide, Robinson)
+ # - Retrieve dummy right and left sides
+ 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):
+ """Retrieve the boundary object in projection coordinates.
+
+ If the CRS of the AreaDefinition is geographic, the returned boundary
+ object is a SphericalBoundary, otherwise a PlanarBoundary is returned.
Parameters
----------
vertices_per_side:
- The number of points to provide for each side. By default (None)
- the full width and height will be provided, except for geostationary
- projection where by default only 50 points are selected.
- frequency:
- Deprecated, use vertices_per_side
- force_clockwise:
- Perform 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
- on knowing the inside versus the outside of a polygon. These
- operations assume that coordinates are clockwise.
- Default is False.
+ The number of points to provide for each side.
+ 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..
"""
- from pyresample.boundary import AreaBoundary
- if frequency is not None:
- warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
- PendingDeprecationWarning, stacklevel=2)
- vertices_per_side = vertices_per_side or frequency
- if self.is_geostationary:
- lon_sides, lat_sides = self._get_geo_boundary_sides(vertices_per_side=vertices_per_side)
- else:
- lon_sides, lat_sides = self.get_bbox_lonlats(vertices_per_side=vertices_per_side,
- force_clockwise=force_clockwise)
- boundary = AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
- return boundary
+ from pyresample.boundary import PlanarBoundary
+ if self.crs.is_geographic:
+ return self.boundary(vertices_per_side=vertices_per_side)
+ return PlanarBoundary(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):
@@ -1628,9 +1779,12 @@ def get_edge_bbox_in_projection_coordinates(self, vertices_per_side: Optional[in
if frequency is not None:
warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
PendingDeprecationWarning, stacklevel=2)
+ warnings.warn("The `get_edge_bbox_in_projection_coordinates` method is pending deprecation."
+ "Use `area.projection_boundary().contour()` instead.",
+ PendingDeprecationWarning, stacklevel=2)
vertices_per_side = vertices_per_side or frequency
- x, y = self._get_bbox_elements(self.get_proj_coords, vertices_per_side)
- return np.hstack(x), np.hstack(y)
+ x, y = self.projection_boundary(vertices_per_side=vertices_per_side).contour(closed=True)
+ return x, y
@property
def area_extent(self):
@@ -2704,27 +2858,47 @@ def get_geostationary_angle_extent(geos_area):
def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
"""Get the bbox in geos projection coordinates of the valid pixels inside `geos_area`.
+ NOTE: The area_bbox must be defined with many vertices to avoid loosing
+ vertices at the corners later on in _get_geostationary_boundary_sides.
+
Args:
geos_area: Geostationary area definition to get the bounding box for.
nb_points: Number of points on the polygon.
"""
+ from shapely.geometry import Polygon
+
x, y = get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points)
- ll_x, ll_y, ur_x, ur_y = geos_area.area_extent
- from shapely.geometry import Polygon
geo_bbox = Polygon(np.vstack((x, y)).T)
- area_bbox = Polygon(((ll_x, ll_y), (ll_x, ur_y), (ur_x, ur_y), (ur_x, ll_y)))
+ area_bbox = _create_area_extent_polygon(geos_area.area_extent, nb_points=nb_points)
intersection = area_bbox.intersection(geo_bbox)
try:
x, y = intersection.boundary.xy
except NotImplementedError:
- return [], []
+ # 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])
+def _create_area_extent_polygon(area_extent, nb_points=50):
+ """Create the area_extent polygon with nb_points vertices per side."""
+ from shapely.geometry import Polygon
+
+ # Generate points for each edge
+ ll_x, ll_y, ur_x, ur_y = area_extent
+ bottom = [(x, ll_y) for x in np.linspace(ll_x, ur_x, nb_points + 2)]
+ right = [(ur_x, y) for y in np.linspace(ll_y, ur_y, nb_points + 2)][1:]
+ top = [(x, ur_y) for x in np.linspace(ur_x, ll_x, nb_points + 2)][1:]
+ left = [(ll_x, y) for y in np.linspace(ur_y, ll_y, nb_points + 2)][1:-1]
+
+ # Combine points to form the area_extent polygon
+ return Polygon(bottom + right + top + left)
+
+
def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
- """Get the bbox in geos projection coordinates of the full disk in `geos_area` projection.
+ """Get the valid boundary geos projection coordinates of the full disk.
Args:
geos_area: Geostationary area definition to get the bounding box for.
@@ -2740,7 +2914,24 @@ def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
y = -np.sin(points_around) * (y_max_angle - 0.0001)
x *= h
y *= h
+ return x, y
+
+
+def _get_geostationary_bounding_box(geos_area, coordinates="geographic", nb_points=50):
+ """Get the bounding box coordinates of the valid pixels inside `geos_area`.
+ If coordinates='geographic', it returns the lat/lon coordinates.
+ If coordinates='projection', it returns the projection coordinates.
+
+ Args:
+ geos_area: Geostationary area definition to get the bounding box for.
+ coordinates: Whether to retrieve geographic or projection coordinates.
+ nb_points: Number of points on the polygon
+ """
+ x, y = get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points)
+ if coordinates == "geographic":
+ lons, lats = Proj(geos_area.crs)(x, y, inverse=True)
+ return lons, lats
return x, y
@@ -2751,9 +2942,12 @@ def get_geostationary_bounding_box_in_lonlats(geos_area, nb_points=50):
geos_area: Geostationary area definition to get the bounding box for.
nb_points: Number of points on the polygon
"""
- x, y = get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points)
- lons, lats = Proj(geos_area.crs)(x, y, inverse=True)
- return lons, lats
+ warnings.warn("'get_geostationary_bounding_box_in_lonlats' is deprecated. Please call "
+ "'area.boundary().contour()' instead.",
+ DeprecationWarning, stacklevel=2)
+ return _get_geostationary_bounding_box(geos_area,
+ coordinates="geographic",
+ nb_points=nb_points)
def get_geostationary_bounding_box(geos_area, nb_points=50):
@@ -2764,10 +2958,21 @@ def get_geostationary_bounding_box(geos_area, nb_points=50):
nb_points: Number of points on the polygon
"""
- warnings.warn("'get_geostationary_bounding_box' is deprecated. Please use "
- "'get_geostationary_bounding_box_in_lonlats' instead.",
+ warnings.warn("'get_geostationary_bounding_box' is deprecated. Please call "
+ "'area.boundary().contour()' instead.",
DeprecationWarning, stacklevel=2)
- return get_geostationary_bounding_box_in_lonlats(geos_area, nb_points)
+ return get_geostationary_bounding_box(geos_area,
+ coordinates="geographic",
+ nb_points=nb_points)
+
+
+def _is_any_corner_out_of_earth_disk(area_def):
+ """Check if the area corners are out of the Earth disk."""
+ try:
+ _ = area_def.corners
+ return False
+ except Exception:
+ return True
def combine_area_extents_vertical(area1, area2):
diff --git a/pyresample/gradient/__init__.py b/pyresample/gradient/__init__.py
index ed7fb8270..63d3b2ece 100644
--- a/pyresample/gradient/__init__.py
+++ b/pyresample/gradient/__init__.py
@@ -36,11 +36,7 @@
from shapely.geometry import Polygon
from pyresample import CHUNK_SIZE
-from pyresample.geometry import (
- AreaDefinition,
- SwathDefinition,
- get_geostationary_bounding_box_in_lonlats,
-)
+from pyresample.geometry import AreaDefinition, SwathDefinition
from pyresample.gradient._gradient_search import (
one_step_gradient_indices,
one_step_gradient_search,
@@ -107,7 +103,7 @@ def _get_projection_coordinates(self, datachunks):
chunks=datachunks)
src_crs = self.source_geo_def.crs
self.use_input_coords = True
- except AttributeError:
+ except NotImplementedError:
self.src_x, self.src_y = self.source_geo_def.get_lonlats(
chunks=datachunks)
src_crs = pyproj.CRS.from_string("+proj=longlat")
@@ -116,7 +112,7 @@ def _get_projection_coordinates(self, datachunks):
self.dst_x, self.dst_y = self.target_geo_def.get_proj_coords(
chunks=CHUNK_SIZE)
dst_crs = self.target_geo_def.crs
- except AttributeError as err:
+ except NotImplementedError as err:
if self.use_input_coords is False:
raise NotImplementedError('Cannot resample lon/lat to lon/lat with gradient search.') from err
self.dst_x, self.dst_y = self.target_geo_def.get_lonlats(
@@ -133,17 +129,22 @@ def _get_projection_coordinates(self, datachunks):
src_prj=src_crs, dst_prj=dst_crs)
self.prj = pyproj.Proj(self.target_geo_def.crs)
+ def _get_prj_poly(self, geo_def):
+ # - None if out of Earth Disk
+ # - False is SwathDefinition
+ if isinstance(geo_def, SwathDefinition):
+ return False
+ try:
+ poly = _get_polygon(self.prj, geo_def)
+ except Exception:
+ poly = None
+ return poly
+
def _get_src_poly(self, src_y_start, src_y_end, src_x_start, src_x_end):
"""Get bounding polygon for source chunk."""
geo_def = self.source_geo_def[src_y_start:src_y_end,
src_x_start:src_x_end]
- try:
- src_poly = get_polygon(self.prj, geo_def)
- except AttributeError:
- # Can't create polygons for SwathDefinition
- src_poly = False
-
- return src_poly
+ return self._get_prj_poly(geo_def)
def _get_dst_poly(self, idx, dst_x_start, dst_x_end,
dst_y_start, dst_y_end):
@@ -152,13 +153,8 @@ def _get_dst_poly(self, idx, dst_x_start, dst_x_end,
if dst_poly is None:
geo_def = self.target_geo_def[dst_y_start:dst_y_end,
dst_x_start:dst_x_end]
- try:
- dst_poly = get_polygon(self.prj, geo_def)
- except AttributeError:
- # Can't create polygons for SwathDefinition
- dst_poly = False
+ dst_poly = self._get_prj_poly(geo_def)
self.dst_polys[idx] = dst_poly
-
return dst_poly
def get_chunk_mappings(self):
@@ -166,6 +162,8 @@ def get_chunk_mappings(self):
src_y_chunks, src_x_chunks = self.src_x.chunks
dst_y_chunks, dst_x_chunks = self.dst_x.chunks
+ dst_is_swath = isinstance(self.target_geo_def, SwathDefinition)
+
coverage_status = []
src_slices, dst_slices = [], []
dst_mosaic_locations = []
@@ -179,7 +177,6 @@ def get_chunk_mappings(self):
# Get source chunk polygon
src_poly = self._get_src_poly(src_y_start, src_y_end,
src_x_start, src_x_end)
-
dst_x_start = 0
for x_step_number, dst_x_step in enumerate(dst_x_chunks):
dst_x_end = dst_x_start + dst_x_step
@@ -187,9 +184,18 @@ def get_chunk_mappings(self):
for y_step_number, dst_y_step in enumerate(dst_y_chunks):
dst_y_end = dst_y_start + dst_y_step
# Get destination chunk polygon
- dst_poly = self._get_dst_poly((x_step_number, y_step_number),
- dst_x_start, dst_x_end,
- dst_y_start, dst_y_end)
+ # - Retrieve if source chunk poly is inside Earth Disk
+ # - src_poly = None
+ # - Skips lot polygon creations when source is GEO FD
+ # - Retrieve if dst area is swath
+ # - Currently dst_poly will be False
+ # - check_overlap will return True
+ if src_poly is not None or dst_is_swath:
+ dst_poly = self._get_dst_poly((x_step_number, y_step_number),
+ dst_x_start, dst_x_end,
+ dst_y_start, dst_y_end)
+ else:
+ dst_poly = None
covers = check_overlap(src_poly, dst_poly)
@@ -300,13 +306,15 @@ def compute(self, data, fill_value=None, **kwargs):
def check_overlap(src_poly, dst_poly):
"""Check if the two polygons overlap."""
+ # swath definition case
if dst_poly is False or src_poly is False:
covers = True
+ # area / area case
elif dst_poly is not None and src_poly is not None:
covers = src_poly.intersects(dst_poly)
+ # out of earth disk case
else:
covers = False
-
return covers
@@ -373,21 +381,17 @@ def _check_input_coordinates(dst_x, dst_y,
raise ValueError("Target arrays should all have the same shape")
-def get_border_lonlats(geo_def: AreaDefinition):
+def _get_border_lonlats(geo_def: AreaDefinition, vertices_per_side=None):
"""Get the border x- and y-coordinates."""
if geo_def.is_geostationary:
- lon_b, lat_b = get_geostationary_bounding_box_in_lonlats(geo_def, 3600)
- else:
- lons, lats = geo_def.get_boundary_lonlats()
- lon_b = np.concatenate((lons.side1, lons.side2, lons.side3, lons.side4))
- lat_b = np.concatenate((lats.side1, lats.side2, lats.side3, lats.side4))
-
+ vertices_per_side = 3600
+ lon_b, lat_b = geo_def.boundary(vertices_per_side=vertices_per_side).contour(closed=True)
return lon_b, lat_b
-def get_polygon(prj, geo_def):
+def _get_polygon(prj, geo_def, vertices_per_side=None):
"""Get border polygon from area definition in projection *prj*."""
- lon_b, lat_b = get_border_lonlats(geo_def)
+ lon_b, lat_b = _get_border_lonlats(geo_def, vertices_per_side=vertices_per_side)
x_borders, y_borders = prj(lon_b, lat_b)
boundary = [(x_borders[i], y_borders[i]) for i in range(len(x_borders))
if np.isfinite(x_borders[i]) and np.isfinite(y_borders[i])]
diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py
index fa43001bd..68f060e21 100644
--- a/pyresample/kd_tree.py
+++ b/pyresample/kd_tree.py
@@ -413,15 +413,18 @@ def _get_valid_input_index(source_geo_def,
source_is_coord = isinstance(source_geo_def, geometry.CoordinateDefinition)
if (source_is_coord or source_is_griddish) and target_is_griddish:
# Resampling from swath to grid or from grid to grid
- lonlat_boundary = target_geo_def.get_boundary_lonlats()
-
- # Combine reduced and legal values
- valid_input_index &= \
- data_reduce.get_valid_index_from_lonlat_boundaries(
- lonlat_boundary[0],
- lonlat_boundary[1],
- source_lons, source_lats,
- radius_of_influence)
+ # - If invalid sides, return np.ones
+ try:
+ sides_lons, sides_lats = target_geo_def.boundary().sides
+ # Combine reduced and legal values
+ valid_input_index &= \
+ data_reduce.get_valid_index_from_lonlat_boundaries(
+ sides_lons,
+ sides_lats,
+ source_lons, source_lats,
+ radius_of_influence)
+ except Exception:
+ valid_input_index = np.ones(source_lons.size, dtype=bool)
if isinstance(valid_input_index, np.ma.core.MaskedArray):
# Make sure valid_input_index is not a masked array
@@ -440,15 +443,19 @@ def _get_valid_output_index(source_geo_def, target_geo_def, target_lons,
geometry.AreaDefinition)) and \
isinstance(target_geo_def, geometry.CoordinateDefinition):
# Resampling from grid to swath
- lonlat_boundary = source_geo_def.get_boundary_lonlats()
- valid_output_index = \
- data_reduce.get_valid_index_from_lonlat_boundaries(
- lonlat_boundary[0],
- lonlat_boundary[1],
- target_lons,
- target_lats,
- radius_of_influence)
- valid_output_index = valid_output_index.astype(bool)
+ # - If invalid sides, return np.ones
+ try:
+ sides_lons, sides_lats = source_geo_def.boundary().sides
+ valid_output_index = \
+ data_reduce.get_valid_index_from_lonlat_boundaries(
+ sides_lons,
+ sides_lats,
+ target_lons,
+ target_lats,
+ radius_of_influence)
+ valid_output_index = valid_output_index.astype(bool)
+ except Exception:
+ valid_output_index = np.ones(target_lons.size, dtype=bool)
# Remove illegal values
valid_out = ((target_lons >= -180) & (target_lons <= 180) & (target_lats <= 90) & (target_lats >= -90))
diff --git a/pyresample/slicer.py b/pyresample/slicer.py
index 5306818d8..d4937e58d 100644
--- a/pyresample/slicer.py
+++ b/pyresample/slicer.py
@@ -27,11 +27,7 @@
from pyproj.enums import TransformDirection
from pyresample import AreaDefinition, SwathDefinition
-from pyresample.geometry import (
- IncompatibleAreas,
- InvalidArea,
- get_geostationary_bounding_box_in_proj_coords,
-)
+from pyresample.geometry import IncompatibleAreas, InvalidArea
try:
import dask.array as da
@@ -99,7 +95,7 @@ class SwathSlicer(Slicer):
def get_polygon_to_contain(self):
"""Get the shapely Polygon corresponding to *area_to_contain* in lon/lat coordinates."""
from shapely.geometry import Polygon
- x, y = self.area_to_contain.get_edge_bbox_in_projection_coordinates(10)
+ x, y = self.area_to_contain.projection_boundary(vertices_per_side=10).contour(closed=True)
poly = Polygon(zip(*self._transformer.transform(x, y)))
return poly
@@ -127,16 +123,12 @@ def _assemble_slices(chunk_slices):
def _get_chunk_polygons_for_swath_to_crop(swath_to_crop):
"""Get the polygons for each chunk of the area_to_crop."""
res = []
- from shapely.geometry import Polygon
src_chunks = swath_to_crop.lons.chunks
for _position, (line_slice, col_slice) in _enumerate_chunk_slices(src_chunks):
line_slice = expand_slice(line_slice)
col_slice = expand_slice(col_slice)
smaller_swath = swath_to_crop[line_slice, col_slice]
- lons, lats = smaller_swath.get_edge_lonlats(10)
- lons = np.hstack(lons)
- lats = np.hstack(lats)
- smaller_poly = Polygon(zip(lons, lats))
+ smaller_poly = smaller_swath.boundary(vertices_per_side=10).to_shapely_polygon()
res.append((smaller_poly, (line_slice, col_slice)))
return res
@@ -152,9 +144,13 @@ class AreaSlicer(Slicer):
def get_polygon_to_contain(self):
"""Get the shapely Polygon corresponding to *area_to_contain* in projection coordinates of *area_to_crop*."""
from shapely.geometry import Polygon
- x, y = self.area_to_contain.get_edge_bbox_in_projection_coordinates(frequency=10)
+ x, y = self.area_to_contain.projection_boundary(vertices_per_side=10).contour(closed=True)
if self.area_to_crop.is_geostationary:
- x_geos, y_geos = get_geostationary_bounding_box_in_proj_coords(self.area_to_crop, 360)
+ geo_boundary = self.area_to_crop.projection_boundary(vertices_per_side=360)
+ x_geos, y_geos = geo_boundary.contour(closed=True)
+ # 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))
@@ -163,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:
@@ -175,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
- from shapely.geometry import Polygon
- poly_to_crop = Polygon(zip(*self.area_to_crop.get_edge_bbox_in_projection_coordinates(frequency=10)))
- 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]),
@@ -201,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_area_config.py b/pyresample/test/test_area_config.py
index 70f615f81..8f17bdb0b 100644
--- a/pyresample/test/test_area_config.py
+++ b/pyresample/test/test_area_config.py
@@ -24,6 +24,10 @@
import numpy as np
+from pyresample import _root_path
+
+TEST_FILES_PATH = os.path.join(_root_path, "test", 'test_files')
+
class TestLegacyAreaParser(unittest.TestCase):
"""Test legacy .cfg parsing."""
@@ -31,7 +35,7 @@ class TestLegacyAreaParser(unittest.TestCase):
def test_area_parser_legacy(self):
"""Test legacy area parser."""
from pyresample import parse_area_file
- ease_nh, ease_sh = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'),
+ ease_nh, ease_sh = parse_area_file(os.path.join(TEST_FILES_PATH, 'areas.cfg'),
'ease_nh', 'ease_sh')
# pyproj 2.0+ adds some extra parameters
@@ -63,7 +67,7 @@ def test_area_parser_legacy(self):
def test_load_area(self):
from pyresample import load_area
- ease_nh = load_area(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh')
+ ease_nh = load_area(os.path.join(TEST_FILES_PATH, 'areas.cfg'), 'ease_nh')
# pyproj 2.0+ adds some extra parameters
projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', "
"'no_defs': 'None', 'proj': 'laea', 'type': 'crs', "
@@ -87,11 +91,11 @@ def test_area_file_not_found_exception(self):
def test_not_found_exception(self):
from pyresample.area_config import AreaNotFound, parse_area_file
self.assertRaises(AreaNotFound, parse_area_file,
- os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'no_area')
+ os.path.join(TEST_FILES_PATH, 'areas.cfg'), 'no_area')
def test_commented(self):
from pyresample import parse_area_file
- areas = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'))
+ areas = parse_area_file(os.path.join(TEST_FILES_PATH, 'areas.cfg'))
self.assertNotIn('commented', [area.name for area in areas])
@@ -101,7 +105,7 @@ class TestYAMLAreaParser(unittest.TestCase):
def test_area_parser_yaml(self):
"""Test YAML area parser."""
from pyresample import parse_area_file
- test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area_file = os.path.join(TEST_FILES_PATH, 'areas.yaml')
test_areas = parse_area_file(test_area_file, 'ease_nh', 'ease_sh', 'test_meters', 'test_degrees',
'test_latlong')
ease_nh, ease_sh, test_m, test_deg, test_latlong = test_areas
@@ -158,7 +162,7 @@ def test_dynamic_area_parser_yaml(self):
"""Test YAML area parser on dynamic areas."""
from pyresample import parse_area_file
from pyresample.geometry import DynamicAreaDefinition
- test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area_file = os.path.join(TEST_FILES_PATH, 'areas.yaml')
test_area = parse_area_file(test_area_file, 'test_dynamic_resolution')[0]
self.assertIsInstance(test_area, DynamicAreaDefinition)
@@ -168,7 +172,7 @@ def test_dynamic_area_parser_yaml(self):
# lat/lon
from pyresample import parse_area_file
from pyresample.geometry import DynamicAreaDefinition
- test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area_file = os.path.join(TEST_FILES_PATH, 'areas.yaml')
test_area = parse_area_file(test_area_file, 'test_dynamic_resolution_ll')[0]
self.assertIsInstance(test_area, DynamicAreaDefinition)
@@ -178,7 +182,7 @@ def test_dynamic_area_parser_yaml(self):
def test_dynamic_area_parser_passes_resolution(self):
"""Test that the resolution from the file is passed to a dynamic area."""
from pyresample import parse_area_file
- test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area_file = os.path.join(TEST_FILES_PATH, 'areas.yaml')
test_area = parse_area_file(test_area_file, 'omerc_bb_1000')[0]
assert test_area.resolution == (1000, 1000)
diff --git a/pyresample/test/test_boundary.py b/pyresample/test/test_boundary/test_area_boundary.py
similarity index 72%
rename from pyresample/test/test_boundary.py
rename to pyresample/test/test_boundary/test_area_boundary.py
index e90fa9b03..77fe66f74 100644
--- a/pyresample/test/test_boundary.py
+++ b/pyresample/test/test_boundary/test_area_boundary.py
@@ -16,7 +16,7 @@
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see .
-"""Test the boundary objects."""
+"""Test the AreaBoundary objects."""
import unittest
import numpy as np
@@ -30,22 +30,22 @@ class TestAreaBoundary(unittest.TestCase):
def test_creation_from_lonlat_sides(self):
"""Test AreaBoundary creation from sides."""
- lon_sides = [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])]
- lat_sides = [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])]
+ 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 AreaBoundary
- boundary = AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
+ boundary = AreaBoundary.from_lonlat_sides(sides_lons, sides_lats)
# Assert sides coincides
- for b_lon, src_lon in zip(boundary.sides_lons, lon_sides):
+ for b_lon, src_lon in zip(boundary.sides_lons, sides_lons):
assert np.allclose(b_lon, src_lon)
- for b_lat, src_lat in zip(boundary.sides_lats, lat_sides):
+ for b_lat, src_lat in zip(boundary.sides_lats, sides_lats):
assert np.allclose(b_lat, src_lat)
def test_creation(self):
@@ -54,17 +54,17 @@ def test_creation(self):
(np.array([2., 3.]), np.array([7., 8.])),
(np.array([3., 3.5, 4.]), np.array([8., 8.5, 9.])),
(np.array([4., 1.]), np.array([9., 6.]))]
- lon_sides = [side[0]for side in list_sides]
- lat_sides = [side[1]for side in list_sides]
+ sides_lons = [side[0]for side in list_sides]
+ sides_lats = [side[1]for side in list_sides]
# Define AreaBoundary
boundary = AreaBoundary(*list_sides)
# Assert sides coincides
- for b_lon, src_lon in zip(boundary.sides_lons, lon_sides):
+ for b_lon, src_lon in zip(boundary.sides_lons, sides_lons):
assert np.allclose(b_lon, src_lon)
- for b_lat, src_lat in zip(boundary.sides_lats, lat_sides):
+ for b_lat, src_lat in zip(boundary.sides_lats, sides_lats):
assert np.allclose(b_lat, src_lat)
def test_number_sides_required(self):
@@ -78,16 +78,16 @@ def test_number_sides_required(self):
def test_vertices_property(self):
"""Test AreaBoundary vertices property."""
- lon_sides = [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])]
- lat_sides = [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])]
+ 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 AreaBoundary
- boundary = AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
+ boundary = AreaBoundary.from_lonlat_sides(sides_lons, sides_lats)
# Assert vertices
expected_vertices = np.array([[1., 6.],
diff --git a/pyresample/test/test_boundary/test_geographic_boundary.py b/pyresample/test/test_boundary/test_geographic_boundary.py
new file mode 100644
index 000000000..3990faada
--- /dev/null
+++ b/pyresample/test/test_boundary/test_geographic_boundary.py
@@ -0,0 +1,98 @@
+#!/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 SphericalBoundary objects."""
+
+import numpy as np
+import pytest
+
+from pyresample import SwathDefinition
+from pyresample.boundary import SphericalBoundary
+
+
+class TestSwathSphericalBoundary():
+ """Test 'SphericalBoundary' class for SwathDefinition."""
+
+ 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]),
+ 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])]
+ 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 SphericalBoundary creation."""
+ # Define SphericalBoundary
+ boundary = SphericalBoundary(self.area)
+
+ # Assert sides coincides
+ 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, self.sides_lats):
+ assert np.allclose(b_lat, src_lat)
+
+ def test_minimum_swath_size(self):
+ """Test swath has minimum 2 dimensions."""
+ area = SwathDefinition(self.lons[0], self.lats[1])
+ with pytest.raises(ValueError):
+ SphericalBoundary(area)
+
+ def test_vertices_property(self):
+ """Test SphericalBoundary vertices property."""
+ # Define SphericalBoundary
+ boundary = SphericalBoundary(self.area)
+ # Assert vertices
+ expected_vertices = np.array([[1., 6.],
+ [1.5, 6.5],
+ [2., 7.],
+ [3., 8.],
+ [3.5, 8.5],
+ [4., 9.]])
+ assert np.allclose(boundary.vertices, expected_vertices)
+
+ def test_contour(self):
+ """Test that SphericalBoundary.contour(closed=False) returns the correct (lon,lat) tuple."""
+ # Define SphericalBoundary
+ boundary = SphericalBoundary(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 SphericalBoundary.contour(closed=True) returns the correct (lon,lat) tuple."""
+ # Define SphericalBoundary
+ boundary = SphericalBoundary(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..2a8fdd863
--- /dev/null
+++ b/pyresample/test/test_boundary/test_order.py
@@ -0,0 +1,144 @@
+#!/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 ordering checks."""
+
+import pytest
+
+from pyresample.boundary.spherical_boundary import _is_clockwise_order
+
+
+class Test_Is_Clockwise_Order:
+ """Test clockwise check."""
+
+ 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)
diff --git a/pyresample/test/test_boundary/test_sides.py b/pyresample/test/test_boundary/test_sides.py
new file mode 100644
index 000000000..29de28265
--- /dev/null
+++ b/pyresample/test/test_boundary/test_sides.py
@@ -0,0 +1,107 @@
+#!/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 BoundarySides objects."""
+
+import numpy as np
+import pytest
+
+from pyresample.boundary.sides import BoundarySides
+
+
+class TestBoundarySides:
+ """Test suite for the BoundarySides class with 1D numpy arrays for sides."""
+
+ def test_initialization_valid_input(self):
+ """Test initialization with valid 1D numpy array inputs."""
+ sides = [np.array([1, 2, 3]), # top
+ np.array([3, 4, 5]), # right
+ np.array([5, 6, 7]), # bottom
+ np.array([7, 8, 1])] # left
+ boundary = BoundarySides(sides)
+ assert all(np.array_equal(boundary[i], sides[i]) for i in range(4))
+
+ def test_initialization_invalid_input(self):
+ """Test initialization with invalid inputs, such as wrong number of sides or non-1D arrays."""
+ with pytest.raises(ValueError):
+ BoundarySides([np.array([1, 2]), # Invalid number of sides
+ np.array([2, 3])])
+
+ with pytest.raises(ValueError):
+ BoundarySides([np.array([1, 2]), # Non-1D arrays
+ np.array([[2, 3], [4, 5]]),
+ np.array([5, 6]),
+ np.array([6, 7])])
+
+ with pytest.raises(ValueError):
+ BoundarySides([np.array([1, 2]), # Invalid side connection
+ np.array([3, 4]),
+ np.array([4, 6]),
+ np.array([6, 1])])
+
+ def test_property_accessors(self):
+ """Test property accessors with 1D numpy arrays."""
+ sides = [np.array([1, 2, 3]), # top
+ np.array([3, 4, 5]), # right
+ np.array([5, 6, 7]), # bottom
+ np.array([7, 8, 1])] # left
+ boundary = BoundarySides(sides)
+ assert np.array_equal(boundary.top, sides[0])
+ assert np.array_equal(boundary.right, sides[1])
+ assert np.array_equal(boundary.bottom, sides[2])
+ assert np.array_equal(boundary.left, sides[3])
+
+ def test_vertices_property(self):
+ """Test the vertices property with concatenated 1D numpy arrays."""
+ sides = [np.array([1, 2, 3]), # top
+ np.array([3, 4, 5]), # right
+ np.array([5, 6, 7]), # bottom
+ np.array([7, 8, 1])] # left
+ boundary = BoundarySides(sides)
+ expected_vertices = np.array([1, 2, 3, 4, 5, 6, 7, 8])
+ assert np.array_equal(boundary.vertices, expected_vertices)
+
+ def test_iteration(self):
+ """Test iteration over the 1D numpy array sides."""
+ sides = [np.array([1, 2, 3]), # top
+ np.array([3, 4, 5]), # right
+ np.array([5, 6, 7]), # bottom
+ np.array([7, 8, 1])] # left
+ boundary = BoundarySides(sides)
+ for i, side in enumerate(boundary):
+ assert np.array_equal(side, sides[i])
+
+ def test_indexing_valid(self):
+ """Test valid indexing with 1D numpy arrays."""
+ sides = [np.array([1, 2, 3]), # top
+ np.array([3, 4, 5]), # right
+ np.array([5, 6, 7]), # bottom
+ np.array([7, 8, 1])] # left
+ boundary = BoundarySides(sides)
+ for i in range(4):
+ assert np.array_equal(boundary[i], sides[i])
+
+ def test_indexing_invalid(self):
+ """Test indexing with invalid indices."""
+ sides = [np.array([1, 2, 3]), # top
+ np.array([3, 4, 5]), # right
+ np.array([5, 6, 7]), # bottom
+ np.array([7, 8, 1])] # left
+ boundary = BoundarySides(sides)
+ with pytest.raises(IndexError):
+ boundary[4] # Invalid index
diff --git a/pyresample/test/test_boundary/test_utils.py b/pyresample/test/test_boundary/test_utils.py
new file mode 100644
index 000000000..7e4558ac6
--- /dev/null
+++ b/pyresample/test/test_boundary/test_utils.py
@@ -0,0 +1,146 @@
+#!/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_contour_indices,
+ find_boundary_mask,
+ 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}")
diff --git a/pyresample/test/test_data_reduce.py b/pyresample/test/test_data_reduce.py
index 478a7ad51..5f3adc10c 100644
--- a/pyresample/test/test_data_reduce.py
+++ b/pyresample/test/test_data_reduce.py
@@ -73,9 +73,9 @@ def test_reduce_boundary(self):
lambda y, x: -180 + (360.0 / 1000) * x, (1000, 1000))
lats = np.fromfunction(
lambda y, x: -90 + (180.0 / 1000) * y, (1000, 1000))
- boundary_lonlats = self.area_def.get_boundary_lonlats()
- lons, lats, data = swath_from_lonlat_boundaries(boundary_lonlats[0],
- boundary_lonlats[1],
+ sides_lons, sides_lats = self.area_def.boundary().sides
+ lons, lats, data = swath_from_lonlat_boundaries(sides_lons,
+ sides_lats,
lons, lats, data, 7000)
cross_sum = data.sum()
expected = 20685125.0
diff --git a/pyresample/test/test_geometry/test_area.py b/pyresample/test/test_geometry/test_area.py
index d1aeda7ee..fc005748a 100644
--- a/pyresample/test/test_geometry/test_area.py
+++ b/pyresample/test/test_geometry/test_area.py
@@ -25,13 +25,13 @@
from pyproj import CRS, Proj
import pyresample
+import pyresample.geometry
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,
get_full_geostationary_bounding_box_in_proj_coords,
get_geostationary_angle_extent,
- get_geostationary_bounding_box_in_lonlats,
- get_geostationary_bounding_box_in_proj_coords,
ignore_pyproj_proj_warnings,
)
from pyresample.future.geometry.base import get_array_hashable
@@ -39,49 +39,265 @@
from pyresample.test.utils import assert_future_geometry
+# BUG in 'area_class' fixture
+# --> Overwrite create_test_area here
+def create_test_area(crs, shape, area_extent, **kwargs):
+ """Create an AreaDefinition object for testing."""
+ area = AreaDefinition(crs=crs, shape=shape, area_extent=area_extent, **kwargs)
+ return area
+
+
@pytest.fixture
-def stere_area(create_test_area):
+def stere_area():
"""Create basic polar-stereographic area definition."""
+ proj_dict = {
+ 'a': '6378144.0',
+ 'b': '6356759.0',
+ 'lat_0': '50.00',
+ 'lat_ts': '50.00',
+ 'lon_0': '8.00',
+ 'proj': 'stere'
+ }
+ shape = (800, 800)
+ area_extent = (-1370912.72, -909968.64000000001, 1029087.28, 1490031.3600000001)
return create_test_area(
- {
- 'a': '6378144.0',
- 'b': '6356759.0',
- 'lat_0': '50.00',
- 'lat_ts': '50.00',
- 'lon_0': '8.00',
- 'proj': 'stere'
- },
- 800,
- 800,
- (-1370912.72, -909968.64000000001, 1029087.28, 1490031.3600000001),
+ proj_dict,
+ shape,
+ area_extent,
attrs={"name": 'areaD'},
)
@pytest.fixture
-def geos_src_area(create_test_area):
+def global_lonlat_antimeridian_centered_area():
+ """Create global lonlat projection area centered on the -180 antimeridian."""
+ shape = (4, 4)
+ area_extent = (0, -90.0, 360, 90.0)
+ proj_dict = '+proj=longlat +pm=180'
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def global_platee_caree_area():
+ """Create global platee projection area."""
+ shape = (4, 4)
+ area_extent = (-180.0, -90.0, 180.0, 90.0)
+ proj_dict = 'EPSG:4326'
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def global_platee_caree_minimum_area():
+ """Create minimum size global platee projection area."""
+ """Create global platee projection area."""
+ shape = (2, 2)
+ area_extent = (-180.0, -90.0, 180.0, 90.0)
+ proj_dict = 'EPSG:4326'
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def local_platee_caree_area():
+ """Create local platee projection area."""
+ shape = (4, 4)
+ area_extent = (100, 20, 120, 40)
+ proj_dict = 'EPSG:4326'
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def local_lonlat_antimeridian_centered_area():
+ """Create local lonlat projection area centered on the -180 antimeridian."""
+ shape = (4, 4)
+ area_extent = (100, 20, 120, 40)
+ proj_dict = '+proj=longlat +pm=180'
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def local_meter_area():
+ """Create local meter projection area."""
+ shape = (2, 2)
+ area_extent = (2_600_000.0, 1_050_000, 2_800_000.0, 1_170_000)
+ proj_dict = 'EPSG:2056'
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def south_pole_area():
+ """Create projection area centered on south pole."""
+ shape = (2, 2)
+ area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+ proj_dict = {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def north_pole_area():
+ """Create projection area centered on north pole."""
+ shape = (2, 2)
+ area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+ proj_dict = {'proj': 'laea', 'lat_0': 90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def geos_src_area():
"""Create basic geostationary area definition."""
- x_size = 3712
- y_size = 3712
+ shape = (3712, 3712)
area_extent = (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)
proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'h': 35785831.0,
'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}
return create_test_area(
- proj_dict,
- x_size,
- y_size,
- area_extent,
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
)
@pytest.fixture
-def laea_area(create_test_area):
+def geos_fd_area():
+ """Create full disc geostationary area definition."""
+ shape = (100, 100)
+ area_extent = (-5500000., -5500000., 5500000., 5500000.)
+ proj_dict = {'a': 6378169.00, 'b': 6356583.80, 'h': 35785831.0,
+ 'lon_0': 0, 'proj': 'geos', 'units': 'm'}
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def geos_out_disk_area():
+ """Create out of Earth diskc geostationary area definition."""
+ shape = (10, 10)
+ area_extent = (-5500000., -5500000., -5300000., -5300000.)
+ proj_dict = {'a': 6378169.00, 'b': 6356583.80, 'h': 35785831.0,
+ 'lon_0': 0, 'proj': 'geos', 'units': 'm'}
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def geos_half_out_disk_area():
+ """Create geostationary area definition with portion of boundary out of earth_disk."""
+ shape = (100, 100)
+ area_extent = (-5500000., -10000., 0, 10000.)
+ proj_dict = {'a': 6378169.00, 'b': 6356583.80, 'h': 35785831.0,
+ 'lon_0': 0, 'proj': 'geos', 'units': 'm'}
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def geos_conus_area():
+ """Create CONUS geostationary area definition (portion is out-of-Earth disk)."""
+ shape = (30, 50) # (3000, 5000) for GOES-R CONUS/PACUS
+ proj_dict = {'h': 35786023, 'sweep': 'x', 'x_0': 0, 'y_0': 0,
+ 'ellps': 'GRS80', 'no_defs': None, 'type': 'crs',
+ 'lon_0': -75, 'proj': 'geos', 'units': 'm'}
+ area_extent = (-3627271.29128, 1583173.65752, 1382771.92872, 4589199.58952)
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def geos_mesoscale_area():
+ """Create CONUS geostationary area definition."""
+ shape = (10, 10) # (1000, 1000) for GOES-R mesoscale
+ proj_dict = {'h': 35786023, 'sweep': 'x', 'x_0': 0, 'y_0': 0,
+ 'ellps': 'GRS80', 'no_defs': None, 'type': 'crs',
+ 'lon_0': -75, 'proj': 'geos', 'units': 'm'}
+ area_extent = (-501004.322, 3286588.35232, 501004.322, 4288596.99632)
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def truncated_geos_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)
+ shape = (1392, 3712)
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def truncated_geos_area_in_space():
+ """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)
+ shape = (10, 10)
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
+
+
+@pytest.fixture
+def laea_area():
"""Create basic LAEA area definition."""
- x_size = 10
- y_size = 10
+ shape = (10, 10)
area_extent = [1000000, 0, 1050000, 50000]
proj_dict = {"proj": 'laea', 'lat_0': '60', 'lon_0': '0', 'a': '6371228.0', 'units': 'm'}
- return create_test_area(proj_dict, x_size, y_size, area_extent)
+ return create_test_area(
+ crs=proj_dict,
+ shape=shape,
+ area_extent=area_extent,
+ )
class TestAreaHashability:
@@ -1423,30 +1639,6 @@ def test_create_area_def_nonpole_center(self):
assert area_def.shape == (101, 90)
-@pytest.fixture
-def truncated_geos_area(create_test_area):
- """Create a truncated geostationary area."""
- projection = {'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)
- width = 3712
- height = 1392
- geos_area = create_test_area(projection, width, height, area_extent)
- return geos_area
-
-
-@pytest.fixture
-def truncated_geos_area_in_space(create_test_area):
- """Create a truncated geostationary area."""
- projection = {'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)
- width = 10
- height = 10
- geos_area = create_test_area(projection, width, height, area_extent)
- return geos_area
-
-
class TestGeostationaryTools:
"""Test the geostationary bbox tools."""
@@ -1476,44 +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
@@ -1530,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,
@@ -1555,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):
@@ -1596,6 +1785,23 @@ def test_get_geostationary_angle_extent(self):
np.testing.assert_allclose(expected,
get_geostationary_angle_extent(geos_area))
+ @pytest.mark.parametrize('area_def_name,corner_out_of_disk', [
+ ("geos_fd_area", True),
+ ("geos_out_disk_area", True),
+ ("geos_half_out_disk_area", True),
+ ("geos_conus_area", True),
+ ("geos_mesoscale_area", False),
+ ])
+ def test_is_any_corner_out_of_earth_disk(self, request, area_def_name, corner_out_of_disk):
+ """Test if corner area is out of Earth disk."""
+ from pyresample.geometry import _is_any_corner_out_of_earth_disk
+
+ area_def = request.getfixturevalue(area_def_name)
+ if corner_out_of_disk:
+ assert _is_any_corner_out_of_earth_disk(area_def)
+ else:
+ assert not _is_any_corner_out_of_earth_disk(area_def)
+
class TestCrop:
"""Test the area helpers."""
@@ -1881,14 +2087,30 @@ def test_unsupported_slice_inputs(self, create_test_area, create_test_swath, swa
class TestBoundary:
"""Test 'boundary' method for AreaDefinition classes."""
- def test_polar_south_pole_projection(self, create_test_area):
+ @pytest.mark.parametrize('area_def_name,assert_is_called', [
+ ("geos_fd_area", True),
+ ("geos_out_disk_area", True),
+ ("geos_half_out_disk_area", True),
+ ("geos_conus_area", True),
+ ("geos_mesoscale_area", False),
+ ])
+ def test_get_geographic_sides_call_geostationary_utility(self, request, area_def_name, assert_is_called):
+ area_def = request.getfixturevalue(area_def_name)
+
+ with patch.object(area_def, '_get_geostationary_boundary_sides') as mock_get_geo:
+
+ # Call the method that could trigger the geostationary _get_geostationary_boundary_sides
+ _ = area_def._get_geographic_sides(vertices_per_side=None)
+ # Assert _get_geostationary_boundary_sides was not called
+ if assert_is_called:
+ mock_get_geo.assert_called_once()
+ else:
+ mock_get_geo.assert_not_called()
+
+ def test_polar_south_pole_projection(self, south_pole_area):
"""Test boundary for polar projection around the South Pole."""
- areadef = create_test_area(
- {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'},
- 2, 2,
- (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625),
- )
- boundary = areadef.boundary(force_clockwise=False)
+ areadef = south_pole_area
+ boundary = areadef.boundary()
# Check boundary shape
height, width = areadef.shape
@@ -1902,14 +2124,11 @@ def test_polar_south_pole_projection(self, create_test_area):
[-135., -55.61313895]])
np.testing.assert_allclose(expected_vertices, boundary.vertices)
- def test_nort_pole_projection(self, create_test_area):
+ def test_north_pole_projection(self, north_pole_area):
"""Test boundary for polar projection around the North Pole."""
- areadef = create_test_area(
- {'proj': 'laea', 'lat_0': 90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'},
- 2, 2,
- (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625),
- )
- boundary = areadef.boundary(force_clockwise=False)
+ areadef = north_pole_area
+
+ boundary = areadef.boundary()
# Check boundary shape
height, width = areadef.shape
@@ -1923,34 +2142,30 @@ def test_nort_pole_projection(self, create_test_area):
[-45., 55.61313895]])
np.testing.assert_allclose(expected_vertices, boundary.vertices)
- def test_geostationary_projection(self, create_test_area):
+ def test_full_disc_geostationary_projection(self, geos_fd_area):
"""Test boundary for geostationary projection."""
- areadef = create_test_area(
- {'a': 6378169.00, 'b': 6356583.80, 'h': 35785831.00, 'lon_0': 0, 'proj': 'geos'},
- 100, 100,
- (-5500000., -5500000., 5500000., 5500000.),
- )
+ areadef = geos_fd_area
# Check default boundary shape
default_n_vertices = 50
- boundary = areadef.boundary(frequency=None)
+ boundary = areadef.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.boundary(frequency=n_vertices)
+ boundary = areadef.boundary(vertices_per_side=n_vertices, )
assert boundary.vertices.shape == (minimum_n_vertices, 2)
- # Check odd frequency number
+ # Check odd number of vertices per side
# - Rounded to the sequent even number (to construct the sides)
n_odd_vertices = 5
- boundary = areadef.boundary(frequency=n_odd_vertices)
+ boundary = areadef.boundary(vertices_per_side=n_odd_vertices)
assert boundary.vertices.shape == (n_odd_vertices + 1, 2)
# Check boundary vertices
n_vertices = 10
- boundary = areadef.boundary(frequency=n_vertices, force_clockwise=False)
+ boundary = areadef.boundary(vertices_per_side=n_vertices, )
# Check boundary vertices is in correct order
expected_vertices = np.array([[-7.54251621e+01, 3.53432890e+01],
@@ -1965,14 +2180,10 @@ def test_geostationary_projection(self, create_test_area):
[-7.92337283e+01, 6.94302533e-15]])
np.testing.assert_allclose(expected_vertices, boundary.vertices)
- def test_global_platee_caree_projection(self, create_test_area):
+ def test_global_platee_caree_projection(self, global_platee_caree_area):
"""Test boundary for global platee caree projection."""
- areadef = create_test_area(
- 'EPSG:4326',
- 4, 4,
- (-180.0, -90.0, 180.0, 90.0),
- )
- boundary = areadef.boundary(force_clockwise=False)
+ areadef = global_platee_caree_area
+ boundary = areadef.boundary()
# Check boundary shape
height, width = areadef.shape
@@ -1994,14 +2205,10 @@ def test_global_platee_caree_projection(self, create_test_area):
[-135., 22.5]])
np.testing.assert_allclose(expected_vertices, boundary.vertices)
- def test_minimal_global_platee_caree_projection(self, create_test_area):
+ def test_minimal_global_platee_caree_projection(self, global_platee_caree_minimum_area):
"""Test boundary for global platee caree projection."""
- areadef = create_test_area(
- 'EPSG:4326',
- 2, 2,
- (-180.0, -90.0, 180.0, 90.0),
- )
- boundary = areadef.boundary(force_clockwise=False)
+ areadef = global_platee_caree_minimum_area
+ boundary = areadef.boundary()
# Check boundary shape
height, width = areadef.shape
@@ -2015,14 +2222,10 @@ def test_minimal_global_platee_caree_projection(self, create_test_area):
[-90., -45.]])
np.testing.assert_allclose(expected_vertices, boundary.vertices)
- def test_local_area_projection(self, create_test_area):
+ def test_local_area_projection(self, local_meter_area):
"""Test local area projection in meter."""
- areadef = create_test_area(
- 'EPSG:2056',
- 2, 2,
- (2_600_000.0, 1_050_000, 2_800_000.0, 1_170_000),
- )
- boundary = areadef.boundary(force_clockwise=False)
+ areadef = local_meter_area
+ boundary = areadef.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 47a4a81ed..90d795820 100644
--- a/pyresample/test/test_geometry/test_swath.py
+++ b/pyresample/test/test_geometry/test_swath.py
@@ -152,7 +152,7 @@ def test_swath_def_bbox(
def test_swath_def_bbox_decimated(self, create_test_swath):
swath_def = _gen_swath_def_numpy(create_test_swath)
- bbox_lons, bbox_lats = swath_def.get_bbox_lonlats(frequency=None)
+ bbox_lons, bbox_lats = swath_def.get_bbox_lonlats(vertices_per_side=None)
assert len(bbox_lons) == len(bbox_lats)
assert len(bbox_lons) == 4
assert len(bbox_lons[0]) == 10
@@ -160,7 +160,7 @@ def test_swath_def_bbox_decimated(self, create_test_swath):
assert len(bbox_lons[2]) == 10
assert len(bbox_lons[3]) == 50
- bbox_lons, bbox_lats = swath_def.get_bbox_lonlats(frequency=5)
+ bbox_lons, bbox_lats = swath_def.get_bbox_lonlats(vertices_per_side=5)
assert len(bbox_lons) == len(bbox_lats)
assert len(bbox_lons) == 4
assert len(bbox_lons[0]) == 5
@@ -403,14 +403,20 @@ def test_get_edge_lonlats(self, create_test_swath):
[81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
area = create_test_swath(lons, lats)
lons, lats = area.get_edge_lonlats()
- np.testing.assert_allclose(lons, [-90.67900085, 79.11000061, 81.26400757,
- 81.26400757, 29.67200089, 10.26000023,
- 10.26000023, -5.10700035, -21.52500153,
- -21.52500153, -21.56500053, -90.67900085])
- np.testing.assert_allclose(lats, [85.23900604, 80.84000397, 67.07600403,
- 67.07600403, 54.14700317, 30.54700089,
- 30.54700089, 34.0850029, 35.58000183,
- 35.58000183, 62.25600433, 85.23900604])
+ expected_lons = [
+ -90.67900085, 79.11000061, # 81.26400757,
+ 81.26400757, 29.67200089, # 10.26000023,
+ 10.26000023, -5.10700035, # -21.52500153,
+ -21.52500153, -21.56500053, # -90.67900085,
+ ]
+ expected_lats = [
+ 85.23900604, 80.84000397, # 67.07600403,
+ 67.07600403, 54.14700317, # 30.54700089,
+ 30.54700089, 34.0850029, # 35.58000183,
+ 35.58000183, 62.25600433, # 85.23900604,
+ ]
+ np.testing.assert_allclose(lons, expected_lons)
+ np.testing.assert_allclose(lats, expected_lats)
lats = np.array([[80., 80., 80.],
[80., 90., 80],
@@ -420,10 +426,20 @@ def test_get_edge_lonlats(self, create_test_swath):
[-135., -180., 135.]]).T
area = create_test_swath(lons, lats)
lons, lats = area.get_edge_lonlats()
- np.testing.assert_allclose(lons, [-45., -90., -135., -135., -180., 135.,
- 135., 90., 45., 45., 0., -45.])
- np.testing.assert_allclose(lats, [80., 80., 80., 80., 80., 80., 80.,
- 80., 80., 80., 80., 80.])
+ expected_lons = [
+ -45., -90., # -135.,
+ -135., -180., # 135.,
+ 135., 90., # 45.,
+ 45., 0., # -45.
+ ]
+ expected_lats = [
+ 80., 80., # 80.,
+ 80., 80., # 80.,
+ 80., 80., # 80.,
+ 80., 80., # 80.
+ ]
+ np.testing.assert_allclose(lons, expected_lons)
+ np.testing.assert_allclose(lats, expected_lats)
def test_compute_optimal_bb(self, create_test_swath):
"""Test computing the bb area."""
@@ -586,9 +602,9 @@ def test_swath_definition(self, create_test_swath):
lats = np.array([[65.9, 65.86, 65.82, 65.78],
[65.89, 65.86, 65.82, 65.78]])
- # Define SwathDefinition and retrieve AreaBoundary
+ # Define SwathDefinition and retrieve SphericalBoundary
swath_def = create_test_swath(lons, lats)
- boundary = swath_def.boundary(force_clockwise=False)
+ boundary = swath_def.boundary()
# Check boundary shape
height, width = swath_def.shape
diff --git a/pyresample/test/test_gradient.py b/pyresample/test/test_gradient.py
index a283e63d1..addc3be1c 100644
--- a/pyresample/test/test_gradient.py
+++ b/pyresample/test/test_gradient.py
@@ -123,7 +123,7 @@ def test_get_src_poly_area(self):
self.resampler._get_projection_coordinates(chunks)
self.resampler._get_gradients()
poly = self.resampler._get_src_poly(0, 40, 0, 40)
- assert np.allclose(poly.area, 12365358458842.43)
+ assert np.allclose(poly.area, 12365352661227.719)
def test_get_src_poly_swath(self):
"""Test defining source chunk polygon for SwathDefinition."""
@@ -134,24 +134,19 @@ def test_get_src_poly_swath(self):
poly = self.swath_resampler._get_src_poly(0, 40, 0, 40)
assert poly is False
- @mock.patch('pyresample.gradient.get_polygon')
- def test_get_dst_poly(self, get_polygon):
+ @mock.patch('pyresample.gradient._get_polygon')
+ def test_get_dst_poly(self, _get_polygon):
"""Test defining destination chunk polygon."""
chunks = (10, 10)
self.resampler._get_projection_coordinates(chunks)
self.resampler._get_gradients()
- # First call should make a call to get_polygon()
+ # First call should make a call to _get_polygon()
self.resampler._get_dst_poly('idx1', 0, 10, 0, 10)
- assert get_polygon.call_count == 1
+ assert _get_polygon.call_count == 1
assert 'idx1' in self.resampler.dst_polys
# The second call to the same index should come from cache
self.resampler._get_dst_poly('idx1', 0, 10, 0, 10)
- assert get_polygon.call_count == 1
-
- # Swath defs raise AttributeError, and False is returned
- get_polygon.side_effect = AttributeError
- self.resampler._get_dst_poly('idx2', 0, 10, 0, 10)
- assert self.resampler.dst_polys['idx2'] is False
+ assert _get_polygon.call_count == 1
def test_filter_data(self):
"""Test filtering chunks that do not overlap."""
@@ -602,43 +597,42 @@ def test_check_overlap():
assert check_overlap(poly1, poly2) is False
-def test_get_border_lonlats_geos():
- """Test that correct methods are called in get_border_lonlats() with geos inputs."""
- from pyresample.gradient import get_border_lonlats
+def test__get_border_lonlats_geos():
+ """Test that correct methods are called in _get_border_lonlats() with geos inputs."""
+ from pyresample.gradient import _get_border_lonlats
+ sides_lons = [np.array([1, 2]), np.array([2, 3]), np.array([3, 4]), np.array([4, 1])]
+ sides_lats = [np.array([1, 2]), np.array([2, 3]), np.array([3, 4]), np.array([4, 1])]
geo_def = AreaDefinition("", "", "",
"+proj=geos +h=1234567", 2, 2, [1, 2, 3, 4])
- with mock.patch("pyresample.gradient.get_geostationary_bounding_box_in_lonlats") as get_geostationary_bounding_box:
- get_geostationary_bounding_box.return_value = 1, 2
- res = get_border_lonlats(geo_def)
- assert res == (1, 2)
- get_geostationary_bounding_box.assert_called_with(geo_def, 3600)
-
-
-def test_get_border_lonlats():
- """Test that correct methods are called in get_border_lonlats()."""
- from pyresample.boundary import SimpleBoundary
- from pyresample.gradient import get_border_lonlats
- lon_sides = SimpleBoundary(side1=np.array([1]), side2=np.array([2]),
- side3=np.array([3]), side4=np.array([4]))
- lat_sides = SimpleBoundary(side1=np.array([1]), side2=np.array([2]),
- side3=np.array([3]), side4=np.array([4]))
+ with mock.patch.object(geo_def, "_get_geographic_sides") as mock_get_geographic_sides:
+ mock_get_geographic_sides.return_value = sides_lons, sides_lats
+ lon_b, lat_b = _get_border_lonlats(geo_def)
+ np.testing.assert_allclose(lon_b, np.array([1, 2, 3, 4, 1]))
+ np.testing.assert_allclose(lat_b, np.array([1, 2, 3, 4, 1]))
+
+
+def test__get_border_lonlats():
+ """Test that correct methods are called in _get_border_lonlats()."""
+ from pyresample.gradient import _get_border_lonlats
+ sides_lons = [np.array([1, 2]), np.array([2, 3]), np.array([3, 4]), np.array([4, 1])]
+ sides_lats = [np.array([1, 2]), np.array([2, 3]), np.array([3, 4]), np.array([4, 1])]
geo_def = AreaDefinition("", "", "",
"+proj=lcc +lat_1=25 +lat_2=25", 2, 2, [1, 2, 3, 4])
- with mock.patch.object(geo_def, "get_boundary_lonlats") as get_boundary_lonlats:
- get_boundary_lonlats.return_value = lon_sides, lat_sides
- lon_b, lat_b = get_border_lonlats(geo_def)
- assert np.all(lon_b == np.array([1, 2, 3, 4]))
- assert np.all(lat_b == np.array([1, 2, 3, 4]))
+ with mock.patch.object(geo_def, "_get_geographic_sides") as mock_get_geographic_sides:
+ mock_get_geographic_sides.return_value = sides_lons, sides_lats
+ lon_b, lat_b = _get_border_lonlats(geo_def)
+ np.testing.assert_allclose(lon_b, np.array([1, 2, 3, 4, 1]))
+ np.testing.assert_allclose(lat_b, np.array([1, 2, 3, 4, 1]))
@mock.patch('pyresample.gradient.Polygon')
-@mock.patch('pyresample.gradient.get_border_lonlats')
-def test_get_polygon(get_border_lonlats, Polygon):
+@mock.patch('pyresample.gradient._get_border_lonlats')
+def test__get_polygon(_get_border_lonlats, Polygon):
"""Test polygon creation."""
- from pyresample.gradient import get_polygon
+ from pyresample.gradient import _get_polygon
# Valid polygon
- get_border_lonlats.return_value = (1, 2)
+ _get_border_lonlats.return_value = (1, 2)
geo_def = mock.MagicMock()
prj = mock.MagicMock()
x_borders = [0, 0, 1, 1]
@@ -647,8 +641,7 @@ def test_get_polygon(get_border_lonlats, Polygon):
prj.return_value = (x_borders, y_borders)
poly = mock.MagicMock(area=2.0)
Polygon.return_value = poly
- res = get_polygon(prj, geo_def)
- get_border_lonlats.assert_called_with(geo_def)
+ res = _get_polygon(prj, geo_def)
prj.assert_called_with(1, 2)
Polygon.assert_called_with(boundary)
assert res is poly
@@ -658,18 +651,18 @@ def test_get_polygon(get_border_lonlats, Polygon):
y_borders = [-1, 0, np.nan, 1, 1, np.nan, -1]
boundary = [(0, 0), (0, 1), (1, 1), (2, -1)]
prj.return_value = (x_borders, y_borders)
- res = get_polygon(prj, geo_def)
+ res = _get_polygon(prj, geo_def)
Polygon.assert_called_with(boundary)
assert res is poly
# Polygon area is NaN
poly.area = np.nan
- res = get_polygon(prj, geo_def)
+ res = _get_polygon(prj, geo_def)
assert res is None
# Polygon area is 0.0
poly.area = 0.0
- res = get_polygon(prj, geo_def)
+ res = _get_polygon(prj, geo_def)
assert res is None
diff --git a/pyresample/test/test_image.py b/pyresample/test/test_image.py
index dbfb06da4..5c0273e39 100644
--- a/pyresample/test/test_image.py
+++ b/pyresample/test/test_image.py
@@ -22,7 +22,9 @@
import numpy
-from pyresample import geometry, image, utils
+from pyresample import _root_path, geometry, image, utils
+
+TEST_FILES_PATH = os.path.join(_root_path, "test", 'test_files')
class Test(unittest.TestCase):
@@ -108,8 +110,7 @@ def test_masked_image(self):
area_con = msg_con.resample(self.area_def)
res = area_con.image_data
resampled_mask = res.mask.astype('int')
- expected = numpy.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files', 'mask_grid.dat'),
+ expected = numpy.fromfile(os.path.join(TEST_FILES_PATH, 'mask_grid.dat'),
sep=' ').reshape((800, 800))
self.assertTrue(numpy.array_equal(resampled_mask, expected))
@@ -123,9 +124,7 @@ def test_masked_image_fill(self):
area_con = msg_con.resample(self.area_def)
res = area_con.image_data
resampled_mask = res.mask.astype('int')
- expected = numpy.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
- 'mask_grid.dat'),
+ expected = numpy.fromfile(os.path.join(TEST_FILES_PATH, 'mask_grid.dat'),
sep=' ').reshape((800, 800))
self.assertTrue(numpy.array_equal(resampled_mask, expected))
@@ -146,7 +145,7 @@ def test_nearest_resize(self):
area_con = msg_con.resample(self.msg_area_resize)
res = area_con.image_data
cross_sum = res.sum()
- expected = 2212023.0175830
+ expected = 2211981.7030199994
self.assertAlmostEqual(cross_sum, expected)
def test_nearest_neighbour_multi(self):
diff --git a/pyresample/test/test_kd_tree.py b/pyresample/test/test_kd_tree.py
index e26663241..45c84f1b6 100644
--- a/pyresample/test/test_kd_tree.py
+++ b/pyresample/test/test_kd_tree.py
@@ -23,9 +23,11 @@
import numpy as np
import pytest
-from pyresample import geometry, kd_tree, utils
+from pyresample import _root_path, geometry, kd_tree, utils
from pyresample.test.utils import catch_warnings
+TEST_FILES_PATH = os.path.join(_root_path, "test", 'test_files')
+
class Test(unittest.TestCase):
"""Test nearest neighbor resampling on numpy arrays."""
@@ -496,12 +498,10 @@ def test_masked_nearest(self):
masked_data = np.ma.array(data, mask=mask)
res = kd_tree.resample_nearest(swath_def, masked_data.ravel(),
self.area_def, 50000, segments=1)
- expected_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_nearest_mask.dat'),
sep=' ').reshape((800, 800))
- expected_data = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_data = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_nearest_data.dat'),
sep=' ').reshape((800, 800))
self.assertTrue(np.array_equal(expected_mask, res.mask))
@@ -531,12 +531,10 @@ def test_masked_gauss(self):
masked_data = np.ma.array(data, mask=mask)
res = kd_tree.resample_gauss(swath_def, masked_data.ravel(),
self.area_def, 50000, 25000, segments=1)
- expected_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_mask.dat'),
sep=' ').reshape((800, 800))
- expected_data = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_data = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_data.dat'),
sep=' ').reshape((800, 800))
expected = expected_data.sum()
@@ -552,8 +550,7 @@ def test_masked_fill_float(self):
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
res = kd_tree.resample_nearest(swath_def, data.ravel(),
self.area_def, 50000, fill_value=None, segments=1)
- expected_fill_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_fill_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_fill_value.dat'),
sep=' ').reshape((800, 800))
fill_mask = res.mask
@@ -566,8 +563,7 @@ def test_masked_fill_int(self):
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
res = kd_tree.resample_nearest(swath_def, data.ravel(),
self.area_def, 50000, fill_value=None, segments=1)
- expected_fill_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_fill_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_fill_value.dat'),
sep=' ').reshape((800, 800))
fill_mask = res.mask
@@ -586,8 +582,7 @@ def test_masked_full(self):
masked_data.ravel(
), self.area_def, 50000,
fill_value=None, segments=1)
- expected_fill_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_fill_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_full_fill.dat'),
sep=' ').reshape((800, 800))
fill_mask = res.mask
@@ -614,8 +609,7 @@ def test_masked_full_multi(self):
res = kd_tree.resample_nearest(swath_def,
masked_data, self.area_def, 50000,
fill_value=None, segments=1)
- expected_fill_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_fill_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_full_fill_multi.dat'),
sep=' ').reshape((800, 800, 3))
fill_mask = res.mask
@@ -750,8 +744,7 @@ def test_masked_multi_from_sample(self):
valid_input_index,
valid_output_index, index_array,
fill_value=None)
- expected_fill_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
- 'test_files',
+ expected_fill_mask = np.fromfile(os.path.join(TEST_FILES_PATH,
'mask_test_full_fill_multi.dat'),
sep=' ').reshape((800, 800, 3))
fill_mask = res.mask
diff --git a/pyresample/test/test_plot.py b/pyresample/test/test_plot.py
index 3220f88c2..b7b1ef362 100644
--- a/pyresample/test/test_plot.py
+++ b/pyresample/test/test_plot.py
@@ -24,6 +24,8 @@
import numpy as np
+from pyresample import _root_path
+
try:
import matplotlib
matplotlib.use('Agg')
@@ -36,6 +38,8 @@
Basemap = None
+TEST_FILES_PATH = os.path.join(_root_path, "test", 'test_files')
+
MERIDIANS1 = np.array([-180, -170, -160, -150, -140, -130, -120, -110, -100, -90, -80,
-70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30,
40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140,
@@ -48,8 +52,7 @@
class Test(unittest.TestCase):
"""Test the plot utilities."""
- filename = os.path.abspath(os.path.join(os.path.dirname(__file__),
- 'test_files', 'ssmis_swath.npz'))
+ filename = os.path.join(TEST_FILES_PATH, 'ssmis_swath.npz')
data = np.load(filename)['data']
lons = data[:, 0].astype(np.float64)
lats = data[:, 1].astype(np.float64)
@@ -79,7 +82,7 @@ def test_ellps2axis(self):
def test_area_def2basemap(self):
"""Test the area to Basemap object conversion function."""
from pyresample import parse_area_file, plot
- area_def = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml'), 'ease_sh')[0]
+ area_def = parse_area_file(os.path.join(TEST_FILES_PATH, 'areas.yaml'), 'ease_sh')[0]
bmap = plot.area_def2basemap(area_def)
self.assertTrue(bmap.rmajor == bmap.rminor and bmap.rmajor == 6371228.0,
'Failed to create Basemap object')
@@ -158,7 +161,8 @@ def test_translate_coast_res(self):
def test_plate_carreeplot(self):
"""Test the Plate Caree plotting functionality."""
from pyresample import geometry, kd_tree, parse_area_file, plot
- area_def = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml'), 'pc_world')[0]
+
+ area_def = parse_area_file(os.path.join(TEST_FILES_PATH, 'areas.yaml'), 'pc_world')[0]
swath_def = geometry.SwathDefinition(self.lons, self.lats)
result = kd_tree.resample_nearest(swath_def, self.tb37v, area_def,
radius_of_influence=20000,
@@ -171,7 +175,7 @@ def test_plate_carreeplot(self):
def test_easeplot(self):
"""Test the plotting on the ease grid area."""
from pyresample import geometry, kd_tree, parse_area_file, plot
- area_def = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml'), 'ease_sh')[0]
+ area_def = parse_area_file(os.path.join(TEST_FILES_PATH, 'areas.yaml'), 'ease_sh')[0]
swath_def = geometry.SwathDefinition(self.lons, self.lats)
result = kd_tree.resample_nearest(swath_def, self.tb37v, area_def,
radius_of_influence=20000,
@@ -181,7 +185,7 @@ def test_easeplot(self):
def test_orthoplot(self):
"""Test the ortho plotting."""
from pyresample import geometry, kd_tree, parse_area_file, plot
- area_def = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ortho')[0]
+ area_def = parse_area_file(os.path.join(TEST_FILES_PATH, 'areas.cfg'), 'ortho')[0]
swath_def = geometry.SwathDefinition(self.lons, self.lats)
result = kd_tree.resample_nearest(swath_def, self.tb37v, area_def,
radius_of_influence=20000,
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."""
diff --git a/pyresample/test/test_swath.py b/pyresample/test/test_swath.py
index 23c1d6de6..7db5e5a7a 100644
--- a/pyresample/test/test_swath.py
+++ b/pyresample/test/test_swath.py
@@ -23,17 +23,18 @@
import numpy as np
-from pyresample import geometry, kd_tree
+from pyresample import _root_path, geometry, kd_tree
from pyresample.test.utils import catch_warnings
+TEST_FILES_PATH = os.path.join(_root_path, "test", 'test_files')
+
warnings.simplefilter("always")
class Test(unittest.TestCase):
"""Tests for swath definitions."""
- filename = os.path.abspath(os.path.join(os.path.dirname(__file__),
- 'test_files', 'ssmis_swath.npz'))
+ filename = os.path.join(TEST_FILES_PATH, 'ssmis_swath.npz')
data = np.load(filename)['data']
lons = data[:, 0].astype(np.float64)
lats = data[:, 1].astype(np.float64)
diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py
index 34b4948f8..fd6719164 100644
--- a/pyresample/test/test_utils.py
+++ b/pyresample/test/test_utils.py
@@ -29,6 +29,7 @@
from pyproj import CRS
import pyresample
+from pyresample import _root_path
from pyresample.test.utils import (
assert_future_geometry,
create_test_latitude,
@@ -37,6 +38,8 @@
from pyresample.utils import load_cf_area
from pyresample.utils.row_appendable_array import RowAppendableArray
+TEST_FILES_PATH = os.path.join(_root_path, "test", 'test_files')
+
def tmptiff(width=100, height=100, transform=None, crs=None, dtype=np.uint8):
"""Create a temporary in-memory TIFF file of all ones."""
@@ -238,7 +241,7 @@ def test_def2yaml_converter(self):
import tempfile
from pyresample import convert_def_to_yaml, parse_area_file
- def_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg')
+ def_file = os.path.join(TEST_FILES_PATH, 'areas.cfg')
filehandle, yaml_file = tempfile.mkstemp()
os.close(filehandle)
try:
@@ -460,12 +463,12 @@ class TestLoadCFAreaPublic:
"""Test public API load_cf_area() for loading an AreaDefinition from netCDF/CF files."""
def test_load_cf_no_exist(self):
- cf_file = os.path.join(os.path.dirname(__file__), 'test_files', 'does_not_exist.nc')
+ cf_file = os.path.join(TEST_FILES_PATH, 'does_not_exist.nc')
with pytest.raises(FileNotFoundError):
load_cf_area(cf_file)
def test_load_cf_from_not_nc(self):
- cf_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ cf_file = os.path.join(TEST_FILES_PATH, 'areas.yaml')
with pytest.raises((ValueError, OSError)):
load_cf_area(cf_file)
diff --git a/pyresample/test/test_visualization/test_geometries.py b/pyresample/test/test_visualization/test_geometries.py
new file mode 100644
index 000000000..13f0c777a
--- /dev/null
+++ b/pyresample/test/test_visualization/test_geometries.py
@@ -0,0 +1,83 @@
+#!/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 cartopy plotting utilities."""
+
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+import pytest
+from shapely.geometry import Polygon
+
+from pyresample.visualization.geometries import (
+ _add_map_background,
+ _check_subplot_kw,
+ _initialize_plot,
+ plot_geometries,
+)
+
+
+class TestPlotFunctions:
+ """Test suite for the provided plotting functions."""
+
+ def test_add_map_background(self):
+ """Test adding a map background to an axis."""
+ fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
+ result_ax = _add_map_background(ax)
+ assert isinstance(result_ax, plt.Axes)
+
+ def test_check_subplot_kw_valid(self):
+ """Test _check_subplot_kw with valid input."""
+ valid_kw = {'projection': ccrs.PlateCarree()}
+ assert _check_subplot_kw(valid_kw) == valid_kw
+
+ def test_check_subplot_kw_none(self):
+ """Test _check_subplot_kw with None input."""
+ assert 'projection' in _check_subplot_kw(None)
+
+ def test_check_subplot_kw_invalid(self):
+ """Test _check_subplot_kw with invalid input."""
+ with pytest.raises(TypeError):
+ _check_subplot_kw("invalid")
+
+ with pytest.raises(TypeError):
+ _check_subplot_kw(2)
+
+ with pytest.raises(TypeError):
+ _check_subplot_kw([2])
+
+ with pytest.raises(ValueError):
+ _check_subplot_kw({})
+
+ def test_initialize_plot_with_ax(self):
+ """Test _initialize_plot with an existing ax."""
+ fig, ax = plt.subplots()
+ _, result_ax, initialized_here = _initialize_plot(ax=ax)
+ assert result_ax == ax
+ assert not initialized_here
+
+ @pytest.mark.parametrize("ax_provided", [True, False])
+ def test_plot_geometries(self, ax_provided):
+ """Test plot_geometries function returns the correct type based on ax_provided."""
+ import cartopy
+ vertices1 = [(0, 0), (0, 1), (1, 0)]
+ vertices2 = [(0, 0), (0, 2), (2, 0)]
+ geometries = [Polygon(vertices1), Polygon(vertices2)]
+ crs = ccrs.PlateCarree()
+ ax = plt.axes(projection=crs) if ax_provided else None
+ result = plot_geometries(geometries, crs, ax=ax)
+ assert isinstance(result, cartopy.mpl.feature_artist.FeatureArtist)
diff --git a/pyresample/visualization/__init__.py b/pyresample/visualization/__init__.py
new file mode 100644
index 000000000..40fd322f2
--- /dev/null
+++ b/pyresample/visualization/__init__.py
@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2014-2021 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 .
+"""pyresample tools for visualization."""
diff --git a/pyresample/visualization/geometries.py b/pyresample/visualization/geometries.py
new file mode 100644
index 000000000..cd027ceb9
--- /dev/null
+++ b/pyresample/visualization/geometries.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2010-2020 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 .
+"""Define how to plot a shapely geometry."""
+
+
+def _add_map_background(ax):
+ """Add cartopy map background."""
+ ax.stock_img()
+ ax.coastlines()
+ gl = ax.gridlines(draw_labels=True, linestyle="--")
+ gl.top_labels = False
+ gl.right_labels = False
+ return ax
+
+
+def _check_subplot_kw(subplot_kw):
+ """Check subplot_kw arguments."""
+ import cartopy.crs as ccrs
+
+ if subplot_kw is None:
+ subplot_kw = dict(projection=ccrs.PlateCarree())
+ if not isinstance(subplot_kw, dict):
+ raise TypeError("'subplot_kw' must be a dictionary.'")
+ if "projection" not in subplot_kw:
+ raise ValueError("Specify a cartopy 'projection' in subplot_kw.")
+ return subplot_kw
+
+
+def _initialize_plot(ax=None, subplot_kw=None):
+ """Initialize plot."""
+ import matplotlib.pyplot as plt
+
+ if ax is None:
+ subplot_kw = _check_subplot_kw(subplot_kw)
+ fig, ax = plt.subplots(subplot_kw=subplot_kw)
+ return fig, ax, True
+ else:
+ return None, ax, False
+
+
+def plot_geometries(geometries, crs, ax=None, subplot_kw=None, **kwargs):
+ """Plot geometries in cartopy."""
+ # Create figure if ax not provided
+ fig, ax, initialized_here = _initialize_plot(ax=ax, subplot_kw=subplot_kw)
+ # Add map background if ax not provided as input
+ if initialized_here:
+ ax = _add_map_background(ax)
+ # Add geometries
+ p = ax.add_geometries(geometries, crs=crs, **kwargs)
+ return p