diff --git a/pyresample/ewa/_ll2cr.pyx b/pyresample/ewa/_ll2cr.pyx index c1810f7a..eb1ff383 100644 --- a/pyresample/ewa/_ll2cr.pyx +++ b/pyresample/ewa/_ll2cr.pyx @@ -23,11 +23,14 @@ __docformat__ = "restructuredtext en" import numpy -from pyproj import Proj +from pyproj import Transformer +from pyproj.enums import TransformDirection cimport cython cimport numpy +from pyresample.utils.proj4 import get_geodetic_crs_with_no_datum_shift + numpy.import_array() # column and rows can only be doubles for now until the PROJ.4 is linked directly so float->double casting can be done @@ -40,18 +43,18 @@ cdef extern from "numpy/npy_math.h": bint npy_isnan(double x) -def projection_circumference(p): +def projection_circumference(t): """Return the projection circumference if the projection is cylindrical, None otherwise. Projections that are not cylindrical and centered on the globes axis can not easily have data cross the antimeridian of the projection. """ - lon0, lat0 = p(0, 0, inverse=True) + lon0, lat0 = t.transform(0, 0, direction=TransformDirection.INVERSE) lon1 = lon0 + 180.0 lat1 = lat0 + 5.0 - x0, y0 = p(lon0, lat0) # should result in zero or near zero - x1, y1 = p(lon1, lat0) - x2, y2 = p(lon1, lat1) + x0, y0 = t.transform(lon0, lat0) # should result in zero or near zero + x1, y1 = t.transform(lon1, lat0) + x2, y2 = t.transform(lon1, lat1) if y0 != y1 or x1 != x2: return 0.0 return abs(x1 - x0) * 2 @@ -61,7 +64,7 @@ def projection_circumference(p): @cython.wraparound(False) @cython.cdivision(True) def ll2cr_dynamic(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtype, ndim=2] lat_arr, - cr_dtype fill_in, str proj4_definition, + cr_dtype fill_in, object src_crs, object dst_crs, double cell_width, double cell_height, width=None, height=None, origin_x=None, origin_y=None): @@ -98,13 +101,13 @@ def ll2cr_dynamic(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtyp of limitations in pyproj. """ # pure python stuff for now - p = Proj(proj4_definition) + t = Transformer.from_crs(src_crs, dst_crs, always_xy=True) # Pyproj currently makes a copy so we don't have to do anything special here - cdef tuple projected_tuple = p(lon_arr, lat_arr) + cdef tuple projected_tuple = t.transform(lon_arr, lat_arr) cdef cr_dtype[:, ::1] rows_out = projected_tuple[1] cdef cr_dtype[:, ::1] cols_out = projected_tuple[0] - cdef double proj_circum = projection_circumference(p) + cdef double proj_circum = projection_circumference(t) cdef unsigned int w cdef unsigned int h cdef double ox @@ -203,7 +206,7 @@ def ll2cr_dynamic(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtyp @cython.wraparound(False) @cython.cdivision(True) def ll2cr_static(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtype, ndim=2] lat_arr, - cr_dtype fill_in, str proj4_definition, + cr_dtype fill_in, object src_crs, object dst_crs, double cell_width, double cell_height, unsigned int width, unsigned int height, double origin_x, double origin_y): @@ -229,11 +232,11 @@ def ll2cr_static(numpy.ndarray[cr_dtype, ndim=2] lon_arr, numpy.ndarray[cr_dtype of limitations in pyproj. """ # TODO: Rewrite so it is no GIL - # pure python stuff for now - p = Proj(proj4_definition) + # pure python stuff for now + t = Transformer.from_crs(src_crs, dst_crs, always_xy=True) # Pyproj currently makes a copy so we don't have to do anything special here - cdef tuple projected_tuple = p(lon_arr, lat_arr) + cdef tuple projected_tuple = t.transform(lon_arr, lat_arr) cdef cr_dtype[:, ::1] rows_out = projected_tuple[1] cdef cr_dtype[:, ::1] cols_out = projected_tuple[0] cdef cr_dtype[:, ::1] lons_view = lon_arr diff --git a/pyresample/ewa/ewa.py b/pyresample/ewa/ewa.py index 1ffcb763..9b81b491 100644 --- a/pyresample/ewa/ewa.py +++ b/pyresample/ewa/ewa.py @@ -76,7 +76,7 @@ def ll2cr(swath_def, area_def, fill=np.nan, copy=True): ox = area_def.area_extent[0] + cw / 2. oy = area_def.area_extent[3] + ch / 2. swath_points_in_grid = _ll2cr.ll2cr_static(lons, lats, fill, - p, cw, ch, w, h, ox, oy) + swath_def.crs, area_def.crs, cw, ch, w, h, ox, oy) return swath_points_in_grid, lons, lats diff --git a/pyresample/test/test_ewa_ll2cr.py b/pyresample/test/test_ewa_ll2cr.py index ace95c48..b74e798b 100644 --- a/pyresample/test/test_ewa_ll2cr.py +++ b/pyresample/test/test_ewa_ll2cr.py @@ -21,6 +21,7 @@ import unittest import numpy as np +from pyproj import CRS from pyresample.test.utils import create_test_latitude, create_test_longitude @@ -70,14 +71,15 @@ def test_lcc_basic1(self): lat_arr = create_test_latitude(18.0, 40.0, (50, 100), dtype=np.float64) grid_info = static_lcc.copy() fill_in = np.nan - proj_str = grid_info["proj4_definition"] + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_proj4(grid_info["proj4_definition"]) cw = grid_info["cell_width"] ch = grid_info["cell_height"] ox = grid_info["origin_x"] oy = grid_info["origin_y"] w = grid_info["width"] h = grid_info["height"] - points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, proj_str, + points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, src_crs, dst_crs, cw, ch, w, h, ox, oy) self.assertEqual(points_in_grid, lon_arr.size, "all these test points should fall in this grid") @@ -92,14 +94,15 @@ def test_geo_antimeridian(self): grid_info = static_geo_whole_earth.copy() fill_in = np.nan - proj_str = grid_info['proj4_definition'] + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_proj4(grid_info["proj4_definition"]) cw = grid_info['cell_width'] ch = grid_info['cell_height'] ox = grid_info['origin_x'] oy = grid_info['origin_y'] w = grid_info['width'] h = grid_info['height'] - points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, proj_str, + points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, src_crs, dst_crs, cw, ch, w, h, ox, oy) self.assertEqual(points_in_grid, lon_arr.size, 'all these test points should fall in this grid') @@ -110,14 +113,15 @@ def test_lcc_fail1(self): lat_arr = create_test_latitude(18.0, 40.0, (50, 100), dtype=np.float64) grid_info = static_lcc.copy() fill_in = np.nan - proj_str = grid_info["proj4_definition"] + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_proj4(grid_info["proj4_definition"]) cw = grid_info["cell_width"] ch = grid_info["cell_height"] ox = grid_info["origin_x"] oy = grid_info["origin_y"] w = grid_info["width"] h = grid_info["height"] - points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, proj_str, + points_in_grid = _ll2cr.ll2cr_static(lon_arr, lat_arr, fill_in, src_crs, dst_crs, cw, ch, w, h, ox, oy) self.assertEqual(points_in_grid, 0, "none of these test points should fall in this grid") @@ -131,14 +135,16 @@ def test_latlong_basic1(self): lat_arr = create_test_latitude(15.0, 30.0, (50, 100), dtype=np.float64) grid_info = dynamic_wgs84.copy() fill_in = np.nan - proj_str = grid_info["proj4_definition"] + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_proj4(grid_info["proj4_definition"]) cw = grid_info["cell_width"] ch = grid_info["cell_height"] ox = grid_info["origin_x"] oy = grid_info["origin_y"] w = grid_info["width"] h = grid_info["height"] - points_in_grid, lon_res, lat_res, ox, oy, w, h = _ll2cr.ll2cr_dynamic(lon_arr, lat_arr, fill_in, proj_str, + points_in_grid, lon_res, lat_res, ox, oy, w, h = _ll2cr.ll2cr_dynamic(lon_arr, lat_arr, fill_in, + src_crs, dst_crs, cw, ch, w, h, ox, oy) self.assertEqual(points_in_grid, lon_arr.size, "all points should be contained in a dynamic grid") self.assertIs(lon_arr, lon_res) @@ -152,14 +158,16 @@ def test_latlong_basic2(self): lat_arr = create_test_latitude(15.0, 30.0, (50, 100), twist_factor=-0.1, dtype=np.float64) grid_info = dynamic_wgs84.copy() fill_in = np.nan - proj_str = grid_info["proj4_definition"] + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_proj4(grid_info["proj4_definition"]) cw = grid_info["cell_width"] ch = grid_info["cell_height"] ox = grid_info["origin_x"] oy = grid_info["origin_y"] w = grid_info["width"] h = grid_info["height"] - points_in_grid, lon_res, lat_res, ox, oy, w, h = _ll2cr.ll2cr_dynamic(lon_arr, lat_arr, fill_in, proj_str, + points_in_grid, lon_res, lat_res, ox, oy, w, h = _ll2cr.ll2cr_dynamic(lon_arr, lat_arr, fill_in, + src_crs, dst_crs, cw, ch, w, h, ox, oy) self.assertEqual(points_in_grid, lon_arr.size, "all points should be contained in a dynamic grid") self.assertIs(lon_arr, lon_res) @@ -173,14 +181,16 @@ def test_latlong_dateline1(self): lat_arr = create_test_latitude(15.0, 30.0, (50, 100), twist_factor=-0.1, dtype=np.float64) grid_info = dynamic_wgs84.copy() fill_in = np.nan - proj_str = grid_info["proj4_definition"] + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_proj4(grid_info["proj4_definition"]) cw = grid_info["cell_width"] ch = grid_info["cell_height"] ox = grid_info["origin_x"] oy = grid_info["origin_y"] w = grid_info["width"] h = grid_info["height"] - points_in_grid, lon_res, lat_res, ox, oy, w, h = _ll2cr.ll2cr_dynamic(lon_arr, lat_arr, fill_in, proj_str, + points_in_grid, lon_res, lat_res, ox, oy, w, h = _ll2cr.ll2cr_dynamic(lon_arr, lat_arr, fill_in, + src_crs, dst_crs, cw, ch, w, h, ox, oy) self.assertEqual(points_in_grid, lon_arr.size, "all points should be contained in a dynamic grid") self.assertIs(lon_arr, lon_res) diff --git a/versioneer.py b/versioneer.py index a480f5eb..6213dfd8 100644 --- a/versioneer.py +++ b/versioneer.py @@ -1898,7 +1898,7 @@ def get_cmdclass(cmdclass: Optional[Dict[str, Any]] = None): class cmd_version(Command): description = "report generated version string" - user_options: List[Tuple[str, str, str]] = [] + user_options: List[Tuple[str, str, str]] = [] # type: ignore boolean_options: List[str] = [] def initialize_options(self) -> None: