Skip to content

Commit

Permalink
Merge pull request #648 from djhoese/udpate-ewa-pyproj-transformer
Browse files Browse the repository at this point in the history
  • Loading branch information
@djhoese
djhoese authored Feb 18, 2025
2 parents 1f974b5 + 98f20a7 commit 873a511
Show file tree
Hide file tree
Showing 4 changed files with 41 additions and 28 deletions.
31 changes: 17 additions & 14 deletions pyresample/ewa/_ll2cr.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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):
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion pyresample/ewa/ewa.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down
34 changes: 22 additions & 12 deletions pyresample/test/test_ewa_ll2cr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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")

Expand All @@ -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')
Expand All @@ -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")

Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion versioneer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down

0 comments on commit 873a511

Please sign in to comment.