diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 460796d77..e5c81a19a 100644 --- a/pyresample/geometry.py +++ b/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 . """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.""" diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index 4a39574b5..a359e851c 100644 --- a/pyresample/test/test_geometry.py +++ b/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.""" diff --git a/pyresample/utils/proj4.py b/pyresample/utils/proj4.py index 3d4b07eee..b512b0112 100644 --- a/pyresample/utils/proj4.py +++ b/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)