Update to pyproj 2 (#92)

* Rename projection to crs

Follows pyproj in nomenclature. See https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 .

* environment: Remove channel pinning

Channel pinning has been superseed by strict channel_priority as
proposed at https://conda-forge.org/docs/user/tipsandtricks.html.

* gis: Add grid_cell_areas function to compute areas of grid cells

* cutout: Fix forgotten conversion

* gis: Improve grid_cell_areas

* remove area calculation due to geopandas implementation

* update release notes

* gis.py: revise imports

Co-authored-by: Fabian <fab.hof@gmx.de>

RELEASE_NOTES.rst

@@ -30,6 +30,7 @@ Major changes
     which combines and deprecates `grid_cells()` and `grid_coordinates()`
 * The order of coordinates (indices) for `Cutouts` changed: `x` and `y` (e.g. longitude and latitude) are now
     both ascending (before: `x` ascending and `y` descending).
+* Following the lead of geopandas, pyproj, cartopy and rasterio, atlite now uses Coordinate Reference System (`CRS`) instead of the old fashioned projection strings. 
 
 New features
 ------------

atlite/convert.py

@@ -37,7 +37,7 @@
 
 def convert_and_aggregate(cutout, convert_func, windows=None, matrix=None,
                           index=None, layout=None, shapes=None,
-                          shapes_proj='latlong', per_unit=False,
+                          shapes_crs=4326, per_unit=False,
                           return_capacity=False, capacity_factor=False,
                           show_progress=True, **convert_kwds):
     """
@@ -59,10 +59,9 @@ def convert_and_aggregate(cutout, convert_func, windows=None, matrix=None,
     shapes : list or pd.Series of shapely.geometry.Polygon
         If given, matrix is constructed as indicatormatrix of the polygons, its
         index determines the bus index on the time-series.
-    shapes_proj : str or pyproj.Proj
-        Defaults to 'latlong'. If different to the map projection of the
-        cutout, the shapes are reprojected using pyproj.transform to match
-        cutout.projection (defaults to 'latlong').
+    shapes_crs : pyproj.CRS or compatible
+        If different to the map crs of the cutout, the shapes are
+        transformed to match cutout.crs (defaults to EPSG:4326).
     per_unit : boolean
         Returns the time-series in per-unit units, instead of in MW (defaults
         to False).
@@ -110,7 +109,7 @@ def convert_and_aggregate(cutout, convert_func, windows=None, matrix=None,
             geoseries_like = (pd.Series, gpd.GeoDataFrame, gpd.GeoSeries)
             if isinstance(shapes, geoseries_like) and index is None:
                 index = shapes.index
-            matrix = cutout.indicatormatrix(shapes, shapes_proj)
+            matrix = cutout.indicatormatrix(shapes, shapes_crs)
 
         if matrix is not None:
             matrix = csr_matrix(matrix)

atlite/cutout.py

@@ -20,13 +20,14 @@
 import xarray as xr
 import pandas as pd
 import numpy as np
+import geopandas as gpd
 from tempfile import mktemp
 from numpy import atleast_1d
 from warnings import warn
 from shapely.geometry import box
 from pathlib import Path
-import geopandas as gpd
-import pyproj
+from pyproj import CRS
+
 
 from .utils import CachedAttribute
 from .data import cutout_prepare, available_features
@@ -187,10 +188,10 @@ def __init__(self, path, **cutoutparams):
                      **storable_chunks, **cutoutparams}
             data = xr.Dataset(coords=coords, attrs=attrs)
 
-        # Check projections
+        # Check compatibility of CRS
         modules = atleast_1d(data.attrs.get('module'))
-        projection = set(datamodules[m].projection for m in modules)
-        assert len(projection) == 1, f'Projections of {module} not compatible'
+        crs = set(CRS(datamodules[m].crs) for m in modules)
+        assert len(crs) == 1, f'CRS of {module} not compatible'
 
         self.path = path
         self.data = data
@@ -205,8 +206,8 @@ def module(self):
         return self.data.attrs.get('module')
 
     @property
-    def projection(self):
-        return datamodules[atleast_1d(self.module)[0]].projection
+    def crs(self):
+        return CRS(datamodules[atleast_1d(self.module)[0]].crs)
 
     @property
     def available_features(self):
@@ -284,11 +285,8 @@ def grid(self):
         coords = np.asarray((np.ravel(xs), np.ravel(ys))).T
         span = (coords[self.shape[1] + 1] - coords[0]) / 2
         cells = [box(*c) for c in np.hstack((coords - span, coords + span))]
-        # TODO!: crs should be specific by the module (atm all module have the
-        # same crs)
-        crs = pyproj.CRS.from_epsg(4326)
         return gpd.GeoDataFrame({'x': coords[:, 0], 'y': coords[:, 1],
-                                 'geometry': cells,}, crs=crs)
+                                 'geometry': cells,}, crs=self.crs)
 
 
     def sel(self, path=None, bounds=None, buffer=0, **kwargs):
@@ -322,9 +320,8 @@ def __repr__(self):
                         list(self.prepared_features.index.unique('feature'))))
 
 
-    def indicatormatrix(self, shapes, shapes_proj='latlong'):
-        return compute_indicatormatrix(self.grid, shapes,
-                                       self.projection, shapes_proj)
+    def indicatormatrix(self, shapes, shapes_crs=4326):
+        return compute_indicatormatrix(self.grid, shapes, self.crs, shapes_crs)
 
     # Preparation functions
 

atlite/datasets/cordex.py

@@ -23,12 +23,13 @@
 import os
 import glob
 
-from ..gis import RotProj
-
-# Model and Projection Settings
+# Model and CRS Settings
 model = 'MPI-M-MPI-ESM-LR'
-projection = RotProj(dict(proj='ob_tran', o_proj='latlong', lon_0=180,
-                          o_lon_p=-162, o_lat_p=39.25))
+
+crs = 4326 # TODO
+# something like the following is correct
+# crs = pyproj.crs.DerivedGeographicCRS(4326, pcrs.coordinate_operation.RotatedLatitudeLongitudeConversion(??))
+# RotProj(dict(proj='ob_tran', o_proj='latlong', lon_0=180, o_lon_p=-162, o_lat_p=39.25))
 
 
 def rename_and_clean_coords(ds):

atlite/datasets/era5.py

@@ -33,8 +33,8 @@ def nullcontext():
 
 logger = logging.getLogger(__name__)
 
-# Model and Projection Settings
-projection = 'latlong'
+# Model and CRS Settings
+crs = 4326
 
 features = {
     'height': ['height'],

atlite/datasets/gebco.py

@@ -14,7 +14,7 @@
 import logging
 logger = logging.getLogger(__name__)
 
-projection = 'latlong'
+crs = 4326
 features = {'height': ['height']}
 
 

atlite/datasets/ncep.py

@@ -24,7 +24,7 @@
 import shutil
 
 engine = 'pynio'
-projection = 'latlong'
+crs = 4326
 
 
 def convert_lons_lats_ncep(ds, xs, ys):

atlite/datasets/sarah.py

@@ -21,8 +21,8 @@
 logger = logging.getLogger(__name__)
 
 
-# Model, Projection and Resolution Settings
-projection = 'latlong'
+# Model, CRS and Resolution Settings
+crs = 4326
 dx = 0.05
 dy = 0.05
 dt = '30min'

atlite/gis.py

@@ -18,15 +18,12 @@
 import scipy.sparse
 from collections import OrderedDict
 from warnings import warn
-from functools import partial
-import pyproj
+from pyproj import CRS, Transformer
 import geopandas as gpd
 from shapely.ops import transform
 import rasterio as rio
 import rasterio.warp
 from shapely.strtree import STRtree
-from rasterio.warp import Resampling
-from rasterio.crs import CRS
 
 import logging
 logger = logging.getLogger(__name__)
@@ -75,57 +72,16 @@ def spdiag(v):
     return sp.sparse.csr_matrix((v, inds[:-1], inds), (N, N))
 
 
-class RotProj(pyproj.Proj):
-    def __call__(self, x, y, inverse=False, **kw):
-        if inverse:
-            gx, gy = super(RotProj, self).__call__(x, y,
-                                                   inverse=False, **kw)
-            return np.rad2deg(gx), np.rad2deg(gy)
-        else:
-            return super(RotProj, self).__call__(np.deg2rad(x),
-                                                 np.deg2rad(y),
-                                                 inverse=True, **kw)
-
-
-def as_projection(p):
-    if isinstance(p, pyproj.Proj):
-        return p
-    elif isinstance(p, str):
-        return pyproj.Proj(dict(proj=p))
-    else:
-        return pyproj.Proj(p)
-
-
-def reproject_shapes(shapes, p1, p2):
+def reproject_shapes(shapes, crs1, crs2):
     """
-    Project a collection of `shapes` from one projection `p1` to
-    another projection `p2`
-
-    Projections can be given as strings or instances of pyproj.Proj.
-    Special care is taken for the case where the final projection is
-    of type rotated pole as handled by RotProj.
+    Project a collection of `shapes` from one crs `crs1` to
+    another crs `crs2`
     """
 
-    if p1 == p2:
-        return shapes
-
-    if isinstance(p1, RotProj):
-        if p2 == 'latlong':
-            def reproject_points(x, y):
-                return p1(x, y, inverse=True)
-        else:
-            raise NotImplementedError("`p1` can only be a RotProj if `p2` is "
-                                      "latlong!")
-
-    if isinstance(p2, RotProj):
-        shapes = reproject(shapes, p1, 'latlong')
-        reproject_points = p2
-    else:
-        reproject_points = partial(pyproj.transform, as_projection(p1),
-                                   as_projection(p2))
+    transformer = Transformer.from_crs(crs1, crs2)
 
     def _reproject_shape(shape):
-        return transform(reproject_points, shape)
+        return transform(transformer.transform, shape)
 
     if isinstance(shapes, pd.Series):
         return shapes.map(_reproject_shape)
@@ -143,11 +99,8 @@ def reproject(shapes, p1, p2):
 reproject.__doc__ = reproject_shapes.__doc__
 
 
-def compute_indicatormatrix(
-        orig,
-        dest,
-        orig_proj='latlong',
-        dest_proj='latlong'):
+
+def compute_indicatormatrix(orig, dest, orig_crs=4326, dest_crs=4326):
     """
     Compute the indicatormatrix
 
@@ -172,7 +125,7 @@ def compute_indicatormatrix(
     """
     orig = orig.geometry if isinstance(orig, gpd.GeoDataFrame) else orig
     dest = dest.geometry if isinstance(dest, gpd.GeoDataFrame) else dest
-    dest = reproject_shapes(dest, dest_proj, orig_proj)
+    dest = reproject_shapes(dest, dest_crs, orig_crs)
     indicator = sp.sparse.lil_matrix((len(dest), len(orig)), dtype=np.float)
     tree = STRtree(orig)
     idx = dict((id(o), i) for i, o in enumerate(orig))

environment.yaml

@@ -8,6 +8,7 @@ name: atlite
 channels:
   - defaults
   - conda-forge
+  - defaults
 dependencies:
   - python>=3.6
   - numpy
@@ -24,9 +25,10 @@ dependencies:
   - rtree
   - libspatialindex
   - cdsapi
+  - pyproj>=2.0
   - rasterio>=1.0
-  - geopandas
   - shapely
+  - geopandas
 
   # Recommended for pandas and xarray
   - bottleneck

setup.py

@@ -37,6 +37,7 @@
                       'rasterio>=1.0',
                       'shapely',
                       'progressbar2',
+                      'pyproj>=2',
                       'geopandas',
                       'cdsapi'],
     extras_require = {