Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add 'antimeridian_mode' to DynamicAreaDefinition #431

Merged
merged 9 commits into from
Jul 28, 2022
186 changes: 132 additions & 54 deletions pyresample/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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
mraspaud marked this conversation as resolved.
Show resolved Hide resolved
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
Expand All @@ -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.
gerritholl marked this conversation as resolved.
Show resolved Hide resolved
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_projection": 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.
Expand All @@ -1047,76 +1133,68 @@ 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_projection":
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:
if hasattr(xarr, "compute"):
xmin, xmax = self._wrap_x_corners_dask(xarr)
else:
xmin, xmax = self._wrap_x_corners_numpy(xarr)
djhoese marked this conversation as resolved.
Show resolved Hide resolved
return xmin, xmax

def _wrap_x_corners_numpy(self, xarr):
xmin = np.nanmin(xarr[xarr >= 0])
xmax = np.nanmax(xarr[xarr < 0]) + 360
return xmin, xmax

def _wrap_x_corners_dask(self, xarr):
xarr_pos = da.where(xarr >= 0, xarr, np.nan)
xarr_neg = da.where(xarr < 0, xarr, np.nan)
new_xmin = np.nanmin(xarr_pos)
new_xmax = np.nanmax(xarr_neg) + 360
xmin, xmax = da.compute(new_xmin, new_xmax)
return xmin, xmax


def _invproj(data_x, data_y, proj_dict):
"""Perform inverse projection."""
Expand Down
52 changes: 52 additions & 0 deletions pyresample/test/test_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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_projection", (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."""
Expand Down
17 changes: 16 additions & 1 deletion pyresample/utils/proj4.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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,
Expand All @@ -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)