Add 'antimeridian_mode' to DynamicAreaDefinition

pyresample/geometry.py

@@ -16,6 +16,7 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """Classes for geometry operations."""
+from __future__ import annotations
 
 import hashlib
 import math
@@ -24,11 +25,12 @@
 from functools import partial, wraps
 from logging import getLogger
 from pathlib import Path
-from typing import Optional
+from typing import Optional, Sequence, Union
 
 import numpy as np
 import yaml
 from pyproj import Geod, transform
+from pyproj.aoi import AreaOfUse
 
 from pyresample import CHUNK_SIZE
 from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj, Proj_MP
@@ -982,26 +984,59 @@ def pixel_size_y(self):
             return None
         return self.resolution[1]
 
-    def compute_domain(self, corners, resolution=None, shape=None):
+    def compute_domain(
+            self,
+            corners: Sequence,
+            resolution: Optional[Union[float, tuple[float, float]]] = None,
+            shape: Optional[tuple[int, int]] = None,
+            projection: Optional[Union[CRS, dict, str, int]] = None
+    ):
         """Compute shape and area_extent from corners and [shape or resolution] info.
 
-        Corners represents the center of pixels, while area_extent represents the edge of pixels.
+        Args:
+            corners:
+                4-element sequence representing the outer corners of the
+                region. Note that corners represents the center of pixels,
+                while area_extent represents the edge of pixels. The four
+                values are (xmin_corner, ymin_corner, xmax_corner, ymax_corner).
+                If the x corners are ``None`` then the full extent (area of use)
+                of the projection will be used. When needed, area of use is taken
+                from the PROJ library or in the case of a geographic lon/lat
+                projection -180/180 is used. A RuntimeError is raised if the
+                area of use is needed (when x corners are ``None``) and area
+                of use can't be determined.
+            resolution:
+                Spatial resolution in projection units (typically meters or
+                degrees). If not specified then shape must be provided.
+                If a scalar then it is treated as the x and y resolution. If
+                a tuple then x resolution is the first element, y is the
+                second.
+            shape:
+                Number of pixels in the area as a 2-element tuple. The first
+                is number of rows, the second number of columns.
+            projection:
+                PROJ.4 definition string, dictionary, integer EPSG code, or
+                pyproj CRS object.
 
         Note that ``shape`` is (rows, columns) and ``resolution`` is
         (x_size, y_size); the dimensions are flipped.
+
         """
         if resolution is not None and shape is not None:
             raise ValueError("Both resolution and shape can't be provided.")
         elif resolution is None and shape is None:
             raise ValueError("Either resolution or shape must be provided.")
+        if resolution is not None and isinstance(resolution, (int, float)):
+            resolution = (resolution, resolution)
+        if projection is None:
+            projection = self._projection
 
+        corners = self._update_corners_for_full_extent(corners, shape, resolution, projection)
         if shape:
             height, width = shape
             x_resolution = (corners[2] - corners[0]) * 1.0 / (width - 1)
             y_resolution = (corners[3] - corners[1]) * 1.0 / (height - 1)
         else:
-            if isinstance(resolution, (int, float)):
-                resolution = (resolution, resolution)
             x_resolution, y_resolution = resolution
             width = int(np.rint((corners[2] - corners[0]) * 1.0 / x_resolution + 1))
             height = int(np.rint((corners[3] - corners[1]) * 1.0 / y_resolution + 1))
@@ -1012,7 +1047,33 @@ def compute_domain(self, corners, resolution=None, shape=None):
                        corners[3] + y_resolution / 2)
         return area_extent, width, height
 
-    def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None):
+    def _update_corners_for_full_extent(self, corners, shape, resolution, projection):
+        corners = list(corners)
+        if corners[0] is not None:
+            return corners
+        aou = self._get_crs_area_of_use(projection)
+        if shape is not None:
+            width = shape[1]
+            x_resolution = (aou.east - aou.west) / width
+            corners[0] = aou.west + x_resolution / 2.0
+            corners[2] = aou.east - x_resolution / 2.0
+        else:
+            x_resolution = resolution[0]
+            corners[0] = aou.west + x_resolution / 2.0
+            corners[2] = aou.east - x_resolution / 2.0
+        return corners
+
+    def _get_crs_area_of_use(self, projection):
+        crs = CRS(projection)
+        aou = crs.area_of_use
+        if aou is None:
+            if crs.is_geographic:
+                return AreaOfUse(west=-180.0, south=-90.0, east=180.0, north=90.0)
+            raise RuntimeError("Projection has no defined area of use")
+        return aou
+
+    def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None,
+               antimeridian_mode=None):
         """Create an AreaDefinition from this area with help of some extra info.
 
         Parameters
@@ -1026,6 +1087,31 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None):
           the shape of the resulting area.
         proj_info:
           complementing parameters to the projection info.
+        antimeridian_mode:
+            How to handle lon/lat data crossing the anti-meridian of the
+            projection. This currently only affects lon/lat geographic
+            projections and data cases not covering the north or south pole.
+            The possible options are:
+
+            * "modify_extents": Set the X bounds to the edges of the data, but
+                add 360 to the right-most bound. This has the effect of making
+                the area coordinates continuous from the left side to the
+                right side. However, this means that some coordinates will be
+                outside the coordinate space of the projection. Although most
+                PROJ and pyresample functionality can handle this there may be
+                some edge cases.
+            * "modify_crs": Change the prime meridian of the projection
+                from 0 degrees longitude to 180 degrees longitude. This has
+                the effect of putting the data on a continuous coordinate
+                system. However, this means that comparing data resampled to
+                this resulting area and an area not over the anti-meridian
+                would be more difficult.
+            * "global_extents": Ignore the bounds of the data and use -180/180
+                degrees as the west and east bounds of the data. This will
+                generate a large output area, but with the benefit of keeping
+                the data on the original projection. Note that some resampling
+                methods may produce artifacts when resampling on the edge of
+                the area (the anti-meridian).
 
         Shape parameters are ignored if the instance is created
         with the `optimize_projection` flag set to True.
@@ -1047,76 +1133,56 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None):
         shape = None if None in shape else shape
         area_extent = self.area_extent
         if not area_extent or not width or not height:
-            corners = self._compute_bound_centers(proj_dict, lonslats)
-            area_extent, width, height = self.compute_domain(corners, resolution, shape)
+            projection, corners = self._compute_bound_centers(proj_dict, lonslats, antimeridian_mode=antimeridian_mode)
+            area_extent, width, height = self.compute_domain(corners, resolution, shape, projection)
         return AreaDefinition(self.area_id, self.description, '',
                               projection, width, height,
                               area_extent, self.rotation)
 
-    def _compute_bound_centers(self, proj_dict, lonslats):
+    def _compute_bound_centers(self, proj_dict, lonslats, antimeridian_mode):
+        from pyresample.utils.proj4 import DaskFriendlyTransformer
+
         lons, lats = self._extract_lons_lats(lonslats)
-        if hasattr(lons, 'compute'):
-            return self._compute_bound_centers_dask(proj_dict, lons, lats)
-        return self._compute_bound_centers_numpy(proj_dict, lons, lats)
-
-    def _compute_bound_centers_numpy(self, proj_dict, lons, lats):
-        # TODO: Do more dask-friendly things here
-        proj4 = Proj(proj_dict)
-        xarr, yarr = proj4(np.asarray(lons), np.asarray(lats))
+        crs = CRS(proj_dict)
+        transformer = DaskFriendlyTransformer.from_crs(CRS(4326), crs, always_xy=True)
+        xarr, yarr = transformer.transform(lons, lats)
         xarr[xarr > 9e29] = np.nan
         yarr[yarr > 9e29] = np.nan
         xmin = np.nanmin(xarr)
         xmax = np.nanmax(xarr)
         ymin = np.nanmin(yarr)
         ymax = np.nanmax(yarr)
-        x_passes_antimeridian = (xmax - xmin) > 355
-        epsilon = 0.1
-        y_is_pole = (ymax >= 90 - epsilon) or (ymin <= -90 + epsilon)
-        if proj4.crs.is_geographic and x_passes_antimeridian and not y_is_pole:
-            # cross anti-meridian of projection
-            xmin = np.nanmin(xarr[xarr >= 0])
-            xmax = np.nanmax(xarr[xarr < 0]) + 360
-        return xmin, ymin, xmax, ymax
-
-    def _compute_bound_centers_dask(self, proj_dict, lons, lats):
-        import dask.array as da
-
-        from pyresample.utils.proj4 import DaskFriendlyTransformer
-        crs = CRS(proj_dict)
-        transformer = DaskFriendlyTransformer.from_crs(CRS(4326), crs,
-                                                       always_xy=True)
-        xarr, yarr = transformer.transform(lons, lats)
-        xarr = da.where(xarr > 9e29, np.nan, xarr)
-        yarr = da.where(yarr > 9e29, np.nan, yarr)
-        _xmin = np.nanmin(xarr)
-        _xmax = np.nanmax(xarr)
-        _ymin = np.nanmin(yarr)
-        _ymax = np.nanmax(yarr)
-        xmin, xmax, ymin, ymax = da.compute(
-            _xmin,
-            _xmax,
-            _ymin,
-            _ymax)
-
+        if hasattr(lons, "compute"):
+            xmin, xmax, ymin, ymax = da.compute(xmin, xmax, ymin, ymax)
         x_passes_antimeridian = (xmax - xmin) > 355
         epsilon = 0.1
         y_is_pole = (ymax >= 90 - epsilon) or (ymin <= -90 + epsilon)
         if crs.is_geographic and x_passes_antimeridian and not y_is_pole:
             # cross anti-meridian of projection
-            xarr_pos = da.where(xarr >= 0, xarr, np.nan)
-            xarr_neg = da.where(xarr < 0, xarr, np.nan)
-            xmin = np.nanmin(xarr_pos)
-            xmax = np.nanmax(xarr_neg) + 360
-            xmin, xmax = da.compute(xmin, xmax)
-        return xmin, ymin, xmax, ymax
-
-    def _extract_lons_lats(self, lonslats):
+            xmin, xmax = self._compute_new_x_corners_for_antimeridian(xarr, antimeridian_mode)
+            if antimeridian_mode == "modify_crs":
+                proj_dict.update({"pm": 180.0})
+        return proj_dict, (xmin, ymin, xmax, ymax)
+
+    @staticmethod
+    def _extract_lons_lats(lonslats):
         try:
             lons, lats = lonslats
         except (TypeError, ValueError):
             lons, lats = lonslats.get_lonlats()
         return lons, lats
 
+    def _compute_new_x_corners_for_antimeridian(self, xarr, antimeridian_mode):
+        if antimeridian_mode == "global_extents":
+            xmin, xmax = (None, None)
+        else:
+            wrapped_array = xarr % 360
+            xmin = np.nanmin(wrapped_array)
+            xmax = np.nanmax(wrapped_array)
+            if hasattr(wrapped_array, "compute"):
+                xmin, xmax = da.compute(xmin, xmax)
+        return xmin, xmax
+
 
 def _invproj(data_x, data_y, proj_dict):
     """Perform inverse projection."""

pyresample/test/test_geometry.py

@@ -2145,6 +2145,10 @@ def test_append_area_defs(self, adef):
         adef.assert_called_once_with(area1.area_id, area1.description, area1.proj_id,
                                      area1.crs, area1.width, y_size, area_extent)
 
+
+class TestCreateAreaDef:
+    """Test the 'create_area_def' utility function."""
+
     @staticmethod
     def _compare_area_defs(actual, expected, use_proj4=False):
         if use_proj4:
@@ -2497,6 +2501,54 @@ def test_compute_domain(self):
         assert x_size == 5
         assert y_size == 5
 
+    @pytest.mark.parametrize(
+        (
+            "antimeridian_mode",
+            "expected_shape",
+            "expected_extents",
+            "include_proj_components",
+            "exclude_proj_components"
+        ),
+        [
+            (None, (21, 59), (164.75, 24.75, 194.25, 35.25), tuple(), ("+pm=180",)),
+            ("modify_extents", (21, 59), (164.75, 24.75, 194.25, 35.25), tuple(), ("+pm=180",)),
+            ("modify_crs", (21, 59), (164.75, 24.75, 194.25, 35.25), ("+pm=180",), tuple()),
+            ("global_extents", (21, 720), (-180.0, 24.75, 180.0, 35.25), tuple(), ("+pm=180",)),
+        ],
+    )
+    @pytest.mark.parametrize("use_dask", [False, True])
+    def test_antimeridian_mode(self,
+                               use_dask,
+                               antimeridian_mode,
+                               expected_shape,
+                               expected_extents,
+                               include_proj_components,
+                               exclude_proj_components):
+        """Test that antimeridian_mode affects the result."""
+        dyn_area = geometry.DynamicAreaDefinition('test_area', '', {'proj': 'longlat'})
+        lons, lats = _get_fake_antimeridian_lonlats(use_dask)
+        area = dyn_area.freeze(lonslats=(lons, lats), resolution=0.5, antimeridian_mode=antimeridian_mode)
+        proj_str = area.crs.to_proj4()
+
+        assert area.shape == expected_shape
+        np.testing.assert_allclose(area.area_extent, expected_extents)
+        for include_comp in include_proj_components:
+            assert include_comp in proj_str
+        for exclude_comp in exclude_proj_components:
+            assert exclude_comp not in proj_str
+
+
+def _get_fake_antimeridian_lonlats(use_dask: bool) -> tuple:
+    lon_min = 165
+    lon_max = 195
+    lons = np.arange(lon_min, lon_max, dtype=np.float64)
+    lons[lons >= 180] -= 360.0
+    lats = np.linspace(25.0, 35.0, lons.size, dtype=np.float64)
+    if use_dask:
+        lons = da.from_array(lons, chunks=lons.size // 3)
+        lats = da.from_array(lats, chunks=lons.size // 3)
+    return lons, lats
+
 
 class TestCrop(unittest.TestCase):
     """Test the area helpers."""

pyresample/utils/proj4.py

@@ -105,7 +105,12 @@ def _transform_dask_chunk(x, y, crs_from, crs_to, kwargs, transform_kwargs):
 
 
 class DaskFriendlyTransformer:
-    """Wrapper around the pyproj Transformer class that uses dask."""
+    """Wrapper around the pyproj Transformer class that uses dask.
+
+    If the provided arrays are not dask arrays, they are converted to numpy
+    arrays and pyproj will be called directly (dask is not used).
+
+    """
 
     def __init__(self, src_crs, dst_crs, **kwargs):
         """Initialize the transformer with CRS objects.
@@ -128,6 +133,12 @@ def transform(self, x, y, **kwargs):
         import dask.array as da
         crs_from = self.src_crs
         crs_to = self.dst_crs
+
+        if not hasattr(x, "compute"):
+            x = np.asarray(x)
+            y = np.asarray(y)
+            return self._transform_numpy(x, y, **kwargs)
+
         # CRS objects aren't thread-safe until pyproj 3.1+
         # convert to WKT strings to be safe
         result = da.map_blocks(_transform_dask_chunk, x, y,
@@ -139,3 +150,7 @@ def transform(self, x, y, **kwargs):
         x = result[..., 0]
         y = result[..., 1]
         return x, y
+
+    def _transform_numpy(self, x, y, **kwargs):
+        transformer = PROJTransformer.from_crs(self.src_crs, self.dst_crs, **self.kwargs)
+        return transformer.transform(x, y, **kwargs)