From 96cdf68feaf8e41e8c77798100339ffa4031d7db Mon Sep 17 00:00:00 2001 From: Philipp Kraft Date: Thu, 16 Feb 2023 12:54:54 +0100 Subject: [PATCH 1/2] Add projection.py to resolve pyproj.transform calls --- pysheds/io.py | 21 +++----- pysheds/projection.py | 120 ++++++++++++++++++++++++++++++++++++++++++ pysheds/sgrid.py | 13 ++--- pysheds/sview.py | 16 ++---- 4 files changed, 135 insertions(+), 35 deletions(-) create mode 100644 pysheds/projection.py diff --git a/pysheds/io.py b/pysheds/io.py index b91b96a..aa5f6c9 100644 --- a/pysheds/io.py +++ b/pysheds/io.py @@ -1,19 +1,14 @@ import ast import warnings import numpy as np -import pyproj import rasterio import rasterio.features from affine import Affine -from distutils.version import LooseVersion from pysheds.sview import Raster, ViewFinder, View -_OLD_PYPROJ = LooseVersion(pyproj.__version__) < LooseVersion('2.2') -_pyproj_crs = lambda Proj: Proj.crs if not _OLD_PYPROJ else Proj -_pyproj_crs_is_geographic = 'is_latlong' if _OLD_PYPROJ else 'is_geographic' -_pyproj_init = '+init=epsg:4326' if _OLD_PYPROJ else 'epsg:4326' +from . import projection -def read_ascii(data, skiprows=6, mask=None, crs=pyproj.Proj(_pyproj_init), +def read_ascii(data, skiprows=6, mask=None, crs=projection.init(), xll='lower', yll='lower', metadata={}, **kwargs): """ Reads data from an ascii file and returns a Raster. @@ -97,7 +92,7 @@ def read_raster(data, band=1, window=None, window_crs=None, mask_geometry=False, """ mask = None with rasterio.open(data, **kwargs) as f: - crs = pyproj.Proj(f.crs, preserve_units=True) + crs = projection.to_proj(f.crs) if window is None: shape = f.shape if len(f.indexes) > 1: @@ -110,13 +105,9 @@ def read_raster(data, band=1, window=None, window_crs=None, mask_geometry=False, if window_crs is not None: if window_crs.srs != crs.srs: xmin, ymin, xmax, ymax = window - if _OLD_PYPROJ: - extent = pyproj.transform(window_crs, crs, (xmin, xmax), - (ymin, ymax)) - else: - extent = pyproj.transform(window_crs, crs, (xmin, xmax), - (ymin, ymax), errcheck=True, - always_xy=True) + extent = projection.transform(window_crs, crs, (xmin, xmax), + (ymin, ymax), errcheck=True, + always_xy=True) window = (extent[0][0], extent[1][0], extent[0][1], extent[1][1]) # If window crs not specified, assume it is in raster crs ix_window = f.window(*window) diff --git a/pysheds/projection.py b/pysheds/projection.py new file mode 100644 index 0000000..058b613 --- /dev/null +++ b/pysheds/projection.py @@ -0,0 +1,120 @@ +""" +This module handles all pyproj calls of sgrid / sview. (pgrid and pview can be changed either) + +Most of the complexity comes from handling different versions. This file can be dramatically simplified +if pyproj>2.6.1 becomes a prerequisite for the next version of pysheds +""" +import pyproj +from packaging.version import Version +from typing import Any, Optional + +_pyproj_version = Version(pyproj.__version__) +_OLD_PYPROJ = _pyproj_version < Version('2.2') + + +def __transform_before_2_2( + p1: Any, + p2: Any, + x: Any, + y: Any, + z: Optional[Any] = None, + tt: Optional[Any] = None, + radians: bool = False, + errcheck: bool = False, + _: bool = True, +): + """ + A pyproj version compatibility wrapper, fixes pysheds issue #210 + + - For very old pyproj (<2.2), calls pyproj.transform without always_xy + - for old pyproj (<2.6.1), calls pyproj.transform with always_xy=True + - for recent pyproj, the deprecated call of transform is rebuild + + """ + + return pyproj.transform( + p1=p1, p2=p2, x=x, y=y, z=z, tt=tt, radians=radians, errcheck=errcheck + ) + +def __transform_before_2_6( + p1: Any, + p2: Any, + x: Any, + y: Any, + z: Optional[Any] = None, + tt: Optional[Any] = None, + radians: bool = False, + errcheck: bool = False, + always_xy: bool = True, +): + """ + A pyproj version compatibility wrapper, fixes pysheds issue #210 + + - For very old pyproj (<2.2), calls pyproj.transform without always_xy + - for old pyproj (<2.6.1), calls pyproj.transform with always_xy=True + - for recent pyproj, the deprecated call of transform is rebuild + + """ + + return pyproj.transform( + p1=p1, p2=p2, x=x, y=y, z=z, tt=tt, + radians=radians, errcheck=errcheck, always_xy=always_xy + ) + +def __transform_after_2_6( + p1: Any, + p2: Any, + x: Any, + y: Any, + z: Optional[Any] = None, + tt: Optional[Any] = None, + radians: bool = False, + errcheck: bool = False, + always_xy: bool = True, +): + """ + A pyproj version compatibility wrapper, fixes pysheds issue #210 + + - For very old pyproj (<2.2), calls pyproj.transform without always_xy + - for old pyproj (<2.6.1), calls pyproj.transform with always_xy=True + - for recent pyproj, the deprecated call of transform is rebuild + + """ + + # Code copied from pyproj.transform.transform + # see: https://github.com/pyproj4/pyproj/blob/4faef1724b044ee145832591fcd7ed2639477670/pyproj/transformer.py#L1332-L1334 + return pyproj.Transformer.from_proj(p1, p2, always_xy=always_xy).transform( + xx=x, yy=y, zz=z, tt=tt, radians=radians, errcheck=errcheck + ) + + +if _pyproj_version <= Version('2.2'): + transform = __transform_before_2_6 +elif _pyproj_version < Version('2.6.1'): + transform = __transform_before_2_6 +else: + transform = __transform_after_2_6 + +if _OLD_PYPROJ: + def init(): + """ + Init returns a newl + """ + return pyproj.Proj('+init=epsg:4326') + def to_proj(source: pyproj.Proj): + return pyproj.Proj(source, preserve_units=True) +else: + def init(): + return pyproj.Proj('epsg:4326').crs + def to_proj(source: pyproj.Proj): + return pyproj.Proj(source, preserve_units=True).crs + + +def is_geographic(crs: pyproj.CRS): + try: + return crs.is_geographic + except TypeError: + if hasattr(crs, 'is_latlong'): + return crs.is_latlong + else: + raise diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py index f0521aa..54378c5 100644 --- a/pysheds/sgrid.py +++ b/pysheds/sgrid.py @@ -1,12 +1,11 @@ import sys import ast import copy -import pyproj +# import pyproj import numpy as np import pandas as pd import geojson from affine import Affine -from distutils.version import LooseVersion try: import skimage.measure import skimage.morphology @@ -20,11 +19,6 @@ except: _HAS_RASTERIO = False -_OLD_PYPROJ = LooseVersion(pyproj.__version__) < LooseVersion('2.2') -_pyproj_crs = lambda Proj: Proj.crs if not _OLD_PYPROJ else Proj -_pyproj_crs_is_geographic = 'is_latlong' if _OLD_PYPROJ else 'is_geographic' -_pyproj_init = '+init=epsg:4326' if _OLD_PYPROJ else 'epsg:4326' - # Import input/output functions import pysheds.io @@ -34,6 +28,7 @@ # Import numba functions import pysheds._sgrid as _self +from . import projection class sGrid(): """ @@ -133,7 +128,7 @@ def defaults(self): 'affine' : Affine(1.,0.,0.,0.,1.,0.), 'shape' : (1,1), 'nodata' : 0, - 'crs' : pyproj.Proj(_pyproj_init), + 'crs' : projection.init(), } return props @@ -192,7 +187,7 @@ def extent(self): return extent def read_ascii(self, data, skiprows=6, mask=None, - crs=pyproj.Proj(_pyproj_init), xll='lower', yll='lower', + crs=projection.init(), xll='lower', yll='lower', metadata={}, **kwargs): """ Reads data from an ascii file and returns a Raster. diff --git a/pysheds/sview.py b/pysheds/sview.py index a13d13a..8cf5b35 100644 --- a/pysheds/sview.py +++ b/pysheds/sview.py @@ -1,6 +1,6 @@ import copy import numpy as np -import pyproj +from . import projection from affine import Affine from distutils.version import LooseVersion try: @@ -11,8 +11,6 @@ import pysheds._sview as _self -_OLD_PYPROJ = LooseVersion(pyproj.__version__) < LooseVersion('2.2') -_pyproj_init = '+init=epsg:4326' if _OLD_PYPROJ else 'epsg:4326' class Raster(np.ndarray): """ @@ -234,7 +232,7 @@ def to_crs(self, new_crs, **kwargs): boundary = np.vstack([top, bottom, left, right]) xi, yi = boundary[:,0], boundary[:,1] xb, yb = View.affine_transform(self.affine, xi, yi) - xb_p, yb_p = pyproj.transform(old_crs, new_crs, xb, yb, + xb_p, yb_p = projection.transform(old_crs, new_crs, xb, yb, errcheck=True, always_xy=True) x0_p = xb_p.min() if (dx > 0) else xb_p.max() y0_p = yb_p.min() if (dy > 0) else yb_p.max() @@ -331,7 +329,7 @@ class ViewFinder(): dy_dx : Tuple describing the cell size in the y and x directions. """ def __init__(self, affine=Affine(1., 0., 0., 0., 1., 0.), shape=(1,1), - nodata=0, mask=None, crs=pyproj.Proj(_pyproj_init)): + nodata=0, mask=None, crs=projection.init()): self.affine = affine self.crs = crs self.nodata = nodata @@ -411,11 +409,7 @@ def crs(self): @crs.setter def crs(self, new_crs): - try: - assert (isinstance(new_crs, pyproj.Proj)) - except: - raise TypeError('`crs` must be a `pyproj.Proj` object.') - self._crs = new_crs + self._crs = projection.to_proj(new_crs) @property def size(self): @@ -968,7 +962,7 @@ def _view_same_crs(cls, view, data, data_view, target_view, interpolation='neare @classmethod def _view_different_crs(cls, view, data, data_view, target_view, interpolation='nearest'): y, x = target_view.coords.T - xt, yt = pyproj.transform(target_view.crs, data_view.crs, x=x, y=y, + xt, yt = projection.transform(target_view.crs, data_view.crs, x=x, y=y, errcheck=True, always_xy=True) inv_affine = ~data_view.affine x_ix, y_ix = cls.affine_transform(inv_affine, xt, yt) From d26dde38f097c029a8228c9fa9c0cec85404c78b Mon Sep 17 00:00:00 2001 From: Philipp Kraft Date: Thu, 16 Mar 2023 10:20:55 +0100 Subject: [PATCH 2/2] Fixed type error, test is works --- pysheds/__init__.py | 2 +- pysheds/projection.py | 2 ++ tests/test_grid.py | 3 ++- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/pysheds/__init__.py b/pysheds/__init__.py index e19434e..56f502f 100644 --- a/pysheds/__init__.py +++ b/pysheds/__init__.py @@ -1 +1 @@ -__version__ = "0.3.3" +__version__ = "0.3.3pk" diff --git a/pysheds/projection.py b/pysheds/projection.py index 058b613..c3fe39c 100644 --- a/pysheds/projection.py +++ b/pysheds/projection.py @@ -107,6 +107,8 @@ def to_proj(source: pyproj.Proj): def init(): return pyproj.Proj('epsg:4326').crs def to_proj(source: pyproj.Proj): + if isinstance(source, pyproj.Proj): + source = source.crs return pyproj.Proj(source, preserve_units=True).crs diff --git a/tests/test_grid.py b/tests/test_grid.py index 23bc342..b49e9a2 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -5,7 +5,8 @@ from pysheds.grid import Grid from pysheds.view import Raster, ViewFinder from pysheds.rfsm import RFSM - +import pysheds +print(pysheds.__version__) current_dir = os.path.dirname(os.path.realpath(__file__)) data_dir = os.path.abspath(os.path.join(current_dir, '../data')) dir_path = os.path.join(data_dir, 'dir.asc')