From dc3175c1f1b0bb370be697dcc46b7ad1a49b4540 Mon Sep 17 00:00:00 2001 From: Aaron Date: Wed, 25 Mar 2020 10:18:03 -0700 Subject: [PATCH 1/5] Update dem.py --- dem.py | 344 +++++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 285 insertions(+), 59 deletions(-) diff --git a/dem.py b/dem.py index 2f48b6d..92bf7a5 100644 --- a/dem.py +++ b/dem.py @@ -1480,7 +1480,10 @@ def get_indexes_of_upstream_cells_for_location(self, x, y): def search_down_flow_direction(self, start): l = list() - ((row,col),) = self._xy_to_rowscols((start,)) + #comment out below for original + (row,col) = start + #uncomment below for original + #((row,col),) = self._xy_to_rowscols((start,)) while not (row == 0 or row == self._georef_info.ny-1 or col == 0 or col == self._georef_info.nx - 1 or (row,col) in l): l.append((row, col)) @@ -1494,7 +1497,10 @@ def search_down_flow_direction(self, start): def search_down_flow_direction_from_xy_location(self, start): rcs = self.search_down_flow_direction(start) - return self._rowscols_to_xy(rcs) + #comment out below for original + return rcs + #uncomment below for original + #return self._rowscols_to_xy(rcs) def convert_rivertools_directions_to_arc(self): # Function to convert river tools flow directions to arcGisFlowDirections @@ -2397,62 +2403,38 @@ def __calcD8Area(self, *args, **kwargs): class GeographicMainstemValleyArea(GeographicGridMixin, MainstemValleyArea): pass +class AlongFlowSmoothing(object): - -class KsFromChiWithSmoothing(BaseSpatialGrid): - required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',),'_create'), - (('ai_ascii_filename','EPSGprojectionCode'),'_read_ai'), - (('gdal_filename',), '_read_gdal'), - (('elevation', 'area', 'flow_direction', 'theta', 'vertical_interval'), '_create_from_elevation_area_flow_direction'), - (('elevation', 'area', 'flow_direction', 'theta', 'horizontal_interval'), '_create_from_elevation_area_flow_direction'), - ) - - def _upstream_downstream_indexes(self, area, flow_direction): - upstream_i = (np.ones_like(area._griddata)*-1.0).astype(int) - upstream_j = (np.ones_like(area._griddata)*-1.0).astype(int) - downstream_i = (np.ones_like(area._griddata)*-1.0).astype(int) - downstream_j = (np.ones_like(area._griddata)*-1.0).astype(int) + def _find_points_along_path(self, de, **kwargs): + area = kwargs['area'] + flow_direction = kwargs['flow_direction'] + + import time + t1 = time.time() + + upstream_i = (np.ones_like(area._griddata) * -1.0).astype(int) + upstream_j = (np.ones_like(area._griddata) * -1.0).astype(int) + downstream_i = (np.ones_like(area._griddata) * -1.0).astype(int) + downstream_j = (np.ones_like(area._griddata) * -1.0).astype(int) visited = np.zeros_like(area._griddata).astype(bool) shape = upstream_i.shape - indexes = area.sort(reverse = False) + indexes = area.sort(reverse=False) (i_s, j_s) = np.unravel_index(indexes, shape) + def set_usdsindexes(ij): i, j = ij ds_i, ds_j, good = flow_direction.get_flow_to_cell(i, j) - if good and not visited[ds_i,ds_j]: - visited[i,j] = True - (downstream_i[i, j], downstream_j[i,j]) = (ds_i, ds_j) + if good and not visited[ds_i, ds_j]: + visited[i, j] = True + (downstream_i[i, j], downstream_j[i, j]) = (ds_i, ds_j) (upstream_i[ds_i, ds_j], upstream_j[ds_i, ds_j]) = (i, j) - map(set_usdsindexes, zip(i_s, j_s)) - - return upstream_i, upstream_j, downstream_i, downstream_j - - def _create_from_elevation_area_flow_direction(self, *args, **kwargs): - - elevation = kwargs['elevation'] - area = kwargs['area'] - flow_direction = kwargs['flow_direction'] - theta = kwargs['theta'] - vertical_interval = kwargs.get('vertical_interval', None) - de = area._mean_pixel_dimension() - - self._copy_info_from_grid(elevation) - self._griddata = np.zeros_like(elevation._griddata) - self._n = np.zeros_like(self._griddata).astype(int) - self._n_regression = np.zeros_like(self._griddata).astype(int) - self._mse = np.zeros_like(self._griddata) - self._ss = np.zeros_like(self._griddata) - self._r2 = np.zeros_like(self._griddata) - self._pval = np.zeros_like(self._griddata) - self._griddata[:] = np.nan - - import time - t1 = time.time() - upstream_i, upstream_j, downstream_i, downstream_j = self._upstream_downstream_indexes(area, flow_direction) + + list(map(set_usdsindexes, zip(i_s, j_s))) + t2 = time.time() - - print('completed flow graph in: ' + str(t2-t1) + " s") - + print('completed flow graph in: ' + str(t2 - t1) + " s") + + vertical_interval = kwargs.get('vertical_interval', None) if vertical_interval is not None: def find_points_along_path(this_i, this_j): ret = list() @@ -2471,9 +2453,10 @@ def find_points_along_path(this_i, this_j): ret = [(ds_i, ds_j)] + ret + [(ups_i, ups_j)] delta_e = elevation._griddata[ups_i, ups_j] - elevation._griddata[ds_i, ds_j] return ret + return find_points_along_path else: horizontal_interval = kwargs['horizontal_interval'] - + def find_points_along_path(this_i, this_j): horizontal_distance = 0 ret = list() @@ -2483,21 +2466,51 @@ def find_points_along_path(this_i, this_j): if (ups_i < 0) or (ds_i < 0): return None ret = [(ds_i, ds_j)] + ret + [(ups_i, ups_j)] - horizontal_distance += ((1.0 if ((ds_i == this_i) or (ds_j == this_j)) else 1.414) + (1.0 if ((ups_i == this_i) or (ups_j == this_j)) else 1.414))*de[this_i, this_j] + horizontal_distance += ((1.0 if ((ds_i == this_i) or (ds_j == this_j)) else 1.414) + ( + 1.0 if ((ups_i == this_i) or (ups_j == this_j)) else 1.414)) * de[this_i, this_j] while (horizontal_distance < horizontal_interval) & (ups_i >= 0): (ups_i, ups_j) = (upstream_i[ups_i, ups_j], upstream_j[ups_i, ups_j]) (ds_i, ds_j) = (downstream_i[ds_i, ds_j], downstream_j[ds_i, ds_j]) if (ups_i < 0) or (ds_i < 0): return None - horizontal_distance += (1.0 if ((ds_i == ret[0][0]) or (ds_j == ret[0][1])) else 1.414)*de[ds_i, ds_j] + (1.0 if ((ups_i == ret[-1][0]) or (ups_j == ret[-1][1])) else 1.414)*de[ups_i, ups_j] + horizontal_distance += (1.0 if ((ds_i == ret[0][0]) or (ds_j == ret[0][1])) else 1.414) * de[ + ds_i, ds_j] + (1.0 if ((ups_i == ret[-1][0]) or (ups_j == ret[-1][1])) else 1.414) * de[ + ups_i, ups_j] ret = [(ds_i, ds_j)] + ret + [(ups_i, ups_j)] return ret + return find_points_along_path + +class KsFromChiWithSmoothing(BaseSpatialGrid, AlongFlowSmoothing): + required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',),'_create'), + (('ai_ascii_filename','EPSGprojectionCode'),'_read_ai'), + (('gdal_filename',), '_read_gdal'), + (('elevation', 'area', 'flow_direction', 'theta', 'vertical_interval'), '_create_from_elevation_area_flow_direction'), + (('elevation', 'area', 'flow_direction', 'theta', 'horizontal_interval'), '_create_from_elevation_area_flow_direction'), + ) + + def _create_from_elevation_area_flow_direction(self, *args, **kwargs): - def calc_ks(i,j): + elevation = kwargs['elevation'] + area = kwargs['area'] + theta = kwargs['theta'] + de = area._mean_pixel_dimension() + + self._copy_info_from_grid(elevation) + self._griddata = np.zeros_like(elevation._griddata) + self._n = np.zeros_like(self._griddata).astype(int) + self._n_regression = np.zeros_like(self._griddata).astype(int) + self._mse = np.zeros_like(self._griddata) + self._ss = np.zeros_like(self._griddata) + self._r2 = np.zeros_like(self._griddata) + self._pval = np.zeros_like(self._griddata) + self._griddata[:] = np.nan + + find_points_along_path = self._find_points_along_path(de, **kwargs) + def calc_ks(i,j): points = find_points_along_path(i, j) if points is not None: - pts = zip(*(points)) + pts = list(zip(*(points))) points = np.array(pts).astype(int) adjustment = np.ones((len(points[0]))) i = np.where((points[0,1:-1] != points[0,2:]) & (points[1,1:-1] != points[1,2:])) @@ -2519,7 +2532,7 @@ def calc_ks(i,j): return np.nan, np.nan, np.nan, np.nan, [[],[]], np.nan, 0 i = np.where((area._griddata != 0) & ~np.isnan(area._griddata) & ~np.isnan(elevation._griddata)) - ij = zip(i[0],i[1]) + ij = list(zip(i[0],i[1])) totalnumber = len(ij) counter = 0.0 next_readout = 0.1 @@ -2582,10 +2595,222 @@ def get_band(gdal_dataset, band_number): gdal_file = None return return_object - + + +class ThetaFromChiWithSmoothing(BaseSpatialGrid, AlongFlowSmoothing): + required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',), '_create'), + (('ai_ascii_filename', 'EPSGprojectionCode'), '_read_ai'), + (('gdal_filename',), '_read_gdal'), + (('elevation', 'area', 'flow_direction', 'vertical_interval', 'min_area'), + '_create_from_elevation_area_flow_direction'), + (('elevation', 'area', 'flow_direction', 'horizontal_interval', 'min_area'), + '_create_from_elevation_area_flow_direction'), + ) + + def _create_from_elevation_area_flow_direction(self, *args, **kwargs): + + elevation = kwargs['elevation'] + area = kwargs['area'] + min_area = kwargs['min_area'] + de = area._mean_pixel_dimension() + + self._copy_info_from_grid(elevation) + self._griddata = np.zeros_like(elevation._griddata) + self._n = np.zeros_like(self._griddata).astype(int) + self._n_regression = np.zeros_like(self._griddata).astype(int) + self._mse = np.zeros_like(self._griddata) + self._ss = np.zeros_like(self._griddata) + self._r2 = np.zeros_like(self._griddata) + self._pval = np.zeros_like(self._griddata) + self._griddata[:] = np.nan + + find_points_along_path = self._find_points_along_path(de, **kwargs) + + def calc_theta(i, j): + + points = find_points_along_path(i, j) + if points is not None: + pts = list(zip(*(points))) + points = np.array(pts).astype(int) + adjustment = np.ones((len(points[0]))) + i = np.where((points[0, 1:-1] != points[0, 2:]) & (points[1, 1:-1] != points[1, 2:])) + adjustment[i[0] + 1] += 0.414 + area_profile = area._griddata[points[0], points[1]] + elevation_profile = elevation._griddata[points[0], points[1]] + de_profile = de[points[0], points[1]] + + def r2_for_theta(theta): + if theta > 10 or theta < -10: + return np.inf + chi_profile = np.zeros_like(elevation_profile) + chi_profile[1:] = np.cumsum( + 0.25 * (np.power(area_profile[1:], -theta) + np.power(area_profile[0:-1], -theta)) * ( + de_profile[1:] + de_profile[0:-1]) * adjustment[1:]) + X = np.array(chi_profile) + y = np.array(elevation_profile) - elevation_profile[0] + model = sm.OLS(y, X) + res = model.fit() + return res.ssr + from scipy.optimize import fmin + try: + (theta_bf, funval, iter, funcalls, warnflag) = fmin(r2_for_theta, np.array([0.5]), (), 1E-5, 1E-5, 100, 200, True, False, 0, None) + chi_profile = np.zeros_like(elevation_profile) + chi_profile[1:] = np.cumsum( + 0.25 * (np.power(area_profile[1:], -theta_bf) + np.power(area_profile[0:-1], -theta_bf)) * ( + de_profile[1:] + de_profile[0:-1]) * adjustment[1:]) + X = np.array(chi_profile) + y = np.array(elevation_profile) - elevation_profile[0] + model = sm.OLS(y, X) + res = model.fit() + SS = res.ssr + return theta_bf, SS / float(len(chi_profile)), SS, res.rsquared, points, res.pvalues[0], len(chi_profile) + except: + return np.nan, np.nan, np.nan, np.nan, [[], []], np.nan, 0 + else: + + return np.nan, np.nan, np.nan, np.nan, [[],[]], np.nan, 0 + + i = np.where((area._griddata != 0) & ~np.isnan(area._griddata) & ~np.isnan(elevation._griddata) & (area._griddata > min_area)) + ij = list(zip(i[0], i[1])) + totalnumber = len(ij) + counter = 0.0 + next_readout = 0.1 + sys.stdout.write('Percent completion...') + sys.stdout.flush() + for (i, j) in ij: + self._griddata[i, j], self._mse[i, j], self._ss[i, j], self._r2[i, j], pts, self._pval[i, j], \ + self._n_regression[i, j] = calc_theta(i, j) + self._n[pts[0], pts[1]] += 1 + counter += 1.0 / totalnumber + if counter > next_readout: + sys.stdout.write(str(int(next_readout * 100)) + "...") + sys.stdout.flush() + next_readout += 0.1 + + sys.stdout.write('Percent completion...') + sys.stdout.flush() + sys.stdout.write('100') + sys.stdout.flush() + + def save(self, filename): + + self._create_gdal_representation_from_array(self._georef_info, 'GTiff', + [self._griddata, self._n, self._mse, self._ss, self._r2, self._pval, + self._n_regression], self.dtype, filename, + ['COMPRESS=LZW', 'BIGTIFF=YES'], multiple_bands=True) + + @classmethod + def load(cls, filename): + + def get_band(gdal_dataset, band_number): + band = gdal_dataset.GetRasterBand(band_number) + nodata = band.GetNoDataValue() + grid = band.ReadAsArray().astype(cls.dtype) + if nodata is not None: + nodata_elements = np.where(grid == nodata) + from numpy import uint8 + if cls.dtype is not uint8: + grid[nodata_elements] = np.NAN + return grid + + return_object = cls() + gdal_dataset = gdal.Open(filename) + + geoTransform = gdal_dataset.GetGeoTransform() + nx = gdal_dataset.RasterXSize + ny = gdal_dataset.RasterYSize + + return_object._georef_info.geoTransform = geoTransform + return_object._georef_info.dx = return_object._georef_info.geoTransform[1] + return_object._georef_info.xllcenter = return_object._georef_info.geoTransform[ + 0] + return_object._georef_info.dx / 2.0 + return_object._georef_info.yllcenter = return_object._georef_info.geoTransform[3] - ( + return_object._georef_info.dx * (ny - 0.5)) + return_object._georef_info.nx = nx + return_object._georef_info.ny = ny + + return_object._griddata = get_band(gdal_dataset, 1) + return_object._n = get_band(gdal_dataset, 2) + return_object._mse = get_band(gdal_dataset, 3) + return_object._ss = get_band(gdal_dataset, 4) + return_object._r2 = get_band(gdal_dataset, 5) + return_object._pval = get_band(gdal_dataset, 6) + return_object._n_regression = get_band(gdal_dataset, 7) + + gdal_file = None + return return_object + + + + + +class ChannelSlopeWithSmoothing(BaseSpatialGrid, AlongFlowSmoothing): + required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',), '_create'), + (('ai_ascii_filename', 'EPSGprojectionCode'), '_read_ai'), + (('gdal_filename',), '_read_gdal'), + (('elevation', 'area', 'flow_direction', 'vertical_interval'), + '_create_from_elevation_area_flow_direction'), + (('elevation', 'area', 'flow_direction', 'horizontal_interval'), + '_create_from_elevation_flow_direction'), + ) + + + def _create_from_elevation_flow_direction(self, *args, **kwargs): + + elevation = kwargs['elevation'] + de = elevation._mean_pixel_dimension() + + self._copy_info_from_grid(elevation) + self._griddata = np.zeros_like(elevation._griddata) + self._griddata[:] = np.nan + + find_points_along_path = self._find_points_along_path(de, **kwargs) + + def calc_channel_slope(i, j): + + points = find_points_along_path(i, j) + if points is not None: + pts = list(zip(*(points))) + points = np.array(pts).astype(int) + adjustment = np.ones((len(points[0]))) + ind = np.where((points[0, 1:-1] != points[0, 2:]) & (points[1, 1:-1] != points[1, 2:])) + adjustment[ind[0] + 1] += 0.414 + elevation_profile = elevation._griddata[points[0], points[1]] + de_profile = de[points[0], points[1]] + dx = np.sum(de_profile*adjustment) + dy = elevation_profile[-1] - elevation_profile[0] + return dy / dx + + else: + + return np.nan + + i = np.where(~np.isnan(elevation._griddata)) + ij = list(zip(i[0], i[1])) + totalnumber = len(ij) + counter = 0.0 + next_readout = 0.1 + sys.stdout.write('Percent completion...') + sys.stdout.flush() + for (i, j) in ij: + self._griddata[i, j] = calc_channel_slope(i, j) + counter += 1.0 / totalnumber + if counter > next_readout: + sys.stdout.write(str(int(next_readout * 100)) + "...") + sys.stdout.flush() + next_readout += 0.1 + + sys.stdout.write('Percent completion...') + sys.stdout.flush() + sys.stdout.write('100') + sys.stdout.flush() + class GeographicKsFromChiWithSmoothing(GeographicGridMixin, KsFromChiWithSmoothing): pass +class GeographicThetaFromChiWithSmoothing(GeographicGridMixin, KsFromChiWithSmoothing): + pass + class MultiscaleCurvatureValleyWidth(BaseSpatialGrid): required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',),'_create'), (('ai_ascii_filename','EPSGprojectionCode'),'_read_ai'), @@ -2720,14 +2945,15 @@ def _elevation_fit_for_location(cls, x, y, elevation, de, fix_center = False): # Determine location in index space: ((i,j),) = elevation._xy_to_rowscols(((x,y),)) ((xa, ya),) = elevation._rowscols_to_xy(((i,j),)) - a,b,c,d,e,f = cls.Utilities._calc_coefficients_for_scale(elevation._griddata, de) + a,b,c,d,e,f = cls.Utilities._calc_coefficients_for_scale(elevation, de) + elevation_center = elevation[i,j] a = a[i,j] b = b[i,j] c = c[i,j] d = d[i,j] e = e[i,j] - f = f[i,j] - print('a = ' + str(a) + ', b = ' + str(b) + ', c = ' + str(c) + ', d = ' + str(d) + ', e = ' + str(e) + ', f = ' + str(f)) + f = f[i,j] + elevation_center + print('Window size = ' + str(de) + '/n' + 'a = ' + str(a) + ', b = ' + str(b) + ', c = ' + str(c) + ', d = ' + str(d) + ', e = ' + str(e) + ', f = ' + str(f)) (nx, ny, dx) = (elevation._georef_info.nx, elevation._georef_info.ny, elevation._georef_info.dx) (xllcenter, yllcenter) = (elevation._georef_info.xllcenter, elevation._georef_info.xllcenter) From 1248f0faa9810f8f586aa73e627954d308c5fbe6 Mon Sep 17 00:00:00 2001 From: Aaron Date: Thu, 29 Apr 2021 09:38:03 -0700 Subject: [PATCH 2/5] Create generic_preprocessing Jupyter notebook file allowing for the full suite of pre-processing on a DEM to get area, fd, flow length, etc. --- .DS_Store | Bin 0 -> 8196 bytes .../Generic_Preprocessing-checkpoint.ipynb | 91 ++++++++++++++++++ Generic_Preprocessing.ipynb | 91 ++++++++++++++++++ 3 files changed, 182 insertions(+) create mode 100644 .DS_Store create mode 100644 .ipynb_checkpoints/Generic_Preprocessing-checkpoint.ipynb create mode 100644 Generic_Preprocessing.ipynb diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..e7ff40a29cc7c883feedff5926a734bed9c66d4a GIT binary patch literal 8196 zcmeHMTWl0n7(U-pV5SW8v{pL60t=NwXrqhET?tEzB2+?S@1lo7+V+wUa`8*n~(}Y~~Yg!`WR1pQ)GiI?Pte4r$WWhDsF(b}EoPjt4 zaR%ZH#2L65GC+5>Xi9~9UrOUK&On@jTQURU`;ek8Xf~o#l7^oS8p0BQs3-x#M$c3S zgf^*&W+OT!X=qA5W%Yo_6j6x*X-@K#P$!y==#-?=oI#p1L@OhzP!O(8cBxREAuDM- z#uC$&LGCImgS9#@D$q$M%DC`WzCos%Fobt0>h<_5A9G zJ4d|n&cFzYeydT~FJ>J^Fp`U&yFAai{A}cUY-@ao*SC3YVA!rz2#syiB9*~j+cmxM zRT2h%zX0jzUv1wDfCAB`YadJ{s8dj~_w)?>7k)spG zCXc@)BSiF-U^O!$ez|hEj$dk%;U?ZFhkKzMu39}$o3HPBU{|xYr_auNZl`_J8sarR zMZBA)9R;gI(e>s%r#R~RL%cR`IQgPu1lA765z54?@$&`GaXLKT4s6fu?+dKMLD|t& zqP_v&D+Fz!vmwrlzV(nrzU-2V|9<)W-!uk_<74UZFhKp^6UooLE9>*wP;bhTw~MO8V;K*P*)>tNxVUBe&jF=8;{};0vyFLOyW3B;u$=L=Wzal0}R4 z6)TfTE&!*5)27u)je`z^8cv&6lQcmx;aE#2?^X%+d{#J3ZqPO<1aLlEO1EnFs&oW+ zm6UGR2*>RzpA$N2TdHZCsp|D0`=5aQo?T$Svr9-|EX z^bkY`upfi4u7GI?<2Z~HcpOjQNdoJ$1lCh{0Wab-Ud3y8ok068-Xq9\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfilled_dem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFilledElevation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;32mdel\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 539\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mError\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInputError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Input Error'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Inputs not satisfied'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 540\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 541\u001b[0;31m \u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'self.'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mevaluative_action\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'(*args, **kwargs)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 542\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 543\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__setitem__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_create_from_elevation\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2254\u001b[0m \u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'elevation'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2255\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_copy_info_from_grid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2256\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2257\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2258\u001b[0m \u001b[0;32mclass\u001b[0m \u001b[0mArea\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mBaseSpatialGrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_flood\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2171\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misEmpty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2172\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2173\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcol\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2175\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'binary_result'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'clip_to_fill'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 2105\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2106\u001b[0m \u001b[0;31m#Remove an item and its priority from the queue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2107\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheapq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheappop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__pq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2108\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__nItems\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2109\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], + "source": [ + "filled_dem = d.FilledElevation(elevation = dem)\n", + "del dem" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fd = d.FlowDirectionD8(flooded_dem = filled_dem)\n", + "del filled_dem\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "area = d.Area(flow_direction = fd)\n", + "area.save('drainages2016_area.tif')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Generic_Preprocessing.ipynb b/Generic_Preprocessing.ipynb new file mode 100644 index 0000000..7c3b5e1 --- /dev/null +++ b/Generic_Preprocessing.ipynb @@ -0,0 +1,91 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import TopoAnalysis.dem as d\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "dem = d.Elevation.load('/Volumes/Steelqat/DrainageWork/2016dem.tif')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfilled_dem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFilledElevation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;32mdel\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 539\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mError\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInputError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Input Error'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Inputs not satisfied'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 540\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 541\u001b[0;31m \u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'self.'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mevaluative_action\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'(*args, **kwargs)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 542\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 543\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__setitem__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_create_from_elevation\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2254\u001b[0m \u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'elevation'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2255\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_copy_info_from_grid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2256\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2257\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2258\u001b[0m \u001b[0;32mclass\u001b[0m \u001b[0mArea\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mBaseSpatialGrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_flood\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2171\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misEmpty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2172\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2173\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcol\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2175\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'binary_result'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'clip_to_fill'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 2105\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2106\u001b[0m \u001b[0;31m#Remove an item and its priority from the queue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2107\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheapq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheappop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__pq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2108\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__nItems\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2109\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], + "source": [ + "filled_dem = d.FilledElevation(elevation = dem)\n", + "del dem" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fd = d.FlowDirectionD8(flooded_dem = filled_dem)\n", + "del filled_dem\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "area = d.Area(flow_direction = fd)\n", + "area.save('drainages2016_area.tif')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 84e893ffbffce954f3c3bde969858c427ebda32a Mon Sep 17 00:00:00 2001 From: Aaron Date: Thu, 29 Apr 2021 09:56:27 -0700 Subject: [PATCH 3/5] Manually updated dem.py Overwrote old dem.py files --- MovingWindow.py | 10 +-- dem.py | 204 +++++++++++++++++++++++++++--------------------- 2 files changed, 119 insertions(+), 95 deletions(-) diff --git a/MovingWindow.py b/MovingWindow.py index eb61359..2042f9a 100644 --- a/MovingWindow.py +++ b/MovingWindow.py @@ -1,4 +1,4 @@ -import Error +import error import numpy as np class MovingWindow(object): @@ -7,11 +7,11 @@ class MovingWindow(object): def __init__(self, *args, **kwargs): if kwargs.get('window_dimension') is None: - raise Error.InputError('Window dimension','is a required parameter') + raise error.InputError('Window dimension', 'is a required parameter') self.window_dimension = kwargs.get('window_dimension') if self.__class__ is MovingWindow: - raise Error.Error('MovingWindow is an abstract base class') + raise error.Error('MovingWindow is an abstract base class') def _build_search_kernel(self, dx): return None, None @@ -58,7 +58,7 @@ class RectangularMovingWindow(MovingWindow): def __init__(self, *args, **kwargs): super(RectangularMovingWindow, self).__init__(*args, **kwargs) if self.__class__ is RectangularMovingWindow: - raise Error.Error('RectangularMovingWindow has no bound function') + raise error.Error('RectangularMovingWindow has no bound function') def _build_search_kernel(self, dx): # Function to build a search kernel that is either square or circular @@ -73,7 +73,7 @@ class CircularMovingWindow(MovingWindow): def __init__(self, *args, **kwargs): super(CircularMovingWindow, self).__init__(*args, **kwargs) if self.__class__ is CircularMovingWindow: - raise Error.Error('CircularMovingWindow has no bound function') + raise error.Error('CircularMovingWindow has no bound function') def _build_search_kernel(self, dx): # Function to build a search kernel that is either square or circular diff --git a/dem.py b/dem.py index 0ad1fd7..5d6f7b7 100644 --- a/dem.py +++ b/dem.py @@ -4,9 +4,7 @@ #If you somehow have this, and have never spoken to me, please reach out #- I'd be curious to hear what you are doing (maybe I can even help with something!) -from osgeo import gdal # Used to load gis in files -import osr # Used with gis files -from osgeo import ogr +from osgeo import gdal, ogr, osr # Used to load gis in files import os # Used to join file paths, iteract with os in other ways.. import glob # Used for finding files that I want to mosaic (by allowing wildcard searches of the filesystem) import heapq # Used for constructing priority queue, which is used for filling dems @@ -16,6 +14,7 @@ from matplotlib import pyplot as plt import sys from scipy import stats +from . import error as Error try: import statsmodels.api as sm @@ -36,7 +35,7 @@ def _get_projection_from_EPSG_projection_code(self, EPSGprojectionCode): def _get_gdal_type_for_numpy_type(self, numpy_type): from numpy import float64, uint8, uint16, int16, uint32, int32, float32, complex64 - from gdal import GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32, GDT_Int32, GDT_Float32, GDT_Float64, GDT_CFloat64, GDT_Unknown + from osgeo.gdal import GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32, GDT_Int32, GDT_Float32, GDT_Float64, GDT_CFloat64, GDT_Unknown type_map = { uint8: GDT_Byte, uint16: GDT_UInt16, @@ -57,7 +56,7 @@ def _get_gdal_type_for_numpy_type(self, numpy_type): def _get_numpy_type_for_gdal_type(self, gdal_type): from numpy import float64, uint8, uint16, int16, uint32, int32, float32, complex64 - from gdal import GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32, GDT_Int32, GDT_Float32, GDT_Float64, GDT_CFloat64 + from osgeo.gdal import GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32, GDT_Int32, GDT_Float32, GDT_Float64, GDT_CFloat64 type_map = { GDT_Byte: uint8, GDT_UInt16: uint16, @@ -221,12 +220,12 @@ def _getRasterGeoTransformFromAsciiRaster(self, fileName): georef_data = dict() with open(fileName, "r") as ascii_file: - for _ in xrange(6): + for _ in range(6): line = ascii_file.readline() (key, value) = (line.split()[0], float(line.split()[-1])) georef_data[key.lower()] = value - required_values = ('ncols', 'nrows', 'cellsize','nodata_value') + required_values = ('ncols', 'nrows', 'cellsize') if len(set(required_values).difference(set(georef_data.keys()))) != 0: raise Error.InputError('A/I ASCII grid error','The following properties are missing: ' + str(set(required_values).difference(set(georef_data.keys())))) @@ -236,7 +235,7 @@ def _getRasterGeoTransformFromAsciiRaster(self, fileName): if georef_data.get('yllcorner') is None and georef_data.get('yllcenter') is None: raise Error.InputError('A/I ASCII grid error','Neither YLLCorner nor YLLCenter is present.') - + dx = georef_data.get('cellsize') nx = int(georef_data.get('ncols')) ny = int(georef_data.get('nrows')) @@ -666,7 +665,13 @@ def _area_per_pixel(self, *args, **kwargs): def _mean_pixel_dimension(self, *args, **kwargs): return self._georef_info.dx * np.ones_like(self._griddata, dtype = float64) - + + def get_XY_matricies(self): + xllc, yllc, nx, ny, dx = (self._georef_info.xllcenter, self._georef_info.yllcenter, self._georef_info.nx, self._georef_info.ny, self._georef_info.dx) + x = np.arange(xllc, xllc+(nx-1)*dx, dx) + y = np.arange(yllc, yllc+(ny-1)*dx, dx) + return np.meshgrid(x,y) + def resample(self, de): return_grid = self.__class__() return_grid._copy_info_from_grid(self, set_zeros=True) @@ -946,7 +951,7 @@ def principal_curvatures(self): def plot(self, **kwargs): interactive = kwargs.pop('interactive', True) - colorbar = kwargs.pop('colorbar', False) + colorbar = kwargs.pop('colorbar', True) extent = [self._georef_info.xllcenter, self._georef_info.xllcenter+(self._georef_info.nx-0.5)*self._georef_info.dx, self._georef_info.yllcenter, self._georef_info.yllcenter+(self._georef_info.ny-0.5)*self._georef_info.dx] plt.imshow(self._griddata, extent = extent, **kwargs) if interactive: @@ -954,7 +959,6 @@ def plot(self, **kwargs): plt.show(block=False) else: plt.ioff() - plt.show(block=True) if colorbar: plt.colorbar() @@ -1307,48 +1311,54 @@ def get_upstream_cell_indexes(self, i, j): return options - def __get_flow_from_cell(self, i, j): + def __get_flow_from_cell(self, i, j, max_recursion_depth=None, depth=0): i_source = [i] j_source = [j] - + + if max_recursion_depth is not None and depth > max_recursion_depth: + return i_source, j_source + + recursion_depth = depth + 1 + + if self[i, j+1] == 16: - [i_append, j_append] = self.__get_flow_from_cell(i,j+1) + [i_append, j_append] = self.__get_flow_from_cell(i,j+1,max_recursion_depth=max_recursion_depth,depth=recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i+1, j+1] == 32: - [i_append, j_append] = self.__get_flow_from_cell(i+1,j+1) + [i_append, j_append] = self.__get_flow_from_cell(i+1,j+1,max_recursion_depth=max_recursion_depth,depth=recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i+1, j] == 64: - [i_append, j_append] = self.__get_flow_from_cell(i+1,j) + [i_append, j_append] = self.__get_flow_from_cell(i+1,j,max_recursion_depth=max_recursion_depth,depth=recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i+1, j-1] == 128: - [i_append, j_append] = self.__get_flow_from_cell(i+1,j-1) + [i_append, j_append] = self.__get_flow_from_cell(i+1,j-1,max_recursion_depth=max_recursion_depth,depth=recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i, j-1] == 1: - [i_append, j_append] = self.__get_flow_from_cell(i,j-1) + [i_append, j_append] = self.__get_flow_from_cell(i,j-1,max_recursion_depth=max_recursion_depth,depth=recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i-1, j-1] == 2: - [i_append, j_append] = self.__get_flow_from_cell(i-1,j-1) + [i_append, j_append] = self.__get_flow_from_cell(i-1,j-1,max_recursion_depth=max_recursion_depth,depth=recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i-1, j] == 4: - [i_append, j_append] = self.__get_flow_from_cell(i-1,j) + [i_append, j_append] = self.__get_flow_from_cell(i-1,j,max_recursion_depth = max_recursion_depth,depth = recursion_depth) i_source = i_source + i_append j_source = j_source + j_append if self[i-1, j+1] == 8: - [i_append, j_append] = self.__get_flow_from_cell(i-1,j+1) + [i_append, j_append] = self.__get_flow_from_cell(i-1,j+1,max_recursion_depth = max_recursion_depth,depth = recursion_depth) i_source = i_source + i_append j_source = j_source + j_append @@ -1487,19 +1497,24 @@ def get_indexes_of_upstream_cells_for_location(self, x, y): ((i, j),) = self._xy_to_rowscols(v) return self.get_indexes_of_upstream_cells(i, j) - def search_down_flow_direction(self, start, search_length = np.inf): + def search_down_flow_direction(self, start, search_length = np.inf, len_out = False): l = list() - + len = list() ((row,col),) = self._xy_to_rowscols((start,)) length = 0.0 while not (row == 0 or row == self._georef_info.ny-1 or col == 0 or col == self._georef_info.nx - 1 or \ (row,col) in l) and (length < search_length): l.append((row, col)) + if len_out: + len.append(length) row,col,inBounds = self.get_flow_to_cell(row, col) length += self._georef_info.dx * (1.0 if (row == col) else 1.414) if not inBounds: break - return tuple(l) + if len_out: + return tuple(l), tuple(len) + else: + return tuple(l) def search_down_flow_direction_from_rowscols_location(self, start, return_rowscols = False, search_length = np.inf): @@ -2081,7 +2096,7 @@ class GeographicLaplacian(GeographicGridMixin, Laplacian): class PriorityQueueMixIn(object): - aggradation_slope = 0.000000001 + aggradation_slope = 1E-12 class priorityQueue: #Implements a priority queue using heapq. Python has a priority queue module built in, but it @@ -2183,9 +2198,7 @@ def _flood(self, *args, **kwargs): if not closed[neighborRows[i], neighborCols[i]]: #If this was a hole (lower than the cell downstream), fill it if self._griddata[neighborRows[i], neighborCols[i]] <= elevation: - - - + if kwargs.get('maximum_pit_depth'): fill_depth = elevation - self._griddata[neighborRows[i], neighborCols[i]] if fill_depth > kwargs.get('maximum_pit_depth'): @@ -2437,7 +2450,6 @@ def _find_points_along_path(self, de, **kwargs): upstream_j = (np.ones_like(area._griddata) * -1.0).astype(int) downstream_i = (np.ones_like(area._griddata) * -1.0).astype(int) downstream_j = (np.ones_like(area._griddata) * -1.0).astype(int) - shape = upstream_i.shape indexes = area.sort(reverse=False) (i_s, j_s) = np.unravel_index(indexes, shape) @@ -2774,11 +2786,22 @@ class ChannelSlopeWithSmoothing(BaseSpatialGrid, AlongFlowSmoothing): '_create_from_elevation_flow_direction'), ) + scale_factor = 1.0 + + def points_along_path(self, points, i, j): + return points def calc_channel_slope(self, i, j, elevation, de, find_points_along_path): points = find_points_along_path(i, j) + + if points is None or len(points) == 0: + return np.nan + + points = self.points_along_path(points, i, j) + if points is not None: + pts = list(zip(*(points))) points = np.array(pts).astype(int) adjustment = np.ones((len(points[0]))) @@ -2797,6 +2820,11 @@ def calc_channel_slope(self, i, j, elevation, de, find_points_along_path): def _create_from_elevation_flow_direction(self, *args, **kwargs): elevation = kwargs['elevation'] + if kwargs.get('horizontal_interval') is not None: + kwargs['horizontal_interval'] = kwargs['horizontal_interval']*self.scale_factor + if kwargs.get('vertical_interval') is not None: + kwargs['vertical_interval'] = kwargs['vertical_interval']*self.scale_factor + de = elevation._mean_pixel_dimension() self._copy_info_from_grid(elevation) @@ -2805,7 +2833,6 @@ def _create_from_elevation_flow_direction(self, *args, **kwargs): find_points_along_path = self._find_points_along_path(de, **kwargs) - i = np.where(~np.isnan(elevation._griddata)) ij = list(zip(i[0], i[1])) totalnumber = len(ij) @@ -2825,7 +2852,7 @@ def _create_from_elevation_flow_direction(self, *args, **kwargs): sys.stdout.write('100') sys.stdout.flush() -class ChannelDownSlopeWithSmoothing(BaseSpatialGrid, AlongFlowSmoothing): +class ChannelDownSlopeWithSmoothing(ChannelSlopeWithSmoothing): required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',), '_create'), (('ai_ascii_filename', 'EPSGprojectionCode'), '_read_ai'), (('gdal_filename',), '_read_gdal'), @@ -2834,59 +2861,26 @@ class ChannelDownSlopeWithSmoothing(BaseSpatialGrid, AlongFlowSmoothing): (('elevation', 'area', 'flow_direction', 'horizontal_interval'), '_create_from_elevation_flow_direction'), ) + scale_factor = 2.0 - def calc_channel_slope(self, i, j, elevation, de, find_points_along_path): - - points = find_points_along_path(i, j) - - if points is not None: - number_of_points = int(len(points) / 2.0)+1 - points = points[0:number_of_points] - pts = list(zip(*(points))) - points = np.array(pts).astype(int) - adjustment = np.ones((len(points[0]))) - ind = np.where((points[0, 1:-1] != points[0, 2:]) & (points[1, 1:-1] != points[1, 2:])) - adjustment[ind[0] + 1] += 0.414 - elevation_profile = elevation._griddata[points[0], points[1]] - de_profile = de[points[0], points[1]] - dx = np.sum(de_profile * adjustment) - dy = elevation_profile[-1] - elevation_profile[0] - return dy / dx - - else: - - return np.nan + def points_along_path(self, points, i, j): + position_of_center = list([ind[0] for ind in zip(range(len(points)),points) if (ind[1][0] == i and ind[1][1] == j)])[0] + return points[0:position_of_center] if len(points[0:position_of_center]) > 0 else None - def _create_from_elevation_flow_direction(self, *args, **kwargs): - - elevation = kwargs['elevation'] - de = elevation._mean_pixel_dimension() - - self._copy_info_from_grid(elevation) - self._griddata = np.zeros_like(elevation._griddata) - self._griddata[:] = np.nan - - find_points_along_path = self._find_points_along_path(de, **kwargs) - - i = np.where(~np.isnan(elevation._griddata)) - ij = list(zip(i[0], i[1])) - totalnumber = len(ij) - counter = 0.0 - next_readout = 0.1 - sys.stdout.write('Percent completion...') - sys.stdout.flush() - for (i, j) in ij: +class ChannelUpSlopeWithSmoothing(ChannelSlopeWithSmoothing): + required_inputs_and_actions = ((('nx', 'ny', 'projection', 'geo_transform',), '_create'), + (('ai_ascii_filename', 'EPSGprojectionCode'), '_read_ai'), + (('gdal_filename',), '_read_gdal'), + (('elevation', 'area', 'flow_direction', 'vertical_interval'), + '_create_from_elevation_area_flow_direction'), + (('elevation', 'area', 'flow_direction', 'horizontal_interval'), + '_create_from_elevation_flow_direction'), + ) + scale_factor = 2.0 - self._griddata[i, j] = self.calc_channel_slope(i, j, elevation, de, find_points_along_path) - counter += 1.0 / totalnumber - if counter > next_readout: - sys.stdout.write(str(int(next_readout * 100)) + "...") - sys.stdout.flush() - next_readout += 0.1 - sys.stdout.write('Percent completion...') - sys.stdout.flush() - sys.stdout.write('100') - sys.stdout.flush() + def points_along_path(self, points, i, j): + position_of_center = list([ind[0] for ind in zip(range(len(points)),points) if (ind[1][0] == i and ind[1][1] == j)])[0] + return points[position_of_center:] if len(points[position_of_center:]) > 0 else None class GeographicKsFromChiWithSmoothing(GeographicGridMixin, KsFromChiWithSmoothing): pass @@ -2935,10 +2929,10 @@ def _convolve(cls, X, Y, Z, K, fix_center=False): from numpy import fliplr as fliplr from numpy.fft import ifftshift - Xt = fliplr(X) - Yt = flipud(Y) + Xt = X + Yt = Y - FZ = fft2(Z._griddata) + FZ = fft2(fliplr(fliplr(Z._griddata))) FX1 = fft2(np.power(Xt, 2)) FX2 = fft2(np.power(Yt, 2)) FX3 = fft2(Xt * Yt) @@ -3025,24 +3019,54 @@ def _Cmin_for_scale(cls, Z, de, fix_center=False): @classmethod def _elevation_fit_for_location(cls, x, y, elevation, de, fix_center = False): + Z = Elevation() + Z._copy_info_from_grid(elevation, set_zeros=False) + + needs_reshaping_x = ((Z._georef_info.nx % 2) != 1) + needs_reshaping_y = ((Z._georef_info.ny % 2) != 1) + + if needs_reshaping_x: + nx = Z._georef_info.nx + 1 + xllcenter = Z._georef_info.xllcenter - Z._georef_info.dx + Z._georef_info.nx = nx + Z._georef_info.xllcenter = xllcenter + Z._griddata = np.concatenate(((np.reshape(Z._griddata[:, 0], (Z._georef_info.ny, 1)), Z._griddata)), axis=1) + if needs_reshaping_y: + ny = Z._georef_info.ny + 1 + yllcenter = Z._georef_info.yllcenter - Z._georef_info.dx + Z._georef_info.ny = ny + Z._georef_info.yllcenter = yllcenter + Z._griddata = np.concatenate(((np.reshape(Z._griddata[0, :], (1, Z._georef_info.nx)), Z._griddata)), axis=0) + # Determine location in index space: ((i,j),) = elevation._xy_to_rowscols(((x,y),)) ((xa, ya),) = elevation._rowscols_to_xy(((i,j),)) - a,b,c,d,e,f = cls.Utilities._calc_coefficients_for_scale(elevation, de) - elevation_center = elevation[i,j] + a,b,c,d,e,f = cls.Utilities._calc_coefficients_for_scale(Z, de) + elevation_center = Z[i,j] a = a[i,j] b = b[i,j] c = c[i,j] d = d[i,j] e = e[i,j] f = f[i,j] + elevation_center - print('Window size = ' + str(de) + '/n' + 'a = ' + str(a) + ', b = ' + str(b) + ', c = ' + str(c) + ', d = ' + str(d) + ', e = ' + str(e) + ', f = ' + str(f)) + print('Window size = ' + str(de) + '\n' + 'a = ' + str(a) + ', b = ' + str(b) + ', c = ' + str(c) + ', d = ' + str(d) + ', e = ' + str(e) + ', f = ' + str(f)) - (nx, ny, dx) = (elevation._georef_info.nx, elevation._georef_info.ny, elevation._georef_info.dx) - (xllcenter, yllcenter) = (elevation._georef_info.xllcenter, elevation._georef_info.xllcenter) + (nx, ny, dx) = (Z._georef_info.nx, Z._georef_info.ny, Z._georef_info.dx) + (xllcenter, yllcenter) = (Z._georef_info.xllcenter, Z._georef_info.yllcenter) [X,Y] = np.meshgrid(np.arange(xllcenter+dx/2, xllcenter + (nx-0.5)*dx, dx),np.arange(yllcenter+dx/2, yllcenter + (ny-0.5)*dx, dx)) - return a*np.power(X,2) + b*np.power(Y,2) + c*X*Y + d*X + e*Y + f + X -= xa + Y -= ya + + Z = Elevation() + Z._copy_info_from_grid(elevation, set_zeros=True) + Z._griddata = a*np.power(X,2) + b*np.power(Y,2) + c*X*Y + d*X + e*Y + f + + start_x_index = 1 if needs_reshaping_x else 0 + start_y_index = 1 if needs_reshaping_y else 0 + Z._griddata = Z._griddata[start_y_index:, start_x_index:] + + return Z def _create_from_inputs(self, *args, **kwargs): From 34398b93634246cb4968e6ce65398d673124f018 Mon Sep 17 00:00:00 2001 From: Aaron Date: Mon, 3 May 2021 09:16:09 -0700 Subject: [PATCH 4/5] Cleaned up branch Removed old files no longer needed --- .DS_Store | Bin 8196 -> 6148 bytes Generic_Preprocessing.ipynb | 91 ------------------------------------ 2 files changed, 91 deletions(-) delete mode 100644 Generic_Preprocessing.ipynb diff --git a/.DS_Store b/.DS_Store index e7ff40a29cc7c883feedff5926a734bed9c66d4a..9439717ef7222aed0a7a098e026472e68ef7e447 100644 GIT binary patch delta 399 zcmZp1XfcprU|?W$DortDU=RQ@Ie-{Mvv5r;6q~50D9Q?w2a6>$Br+s3WH6*M#BVHI z&a4iS;bw>j!UBd$pbT7=AwH?RxF9JfKMAPf#AG9Z1}^dHYBLiH9R*YK$$JH48D~ts zEg&C06So0n!9{sF`FZIO%NfKNf`HcMGUPLq0O??$b;V#-E|3lZ%M>uUF=R620NKTp z)dc!E3@y!c6bvm*fz~Fn7XX_{&cHCH4z_^&5gF}!R=zJg$;06+|pfKH7_?>w&zl}o6H#pd literal 8196 zcmeHMTWl0n7(U-pV5SW8v{pL60t=NwXrqhET?tEzB2+?S@1lo7+V+wUa`8*n~(}Y~~Yg!`WR1pQ)GiI?Pte4r$WWhDsF(b}EoPjt4 zaR%ZH#2L65GC+5>Xi9~9UrOUK&On@jTQURU`;ek8Xf~o#l7^oS8p0BQs3-x#M$c3S zgf^*&W+OT!X=qA5W%Yo_6j6x*X-@K#P$!y==#-?=oI#p1L@OhzP!O(8cBxREAuDM- z#uC$&LGCImgS9#@D$q$M%DC`WzCos%Fobt0>h<_5A9G zJ4d|n&cFzYeydT~FJ>J^Fp`U&yFAai{A}cUY-@ao*SC3YVA!rz2#syiB9*~j+cmxM zRT2h%zX0jzUv1wDfCAB`YadJ{s8dj~_w)?>7k)spG zCXc@)BSiF-U^O!$ez|hEj$dk%;U?ZFhkKzMu39}$o3HPBU{|xYr_auNZl`_J8sarR zMZBA)9R;gI(e>s%r#R~RL%cR`IQgPu1lA765z54?@$&`GaXLKT4s6fu?+dKMLD|t& zqP_v&D+Fz!vmwrlzV(nrzU-2V|9<)W-!uk_<74UZFhKp^6UooLE9>*wP;bhTw~MO8V;K*P*)>tNxVUBe&jF=8;{};0vyFLOyW3B;u$=L=Wzal0}R4 z6)TfTE&!*5)27u)je`z^8cv&6lQcmx;aE#2?^X%+d{#J3ZqPO<1aLlEO1EnFs&oW+ zm6UGR2*>RzpA$N2TdHZCsp|D0`=5aQo?T$Svr9-|EX z^bkY`upfi4u7GI?<2Z~HcpOjQNdoJ$1lCh{0Wab-Ud3y8ok068-Xq9\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfilled_dem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFilledElevation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;32mdel\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 539\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mError\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInputError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Input Error'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Inputs not satisfied'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 540\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 541\u001b[0;31m \u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'self.'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mevaluative_action\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'(*args, **kwargs)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 542\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 543\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__setitem__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_create_from_elevation\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2254\u001b[0m \u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'elevation'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2255\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_copy_info_from_grid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2256\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2257\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2258\u001b[0m \u001b[0;32mclass\u001b[0m \u001b[0mArea\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mBaseSpatialGrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_flood\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2171\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misEmpty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2172\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2173\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcol\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2175\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'binary_result'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'clip_to_fill'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 2105\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2106\u001b[0m \u001b[0;31m#Remove an item and its priority from the queue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2107\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheapq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheappop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__pq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2108\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__nItems\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2109\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mKeyboardInterrupt\u001b[0m: " - ] - } - ], - "source": [ - "filled_dem = d.FilledElevation(elevation = dem)\n", - "del dem" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fd = d.FlowDirectionD8(flooded_dem = filled_dem)\n", - "del filled_dem\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "area = d.Area(flow_direction = fd)\n", - "area.save('drainages2016_area.tif')" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.8" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From 5faa4a80c5ee0bcb09150894cdfee24ab8398c96 Mon Sep 17 00:00:00 2001 From: Aaron Date: Mon, 3 May 2021 09:28:35 -0700 Subject: [PATCH 5/5] Update .gitignore --- .DS_Store | Bin 6148 -> 6148 bytes .gitignore | 2 + .../Generic_Preprocessing-checkpoint.ipynb | 91 ------------------ 3 files changed, 2 insertions(+), 91 deletions(-) delete mode 100644 .ipynb_checkpoints/Generic_Preprocessing-checkpoint.ipynb diff --git a/.DS_Store b/.DS_Store index 9439717ef7222aed0a7a098e026472e68ef7e447..b061ad73257259ea7915bd04dc12ec7be9d70a10 100644 GIT binary patch delta 28 kcmZoMXfc@J&&aYdU^g=(%Vr*y2&RoCiy1ewbNuB80Dhtf#Q*>R delta 80 zcmZoMXfc@J&&awlU^g=(>t-I72qt5GhEj$+h7yKMAkJY(Wk_LAU`S)g2eMOuyrSIv g6qlr&{3Hej29CqpFK6p09CqATSircMo#QV*04kmp#Q*>R diff --git a/.gitignore b/.gitignore index 614013b..788ab00 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ *.pyc .idea .idea/* +.ipynb_checkpoints +*/.ipynb_checkpoints/* \ No newline at end of file diff --git a/.ipynb_checkpoints/Generic_Preprocessing-checkpoint.ipynb b/.ipynb_checkpoints/Generic_Preprocessing-checkpoint.ipynb deleted file mode 100644 index 7c3b5e1..0000000 --- a/.ipynb_checkpoints/Generic_Preprocessing-checkpoint.ipynb +++ /dev/null @@ -1,91 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import TopoAnalysis.dem as d\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "dem = d.Elevation.load('/Volumes/Steelqat/DrainageWork/2016dem.tif')\n" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [ - { - "ename": "KeyboardInterrupt", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfilled_dem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFilledElevation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;32mdel\u001b[0m \u001b[0mdem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 539\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mError\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInputError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Input Error'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Inputs not satisfied'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 540\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 541\u001b[0;31m \u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'self.'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mevaluative_action\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'(*args, **kwargs)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 542\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 543\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__setitem__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_create_from_elevation\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2254\u001b[0m \u001b[0melevation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'elevation'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2255\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_copy_info_from_grid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0melevation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2256\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2257\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2258\u001b[0m \u001b[0;32mclass\u001b[0m \u001b[0mArea\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mBaseSpatialGrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36m_flood\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2171\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misEmpty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2172\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2173\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcol\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpriority_queue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2175\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'binary_result'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'clip_to_fill'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Documents/Github/TopoAnalysis/dem.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 2105\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2106\u001b[0m \u001b[0;31m#Remove an item and its priority from the queue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2107\u001b[0;31m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheapq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheappop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__pq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2108\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__nItems\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2109\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mpriority\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mKeyboardInterrupt\u001b[0m: " - ] - } - ], - "source": [ - "filled_dem = d.FilledElevation(elevation = dem)\n", - "del dem" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fd = d.FlowDirectionD8(flooded_dem = filled_dem)\n", - "del filled_dem\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "area = d.Area(flow_direction = fd)\n", - "area.save('drainages2016_area.tif')" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.8" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -}