From c36d94a1f2ebcbcb462a1635af25d631482254eb Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 10 Apr 2019 15:27:49 -0400 Subject: [PATCH 1/8] Add return index option to hand --- pysheds/grid.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/pysheds/grid.py b/pysheds/grid.py index 7adc729..d016e6f 100644 --- a/pysheds/grid.py +++ b/pysheds/grid.py @@ -1688,7 +1688,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 +1799,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 +1831,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: @@ -2592,6 +2595,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.') From b4494169468e820b24b801aa4ac0ebdc1cdc41be Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 14:18:06 -0400 Subject: [PATCH 2/8] Fix indent --- pysheds/rfsm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysheds/rfsm.py b/pysheds/rfsm.py index d3f13af..3a28284 100644 --- a/pysheds/rfsm.py +++ b/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: From 23de62c897286619363dc52dfd2ded4ba181067a Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 15:52:27 -0400 Subject: [PATCH 3/8] Add support for pyproj >2.2 --- pysheds/grid.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/pysheds/grid.py b/pysheds/grid.py index d016e6f..eb642e1 100644 --- a/pysheds/grid.py +++ b/pysheds/grid.py @@ -30,6 +30,10 @@ 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' + from pysheds.view import Raster from pysheds.view import BaseViewFinder, RegularViewFinder, IrregularViewFinder from pysheds.view import RegularGridViewer, IrregularGridViewer @@ -457,7 +461,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, @@ -1860,11 +1864,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(), @@ -1921,11 +1925,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(), @@ -2224,8 +2228,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 From 647307786e5514686449a6edaa348d50e623a24f Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 16:26:20 -0400 Subject: [PATCH 4/8] Fix issue with pyproj dependency update --- pysheds/grid.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/pysheds/grid.py b/pysheds/grid.py index eb642e1..3c2178e 100644 --- a/pysheds/grid.py +++ b/pysheds/grid.py @@ -428,8 +428,10 @@ 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] + # new_x, new_y = pyproj.transform(target_view.crs, as_crs, + # target_coords[:,1], target_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, @@ -444,8 +446,10 @@ 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] + # new_x, new_y = pyproj.transform(data_view.crs, target_view.crs, + # data_coords[:,1], data_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, From c22993338a2751699a96a98d75af3e595ef6a5e0 Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 17:56:28 -0400 Subject: [PATCH 5/8] Add fix to pyproj dependency in read_raster --- pysheds/grid.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pysheds/grid.py b/pysheds/grid.py index 3c2178e..35e1aa2 100644 --- a/pysheds/grid.py +++ b/pysheds/grid.py @@ -281,8 +281,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) From 8a745bdaa377067757ba88a56ed680f4448e0d06 Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 18:05:17 -0400 Subject: [PATCH 6/8] Fix pyproj deprecation warnings --- pysheds/grid.py | 7 ++++--- pysheds/view.py | 10 +++++++--- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/pysheds/grid.py b/pysheds/grid.py index 35e1aa2..3f692aa 100644 --- a/pysheds/grid.py +++ b/pysheds/grid.py @@ -33,6 +33,7 @@ _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 @@ -95,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 @@ -112,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 @@ -194,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 diff --git a/pysheds/view.py b/pysheds/view.py index 3d4d3ab..eff67fb 100644 --- a/pysheds/view.py +++ b/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 From bce9107cf4a4bfc3f5c4c88213468cf02a693af8 Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 18:24:18 -0400 Subject: [PATCH 7/8] Increase coverage with hand test --- tests/test_grid.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/test_grid.py b/tests/test_grid.py index 275dee7..d35d480 100644 --- a/tests/test_grid.py +++ b/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') From 5faa2ff33a59080232e92a27845fe6da209f0064 Mon Sep 17 00:00:00 2001 From: Matt Bartos Date: Wed, 1 Apr 2020 18:41:47 -0400 Subject: [PATCH 8/8] Remove commented blocks --- pysheds/grid.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pysheds/grid.py b/pysheds/grid.py index 3f692aa..4915171 100644 --- a/pysheds/grid.py +++ b/pysheds/grid.py @@ -436,8 +436,6 @@ def view(self, data, data_view=None, target_view=None, apply_mask=True, target_coords = target_view.coords new_coords = self._convert_grid_indices_crs(target_coords, target_view.crs, as_crs) new_x, new_y = new_coords[:,1], new_coords[:,0] - # new_x, new_y = pyproj.transform(target_view.crs, as_crs, - # target_coords[:,1], target_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, @@ -454,8 +452,6 @@ def view(self, data, data_view=None, target_view=None, apply_mask=True, # TODO: x and y order might be different 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] - # new_x, new_y = pyproj.transform(data_view.crs, target_view.crs, - # data_coords[:,1], data_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,