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.

Sign up for GitHub

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

Update EWA pyproj Proj usage to Transformer #648

Merged
merged 3 commits into from
Feb 18, 2025
Merged

Update EWA pyproj Proj usage to Transformer #648

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5d91b3a
Update EWA from pyproj Proj to Transformer usage
djhoese Jan 30, 2025
63f12c9
Update versioneer to newest version
djhoese Feb 17, 2025
98f20a7
Fix mypy issues in versioneer
djhoese Feb 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 17 additions & 14 deletions pyresample/ewa/_ll2cr.pyx
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading