Merge pull request #108 from mdbartos/rfsm

Fix pyproj dependency

pysheds/grid.py

@@ -30,6 +30,11 @@
 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'
+
 from pysheds.view import Raster
 from pysheds.view import BaseViewFinder, RegularViewFinder, IrregularViewFinder
 from pysheds.view import RegularGridViewer, IrregularGridViewer
@@ -91,7 +96,7 @@ class Grid(object):
     """
 
     def __init__(self, affine=Affine(0,0,0,0,0,0), shape=(1,1), nodata=0,
-                 crs=pyproj.Proj('+init=epsg:4326'),
+                 crs=pyproj.Proj(_pyproj_init),
                  mask=None):
         self.affine = affine
         self.shape = shape
@@ -108,7 +113,7 @@ def defaults(self):
             'affine' : Affine(0,0,0,0,0,0),
             'shape' : (1,1),
             'nodata' : 0,
-            'crs' : pyproj.Proj('+init=epsg:4326'),
+            'crs' : pyproj.Proj(_pyproj_init),
         }
         return props
 
@@ -190,7 +195,7 @@ def add_gridded_data(self, data, data_name, affine=None, shape=None, crs=None,
         self.grids.append(data_name)
         setattr(self, data_name, data)
 
-    def read_ascii(self, data, data_name, skiprows=6, crs=pyproj.Proj('+init=epsg:4326'),
+    def read_ascii(self, data, data_name, skiprows=6, crs=pyproj.Proj(_pyproj_init),
                    xll='lower', yll='lower', metadata={}, **kwargs):
         """
         Reads data from an ascii file into a named attribute of Grid
@@ -277,8 +282,13 @@ def read_raster(self, data, data_name, band=1, window=None, window_crs=None,
                 if window_crs is not None:
                     if window_crs.srs != crs.srs:
                         xmin, ymin, xmax, ymax = window
-                        extent = pyproj.transform(window_crs, crs, (xmin, xmax),
-                                                  (ymin, ymax))
+                        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)
                         window = (extent[0][0], extent[1][0], extent[0][1], extent[1][1])
                 # If window crs not specified, assume it's in raster crs
                 ix_window = f.window(*window)
@@ -424,8 +434,8 @@ def view(self, data, data_view=None, target_view=None, apply_mask=True,
         if as_crs is not None:
             assert(isinstance(as_crs, pyproj.Proj))
             target_coords = target_view.coords
-            new_x, new_y = pyproj.transform(target_view.crs, as_crs,
-                                            target_coords[:,1], target_coords[:,0])
+            new_coords = self._convert_grid_indices_crs(target_coords, target_view.crs, as_crs)
+            new_x, new_y = new_coords[:,1], new_coords[:,0]
             # TODO: In general, crs conversion will yield irregular grid (though not necessarily)
             target_view = IrregularViewFinder(coords=np.column_stack([new_y, new_x]),
                                             shape=target_view.shape, crs=as_crs,
@@ -440,8 +450,8 @@ def view(self, data, data_view=None, target_view=None, apply_mask=True,
         if not same_crs:
             data_coords = data_view.coords
             # TODO: x and y order might be different
-            new_x, new_y = pyproj.transform(data_view.crs, target_view.crs,
-                                            data_coords[:,1], data_coords[:,0])
+            new_coords = self._convert_grid_indices_crs(data_coords, data_view.crs, target_view.crs)
+            new_x, new_y = new_coords[:,1], new_coords[:,0]
             # TODO: In general, crs conversion will yield irregular grid (though not necessarily)
             data_view = IrregularViewFinder(coords=np.column_stack([new_y, new_x]),
                                             shape=data_view.shape, crs=target_view.crs,
@@ -457,7 +467,7 @@ def view(self, data, data_view=None, target_view=None, apply_mask=True,
             # If spline interpolation is needed, use RectBivariate
             elif interpolation == 'spline':
                 # If latitude/longitude, use RectSphereBivariate
-                if target_view.crs.is_latlong():
+                if getattr(_pyproj_crs(target_view.crs), _pyproj_crs_is_geographic):
                     array_view = RegularGridViewer._view_rectspherebivariate(data, data_view,
                                                                              target_view,
                                                                              x_tolerance=tolerance,
@@ -1688,7 +1698,8 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None, dirmap=None, nodata_in=N
 
     def compute_hand(self, fdir, dem, drainage_mask, out_name='hand', dirmap=None,
                      nodata_in_fdir=None, nodata_in_dem=None, nodata_out=np.nan, routing='d8',
-                     inplace=True, apply_mask=False, ignore_metadata=False, **kwargs):
+                     inplace=True, apply_mask=False, ignore_metadata=False, return_index=False,
+                     **kwargs):
         """
         Computes the height above nearest drainage (HAND), based on a flow direction grid,
         a digital elevation grid, and a grid containing the locations of drainage channels.
@@ -1798,7 +1809,8 @@ def compute_hand(self, fdir, dem, drainage_mask, out_name='hand', dirmap=None,
                     hand[child] = hand[parent]
                     source = np.unique(child)
                 hand = hand.reshape(dem.shape)
-                hand = np.where(hand != -1, dem - dem.flat[hand], nodata_out)
+                if not return_index:
+                    hand = np.where(hand != -1, dem - dem.flat[hand], nodata_out)
             except:
                 raise
             finally:
@@ -1829,7 +1841,8 @@ def compute_hand(self, fdir, dem, drainage_mask, out_name='hand', dirmap=None,
                     hand[child] = hand[parent]
                     source = child
                 hand = hand.reshape(dem.shape)
-                hand = np.where(hand != -1, dem - dem.flat[hand], nodata_out)
+                if not return_index:
+                    hand = np.where(hand != -1, dem - dem.flat[hand], nodata_out)
             except:
                 raise
             finally:
@@ -1857,11 +1870,11 @@ def cell_area(self, out_name='area', nodata_out=0, inplace=True, as_crs=None):
                  CRS at which to compute the area of each cell.
         """
         if as_crs is None:
-            if self.crs.is_latlong():
+            if getattr(_pyproj_crs(self.crs), _pyproj_crs_is_geographic):
                 warnings.warn(('CRS is geographic. Area will not have meaningful '
                             'units.'))
         else:
-            if as_crs.is_latlong():
+            if getattr(_pyproj_crs(as_crs), _pyproj_crs_is_geographic):
                 warnings.warn(('CRS is geographic. Area will not have meaningful '
                             'units.'))
         indices = np.vstack(np.dstack(np.meshgrid(*self.grid_indices(),
@@ -1918,11 +1931,11 @@ def cell_distances(self, data, out_name='cdist', dirmap=None, nodata_in=None, no
         if routing.lower() != 'd8':
             raise NotImplementedError('Only implemented for D8 routing.')
         if as_crs is None:
-            if self.crs.is_latlong():
+            if getattr(_pyproj_crs(self.crs), _pyproj_crs_is_geographic):
                 warnings.warn(('CRS is geographic. Area will not have meaningful '
                             'units.'))
         else:
-            if as_crs.is_latlong():
+            if getattr(_pyproj_crs(as_crs), _pyproj_crs_is_geographic):
                 warnings.warn(('CRS is geographic. Area will not have meaningful '
                             'units.'))
         indices = np.vstack(np.dstack(np.meshgrid(*self.grid_indices(),
@@ -2221,8 +2234,13 @@ def _dy_dx(self):
     #     return new_bbox
 
     def _convert_grid_indices_crs(self, grid_indices, old_crs, new_crs):
-        x2, y2 = pyproj.transform(old_crs, new_crs, grid_indices[:,1],
-                                  grid_indices[:,0])
+        if _OLD_PYPROJ:
+            x2, y2 = pyproj.transform(old_crs, new_crs, grid_indices[:,1],
+                                    grid_indices[:,0])
+        else:
+            x2, y2 = pyproj.transform(old_crs, new_crs, grid_indices[:,1],
+                                      grid_indices[:,0], errcheck=True,
+                                      always_xy=True)
         yx2 = np.column_stack([y2, x2])
         return yx2
 
@@ -2592,6 +2610,8 @@ def extract_profiles(self, fdir, mask, dirmap=None, nodata_in=None, routing='d8'
         -------
         profiles : np.ndarray
                    Array of channel profiles
+        connections : dict
+                      Dictionary containing connections between channel profiles
         """
         if routing.lower() != 'd8':
             raise NotImplementedError('Only implemented for D8 routing.')

pysheds/rfsm.py

@@ -485,7 +485,7 @@ def volume_to_level(self, node, waterlevel):
                 minelev = np.nanmin(self.dem)
             target_vol = node.current_vol
             elev = optimize.bisect(self.compute_vol, minelev, maxelev,
-                                args=(node, target_vol))
+                                   args=(node, target_vol))
             if node.name:
                 mask = self.ws[node.level] == node.name
             else:

pysheds/view.py

@@ -3,6 +3,10 @@
 from scipy import interpolate
 import pyproj
 from affine import Affine
+from distutils.version import LooseVersion
+
+_OLD_PYPROJ = LooseVersion(pyproj.__version__) < LooseVersion('2.2')
+_pyproj_init = '+init=epsg:4326' if _OLD_PYPROJ else 'epsg:4326'
 
 class Raster(np.ndarray):
     def __new__(cls, input_array, viewfinder, metadata=None):
@@ -74,7 +78,7 @@ def dy_dx(self):
 
 class BaseViewFinder():
     def __init__(self, shape=None, mask=None, nodata=None,
-                 crs=pyproj.Proj('+init=epsg:4326'), y_coord_ix=0, x_coord_ix=1):
+                 crs=pyproj.Proj(_pyproj_init), y_coord_ix=0, x_coord_ix=1):
         if shape is not None:
             self.shape = shape
         else:
@@ -121,7 +125,7 @@ def size(self):
 
 class RegularViewFinder(BaseViewFinder):
     def __init__(self, affine, shape, mask=None, nodata=None,
-                 crs=pyproj.Proj('+init=epsg:4326'),
+                 crs=pyproj.Proj(_pyproj_init),
                  y_coord_ix=0, x_coord_ix=1):
         if affine is not None:
             self.affine = affine
@@ -219,7 +223,7 @@ def move_window(self, dxmin, dymin, dxmax, dymax):
 
 class IrregularViewFinder(BaseViewFinder):
     def __init__(self, coords, shape=None, mask=None, nodata=None,
-                 crs=pyproj.Proj('+init=epsg:4326'),
+                 crs=pyproj.Proj(_pyproj_init),
                  y_coord_ix=0, x_coord_ix=1):
         if coords is not None:
             self.coords = coords

tests/test_grid.py

@@ -18,7 +18,7 @@
 
 # Initialize grid
 grid = Grid()
-crs = pyproj.Proj('+init=epsg:4326', preserve_units=True)
+crs = pyproj.Proj('epsg:4326', preserve_units=True)
 grid.read_ascii(dir_path, 'dir', dtype=np.uint8, crs=crs)
 grid.read_raster(dem_path, 'dem')
 grid.read_raster(roi_path, 'roi')
@@ -37,8 +37,8 @@
 cells_in_catch = 11422
 catch_shape = (159, 169)
 max_distance = 209
-new_crs = pyproj.Proj('+init=epsg:3083')
-old_crs = pyproj.Proj('+init=epsg:4326', preserve_units=True)
+new_crs = pyproj.Proj('epsg:3083')
+old_crs = pyproj.Proj('epsg:4326', preserve_units=True)
 x, y = -97.29416666666677, 32.73749999999989
 
 
@@ -138,6 +138,9 @@ def test_accumulation():
     pos = np.where(grid.dinf_acc==grid.dinf_acc.max())
     assert(np.round(grid.dinf_acc[pos] / grid.dinf_acc_eff[pos]) == 4.)
 
+def test_hand():
+    grid.compute_hand('dir', 'dem', grid.acc > 100)
+
 def test_flow_distance():
     grid.clip_to('catch')
     grid.flow_distance(x, y, data='catch', dirmap=dirmap, out_name='dist', xytype='label')