Merge pull request #447 from djhoese/bugfix-lonlat-pm180-nearest

Fix handling of lon/lat coordinates on CRS with prime meridian != 0

pyresample/_spatial_mp.py

@@ -102,43 +102,16 @@ def query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
         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 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 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.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 @@ def get_valid_index(self, geometry_def):
         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

@@ -29,11 +29,11 @@
 
 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 (
@@ -2442,13 +2442,15 @@ def get_lonlats(self, nprocs=None, data_slice=None, cache=False, dtype=None, chu
             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)
 
 
+@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():
     )
 
 
+@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 @@
 from unittest import mock
 
 import numpy as np
+from pyproj import Proj
 
 
 class TestNumpyBilinear(unittest.TestCase):
@@ -226,7 +227,6 @@ def test_solve_quadratic(self):
 
     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 @@ def test_get_input_xy(self):
 
     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 @@ def test_get_index_array(self, qnd, ria):
 
     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 @@ def test_get_four_closest_corners(self):
         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 @@ def test_get_corner(self):
         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 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 @@ def test_get_xy_from_lonlat(self):
                                       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)

pyresample/test/test_resamplers/test_nearest.py

@@ -128,6 +128,21 @@ def test_nearest_area_2d_to_area_1n(self, area_def_stere_source, data_2d_float32
         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."""