[Git][debian-gis-team/pyresample][master] 4 commits: New upstream version 1.25.1
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Fri Aug 5 10:35:28 BST 2022
Antonio Valentino pushed to branch master at Debian GIS Project / pyresample
Commits:
f72c4ee9 by Antonio Valentino at 2022-08-05T09:21:45+00:00
New upstream version 1.25.1
- - - - -
c4bfd04e by Antonio Valentino at 2022-08-05T09:22:06+00:00
Update upstream source from tag 'upstream/1.25.1'
Update to upstream version '1.25.1'
with Debian dir d0876663e2ff6e4ccb058a68c671d5bb0f1ce2ed
- - - - -
23627d37 by Antonio Valentino at 2022-08-05T09:23:10+00:00
New usptream release
- - - - -
ac797b0c by Antonio Valentino at 2022-08-05T09:27:20+00:00
Set distribution to unstable
- - - - -
15 changed files:
- CHANGELOG.md
- debian/changelog
- pyresample/_spatial_mp.py
- pyresample/bilinear/_base.py
- pyresample/bilinear/xarr.py
- pyresample/bucket/__init__.py
- pyresample/geo_filter.py
- pyresample/geometry.py
- pyresample/grid.py
- pyresample/test/conftest.py
- pyresample/test/test_bilinear.py
- pyresample/test/test_geometry.py
- pyresample/test/test_resamplers/test_nearest.py
- pyresample/utils/proj4.py
- pyresample/version.py
Changes:
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,24 @@
+## Version 1.25.0 (2022/07/29)
+
+### Issues Closed
+
+* [Issue 428](https://github.com/pytroll/pyresample/issues/428) - Add more flexible antimeridian handling to DynamicAreaDefinition ([PR 431](https://github.com/pytroll/pyresample/pull/431) by [@djhoese](https://github.com/djhoese))
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 446](https://github.com/pytroll/pyresample/pull/446) - Fix incorrect extents for DynamicAreaDefinition with 'modify_crs' antimeridian mode
+
+#### Features added
+
+* [PR 431](https://github.com/pytroll/pyresample/pull/431) - Add 'antimeridian_mode' to DynamicAreaDefinition ([428](https://github.com/pytroll/pyresample/issues/428))
+
+In this release 2 pull requests were closed.
+
+
## Version 1.24.1 (2022/07/06)
### Pull Requests Merged
=====================================
debian/changelog
=====================================
@@ -1,9 +1,10 @@
-pyresample (1.24.1-2) UNRELEASED; urgency=medium
+pyresample (1.25.1-1) unstable; urgency=medium
+ * New upstream release.
* debian/copyright:
- fix typo.
- -- Antonio Valentino <antonio.valentino at tiscali.it> Sat, 09 Jul 2022 16:26:51 +0000
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Fri, 05 Aug 2022 09:26:55 +0000
pyresample (1.24.1-1) unstable; urgency=medium
=====================================
pyresample/_spatial_mp.py
=====================================
@@ -102,43 +102,16 @@ class cKDTree_MP(object):
return _d.copy(), _i.copy()
-class BaseProj(pyproj.Proj):
- """Helper class for easier backwards compatibility."""
-
- def __init__(self, projparams=None, preserve_units=True, **kwargs):
- # have to have this because pyproj uses __new__
- # subclasses would fail when calling __init__ otherwise
- super(BaseProj, self).__init__(projparams=projparams,
- preserve_units=preserve_units,
- **kwargs)
-
-
-class Proj(BaseProj):
- """Helper class to skip transforming lon/lat projection coordinates."""
-
- def __call__(self, data1, data2, inverse=False, radians=False,
- errcheck=False, nprocs=1):
- """Transform coordinates to coordinate system except for geographic coordinate systems."""
- if self.crs.is_geographic:
- return data1, data2
- return super(Proj, self).__call__(data1, data2, inverse=inverse,
- radians=radians, errcheck=errcheck)
-
-
-class Proj_MP(BaseProj):
+class Proj_MP:
"""Multi-processing version of the pyproj Proj class."""
def __init__(self, *args, **kwargs):
self._args = args
self._kwargs = kwargs
- super(Proj_MP, self).__init__(*args, **kwargs)
def __call__(self, data1, data2, inverse=False, radians=False,
errcheck=False, nprocs=2, chunk=None, schedule='guided'):
"""Transform coordinates to coordinates in the current coordinate system."""
- if self.crs.is_geographic:
- return data1, data2
-
grid_shape = data1.shape
n = data1.size
=====================================
pyresample/bilinear/_base.py
=====================================
@@ -34,9 +34,9 @@ import warnings
import numpy as np
from pykdtree.kdtree import KDTree
+from pyproj import Proj
from pyresample import data_reduce, geometry
-from pyresample._spatial_mp import Proj
from ..future.resamplers._transform_utils import lonlat2xyz
=====================================
pyresample/bilinear/xarr.py
=====================================
@@ -30,10 +30,10 @@ import dask.array as da
import numpy as np
import zarr
from dask import delayed
+from pyproj import Proj
from xarray import DataArray, Dataset
from pyresample import CHUNK_SIZE
-from pyresample._spatial_mp import Proj
from pyresample.bilinear._base import (
BilinearBase,
array_slice_for_multiple_arrays,
=====================================
pyresample/bucket/__init__.py
=====================================
@@ -24,8 +24,7 @@ import dask
import dask.array as da
import numpy as np
import xarray as xr
-
-from pyresample._spatial_mp import Proj
+from pyproj import Proj
LOG = logging.getLogger(__name__)
=====================================
pyresample/geo_filter.py
=====================================
@@ -18,6 +18,7 @@
"""Filters based on geolocation validity."""
import numpy as np
+from pyproj import Proj
from . import _spatial_mp, geometry
@@ -60,12 +61,14 @@ class GridFilter(object):
lats = geometry_def.lats[:]
# Get projection coords
+ proj_kwargs = {}
if self.nprocs > 1:
proj = _spatial_mp.Proj_MP(**self.area_def.proj_dict)
+ proj_kwargs["nprocs"] = self.nprocs
else:
- proj = _spatial_mp.Proj(**self.area_def.proj_dict)
+ proj = Proj(**self.area_def.proj_dict)
- x_coord, y_coord = proj(lons, lats, nprocs=self.nprocs)
+ x_coord, y_coord = proj(lons, lats, **proj_kwargs)
# Find array indices of coordinates
target_x = ((x_coord / self.area_def.pixel_size_x) +
=====================================
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,14 +25,15 @@ from collections import OrderedDict
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 import Geod, Proj, transform
+from pyproj.aoi import AreaOfUse
from pyresample import CHUNK_SIZE
-from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj, Proj_MP
+from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP
from pyresample.area_config import create_area_def
from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary
from pyresample.utils import (
@@ -1017,26 +1019,59 @@ class DynamicAreaDefinition(object):
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))
@@ -1047,7 +1082,33 @@ class DynamicAreaDefinition(object):
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
@@ -1061,6 +1122,31 @@ class DynamicAreaDefinition(object):
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.
@@ -1082,76 +1168,59 @@ class DynamicAreaDefinition(object):
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)
+ if antimeridian_mode == "modify_crs":
+ xmin -= 180
+ xmax -= 180
+ return xmin, xmax
+
def _invproj(data_x, data_y, proj_dict):
"""Perform inverse projection."""
@@ -2373,13 +2442,15 @@ class AreaDefinition(_ProjectionDefinition):
lons, lats = res[0], res[1]
return lons, lats
+ proj_kwargs = {}
if nprocs > 1:
target_proj = Proj_MP(self.crs)
+ proj_kwargs["nprocs"] = nprocs
else:
target_proj = Proj(self.crs)
# Get corresponding longitude and latitude values
- lons, lats = target_proj(target_x, target_y, inverse=True, nprocs=nprocs)
+ lons, lats = target_proj(target_x, target_y, inverse=True, **proj_kwargs)
lons = np.asanyarray(lons, dtype=dtype)
lats = np.asanyarray(lats, dtype=dtype)
=====================================
pyresample/grid.py
=====================================
@@ -20,6 +20,7 @@
from __future__ import absolute_import
import numpy as np
+from pyproj import Proj
from pyresample import _spatial_mp, geometry
@@ -107,13 +108,15 @@ def get_linesample(lons, lats, source_area_def, nprocs=1):
Arrays for resampling area by array indexing
"""
# Proj.4 definition of source area projection
+ proj_kwargs = {}
if nprocs > 1:
source_proj = _spatial_mp.Proj_MP(**source_area_def.proj_dict)
+ proj_kwargs["nprocs"] = nprocs
else:
- source_proj = _spatial_mp.Proj(**source_area_def.proj_dict)
+ source_proj = Proj(**source_area_def.proj_dict)
# get cartesian projection values from longitude and latitude
- source_x, source_y = source_proj(lons, lats, nprocs=nprocs)
+ source_x, source_y = source_proj(lons, lats, **proj_kwargs)
# Find corresponding pixels (element by element conversion of ndarrays)
source_pixel_x = (source_area_def.pixel_offset_x +
=====================================
pyresample/test/conftest.py
=====================================
@@ -51,7 +51,7 @@ def _euro_lonlats_dask():
def _antimeridian_lonlats():
lons = create_test_longitude(172.0, 190.0, SRC_SWATH_2D_SHAPE)
- lons[lons > 180.0] = lons - 360.0
+ lons[lons > 180.0] -= 360.0
lats = create_test_latitude(25.0, 33.0, SRC_SWATH_2D_SHAPE)
return lons, lats
@@ -99,6 +99,19 @@ def swath_def_2d_numpy_antimeridian():
return SwathDefinition(lons, lats)
+ at pytest.fixture(scope="session")
+def swath_def_2d_xarray_dask_antimeridian():
+ """Create a SwathDefinition with DataArrays(dask) arrays (200, 1500) over the antimeridian.
+
+ Longitude values go from positive values to negative values as they cross -180/180.
+
+ """
+ lons, lats = _antimeridian_lonlats()
+ lons = xr.DataArray(lons, dims=("y", "x"))
+ lats = xr.DataArray(lats, dims=("y", "x"))
+ return SwathDefinition(lons, lats)
+
+
@pytest.fixture(scope="session")
def area_def_lcc_conus_1km():
"""Create an AreaDefinition with an LCC projection over CONUS (1500, 2000)."""
@@ -150,6 +163,22 @@ def area_def_stere_target():
)
+ at pytest.fixture(scope="session")
+def area_def_lonlat_pm180_target():
+ """Create an AreaDefinition with a geographic lon/lat projection with prime meridian at 180 (800, 850)."""
+ return AreaDefinition(
+ 'lonlat_pm180', '', '',
+ {
+ 'proj': 'longlat',
+ 'pm': '180.0',
+ 'datum': 'WGS84',
+ 'no_defs': None,
+ },
+ DST_AREA_SHAPE[1], DST_AREA_SHAPE[0],
+ [-20.0, 20.0, 20.0, 35.0]
+ )
+
+
@pytest.fixture(scope="session")
def coord_def_2d_float32_dask():
"""Create a 2D CoordinateDefinition of dask arrays (4, 3)."""
=====================================
pyresample/test/test_bilinear.py
=====================================
@@ -27,6 +27,7 @@ import unittest
from unittest import mock
import numpy as np
+from pyproj import Proj
class TestNumpyBilinear(unittest.TestCase):
@@ -226,7 +227,6 @@ class TestNumpyBilinear(unittest.TestCase):
def test_get_input_xy(self):
"""Test calculation of input xy-coordinates."""
- from pyresample._spatial_mp import Proj
from pyresample.bilinear._base import _get_input_xy
proj = Proj(self.target_def.proj_str)
@@ -237,7 +237,6 @@ class TestNumpyBilinear(unittest.TestCase):
def test_get_four_closest_corners(self):
"""Test calculation of bounding corners."""
- from pyresample._spatial_mp import Proj
from pyresample.bilinear._base import (
_get_four_closest_corners,
_get_input_xy,
@@ -842,7 +841,6 @@ class TestXarrayBilinear(unittest.TestCase):
def test_get_input_xy(self):
"""Test computation of input X and Y coordinates in target proj."""
- from pyresample._spatial_mp import Proj
from pyresample.bilinear.xarr import _get_input_xy
proj = Proj(self.target_def.proj_str)
@@ -860,7 +858,6 @@ class TestXarrayBilinear(unittest.TestCase):
import dask.array as da
from pyresample import CHUNK_SIZE
- from pyresample._spatial_mp import Proj
from pyresample.bilinear._base import _get_four_closest_corners
from pyresample.bilinear.xarr import _get_input_xy
@@ -891,7 +888,6 @@ class TestXarrayBilinear(unittest.TestCase):
import dask.array as da
from pyresample import CHUNK_SIZE
- from pyresample._spatial_mp import Proj
from pyresample.bilinear._base import _get_corner, _get_input_xy
proj = Proj(self.target_def.proj_str)
=====================================
pyresample/test/test_geometry.py
=====================================
@@ -26,7 +26,7 @@ import dask.array as da
import numpy as np
import pytest
import xarray as xr
-from pyproj import CRS
+from pyproj import CRS, Proj
from pyresample import geo_filter, geometry, parse_area_file
from pyresample.geometry import (
@@ -1134,7 +1134,6 @@ class Test(unittest.TestCase):
proj_dict,
x_size, y_size,
area_extent)
- from pyresample._spatial_mp import Proj
p__ = Proj(proj_dict)
lon_ul, lat_ul = p__(1000000, 50000, inverse=True)
lon_ur, lat_ur = p__(1050000, 50000, inverse=True)
@@ -2172,6 +2171,10 @@ class TestStackedAreaDefinition:
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:
@@ -2524,6 +2527,54 @@ class TestDynamicAreaDefinition:
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 - 180.0, 24.75, 194.25 - 180.0, 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/test/test_resamplers/test_nearest.py
=====================================
@@ -128,6 +128,21 @@ class TestNearestNeighborResampler:
assert cross_sum == expected
assert res.shape == resampler.target_geo_def.shape
+ def test_nearest_swath_2d_to_area_1n_pm180(self, swath_def_2d_xarray_dask_antimeridian, data_2d_float32_xarray_dask,
+ area_def_lonlat_pm180_target):
+ """Test 2D swath definition to 2D area definition; 1 neighbor; output prime meridian at 180 degrees."""
+ resampler = KDTreeNearestXarrayResampler(
+ swath_def_2d_xarray_dask_antimeridian, area_def_lonlat_pm180_target)
+ res = resampler.resample(data_2d_float32_xarray_dask, radius_of_influence=50000)
+ assert isinstance(res, xr.DataArray)
+ assert isinstance(res.data, da.Array)
+ _check_common_metadata(res, isinstance(area_def_lonlat_pm180_target, AreaDefinition))
+ res = res.values
+ cross_sum = float(np.nansum(res))
+ expected = 115591.0
+ assert cross_sum == expected
+ assert res.shape == resampler.target_geo_def.shape
+
def test_nearest_area_2d_to_area_1n_no_roi(self, area_def_stere_source, data_2d_float32_xarray_dask,
area_def_stere_target):
"""Test 2D area definition to 2D area definition; 1 neighbor, no radius of influence."""
=====================================
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 @@ class DaskFriendlyTransformer:
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 @@ class DaskFriendlyTransformer:
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)
=====================================
pyresample/version.py
=====================================
@@ -25,9 +25,9 @@ def get_keywords():
# setup.py/versioneer.py will grep for the variable names, so they must
# each be defined on a line of their own. _version.py will just call
# get_keywords().
- git_refnames = " (HEAD -> main, tag: v1.24.1)"
- git_full = "c930465a94f63c9d6a9db1b0c3db033d0369a396"
- git_date = "2022-07-06 19:04:53 -0500"
+ git_refnames = " (tag: v1.25.1)"
+ git_full = "164b714a9d00a2b5e103a2214400248047dc8424"
+ git_date = "2022-08-02 07:46:40 -0500"
keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
return keywords
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyresample/-/compare/2f17b9dac698094405f3cbc0bc401e82036d5611...ac797b0c33cd6cf1a15a42a00739070462a43bca
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyresample/-/compare/2f17b9dac698094405f3cbc0bc401e82036d5611...ac797b0c33cd6cf1a15a42a00739070462a43bca
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20220805/6f52d101/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list