diff --git a/new_scripts/hillslopes.py b/new_scripts/hillslopes.py new file mode 100644 index 0000000..81d680b --- /dev/null +++ b/new_scripts/hillslopes.py @@ -0,0 +1,327 @@ +import os +import numpy as np + +import scipy.stats +import collections +import fiona, rasterio, shapely +import rasterio.warp + +import workflow +import workflow.crs +import workflow.warp + +import landcover +import datetime + + +def get_filenames(huc, huc_directory, raster_extension='tif'): + """Set up package directory for one huc to return a dictionary of filenames""" + ppm_dir = os.path.join(huc_directory,'data_preprocessed-meshing') + + filenames = dict() + def register_shapefile(key, filename): + filenames[key] = os.path.join(ppm_dir, filename)+'.shp' + + def register_raster(key, filename): + filenames[key] = os.path.join(ppm_dir, filename)+'.'+raster_extension + + # the HUC shapefile, in native and projected CRS + register_shapefile('huc', f'huc_{huc}') + register_shapefile('huc_proj', f'huc_{huc}_proj') + + # the HUC DEM, land_cover, slope, and aspect rasters + register_raster('dem', f'huc_{huc}_dem') # in units of [m] above sea level + #register_raster('dem_filled', f'huc_{huc}_dem_filled') + #register_raster('d8', f'huc_{huc}_d8') + #register_raster('d8i', f'huc_{huc}_d8i') + #register_raster('weights', f'huc_{huc}_flowpath_weights') + register_raster('land_cover', f'huc_{huc}_landcover') + register_raster('slope', f'huc_{huc}_slope') # in units of [-], e.g. rise over run, NOT % slope + register_raster('aspect', f'huc_{huc}_aspect') # in units of degrees clockwise from North (0 = N, 90 = E, 180 = S) + + # raster and shapefiles of the stream network + register_raster('streams_raster', f'huc_{huc}_streams') + register_shapefile('streams', f'huc_{huc}_streams') # the shapefile extracted from the above raster + register_shapefile('streams_network', f'huc_{huc}_streams_network') # simplified to a network for MOSART + register_shapefile('streams_network_proj', f'huc_{huc}_streams_network_proj') # projected form of the above + + # delineated subcatchments within the HUC + register_shapefile('subcatchments', f'huc_{huc}_subcatchments') # delineated subcatchments + + filenames['flowpaths'] = os.path.join(ppm_dir, 'flowpaths', f'hs_{{}}_flowpaths.pkl') + register_raster('flowpath_length', f'huc_{huc}_flowpath_lengths') # flowpaths within the delineated subcatchments + register_raster('elev_above_streams', f'huc_{huc}_elev_above_streams') + + return filenames + + +def save_shapefile(filename, shps, crs, extra_properties=None): + if len(shps) == 0: + return + + schema = dict() + if type(shps[0]) is shapely.geometry.Polygon: + schema['geometry'] = 'Polygon' + elif type(shps[0]) is shapely.geometry.LineString: + schema['geometry'] = 'LineString' + else: + raise RuntimeError('Currently this function only writes Polygon or LineString types') + + # set up the properties' schema, used for open and save geodata by fiona + schema['properties'] = collections.OrderedDict() + def register_type(key,atype): + if atype is int: + schema['properties'][key] = 'int' + elif atype is str: + schema['properties'][key] = 'str' + elif atype is float: + schema['properties'][key] = 'float' + else: + pass + if extra_properties is None: + extra_properties = dict() + for key, val in extra_properties.items(): + register_type(key, type(val)) + + try: + shp_property = shps[0].properties + except AttributeError: + pass + else: + for key, val in shp_property.items(): + register_type(key, type(val)) + + + with fiona.open(filename, 'w', + driver='ESRI Shapefile', + schema=schema, + crs=workflow.crs.to_fiona(crs), + crs_wkt=workflow.crs.to_wkt(crs)) as c: + for shp in shps: + props = extra_properties.copy() + try: + props.update(shp.properties) + except AttributeError: + pass + + for key in list(props.keys()): + if key not in schema['properties']: + props.pop(key) + + c.write({'geometry': shapely.geometry.mapping(shp), + 'properties': props}) + + +def save_demfile(filename, dem_profile, dem_raster): + rio_profile = dict(dem_profile).copy() + rio_profile.pop('blockxsize') + rio_profile.pop('blockysize') + rio_profile.pop('tiled') + rio_profile['nodata'] = -9999.0 + + rio_dem = np.where(np.isnan(dem_raster), rio_profile['nodata'], dem_raster).astype(rio_profile['dtype']) + with rasterio.open(filename, 'w', **rio_profile) as dst: + dst.write(rio_dem,1) + + +# Average aspect across the domain +def meanAspect(aspect): + ''' + # tests... + assert(meanAspect(np.array([90,90,90])) == 90) + assert(meanAspect(np.array([0,0,0])) == 0) + assert(meanAspect(np.array([180,180,180])) == 180) + assert(meanAspect(np.array([270,270,270])) == 270) + assert(meanAspect(np.array([89,90,91])) == 90) + assert(meanAspect(np.array([179,180,181])) == 180) + assert(meanAspect(np.array([269,270,271])) == 270) + assert(meanAspect(np.array([359,0,1])) == 0) + ''' + + a = aspect[~np.isnan(aspect)] + a = np.where(a > 180, a - 360, a) + sina = np.sin(np.radians(a)) + cosa = np.cos(np.radians(a)) + avg_aspect_radians = np.arctan2(sum(sina), sum(cosa)) + if avg_aspect_radians < 0: + avg_aspect_radians = 2*np.pi + avg_aspect_radians + return np.degrees(avg_aspect_radians) + + + +# Load a subcatchment shape +def loadSubcatchmentShape(filename, subcatch_id): + subcatch_crs, subcatchments = workflow.get_shapes(filename) + matches = [sc for sc in subcatchments if sc.properties['hs_id'] == subcatch_id] + if len(matches) != 1: + raise RuntimeError("Found invalid subcatchments file or no subcatchment of this id.") + return subcatch_crs, matches[0] + + +# Load a subcatchment raster +def loadSubcatchmentRaster(filename, subcatch, subcatch_crs, nanit=True): + profile, raster = workflow.get_raster_on_shape(filename, subcatch, subcatch_crs, mask=True) + if nanit: + raster[raster==profile['nodata']] = np.nan + return profile, raster + + +# Determine hillslope properties for a given subcatchment. +def parameterizeSubcatchment(filenames, huc, subcatch_id, + target_crs=workflow.crs.default_alaska_crs(), + hillslope_keep_fraction=0.95, + hillslope_bin_dx=100, + plot_smoothed=False): + + # Find the given subcatchment shape + subcatch_crs, subcatch = loadSubcatchmentShape(filenames['subcatchments'], subcatch_id) + + # Create a dictionary icluding all hillslope info and parameters + hillslope = dict() + hillslope['huc'] = huc + hillslope['subcatchment_id'] = subcatch_id + hillslope['centroid'] = subcatch.centroid.coords[0] # lat/long, required to get meteorological data + hillslope['subcatchment'] = workflow.warp.shply(subcatch, subcatch_crs, target_crs) + hillslope['subcatchment_native_crs'] = subcatch_crs + hillslope['subcatchment_target_crs'] = target_crs + hillslope['total_area'] = hillslope['subcatchment'].area # m^2 + + + # Procedures to determine a single hillslope profile geometry for each subcatchment + # Most of the parameters are based on bins in length of flowpath to the stream network + + # 1. Load raster of flowpath lengths for the given subcatchment + fp_profile, fp_lengths = loadSubcatchmentRaster(filenames['flowpath_length'], subcatch, subcatch_crs) + subcatch_mask = ~np.isnan(fp_lengths) + hillslope['raster_profile'] = fp_profile + + # 2. Ensure raster of flowpath lengths >= 0, and sort flowpath lengths + assert(fp_lengths[subcatch_mask].min() >= -1.e-3) + fp_lengths[fp_lengths < 0] = 0. + fp_lengths_masked = fp_lengths[subcatch_mask] + fp_lengths_masked_sorted = np.sort(fp_lengths_masked) + + # 3. Set bin parameters based on flowpath length + hillslope['total_length'] = fp_lengths_masked_sorted[-1] * hillslope_keep_fraction + hillslope['num_bins'] = int(np.round(hillslope['total_length'] / hillslope_bin_dx)) + hillslope['bin_dx'] = hillslope['total_length'] / hillslope['num_bins'] + assert(hillslope['num_bins'] > 3) + + # 4. Bin by flowpath length and save to dictionary + pixel_bin_raster = np.nan * np.ones(fp_lengths.shape, 'd') + pixel_bin_counts = np.zeros((hillslope['num_bins'],),'i') + + i = 0 + interval = 1 + while i < hillslope['num_bins']: + bin_start = i * hillslope['bin_dx'] + bin_end = (i+interval) * hillslope['bin_dx'] + local_mask = (fp_lengths >= bin_start) & (fp_lengths < bin_end) + pixel_bin_raster[local_mask] = i + pixel_bin_counts[i] = np.count_nonzero(local_mask) + + if pixel_bin_counts[i] > 0: + i += 1 + if i-1+interval>hillslope['num_bins']: + break + else: + interval += 1 + continue + + hillslope['bin_counts'] = pixel_bin_counts[np.nonzero(pixel_bin_counts)] + hillslope['num_bins'] = len(hillslope['bin_counts']) + hillslope['bins'] = pixel_bin_raster + +# assert(hillslope['bin_counts'].min() > 0) +# hillslope['bin_counts'] = pixel_bin_counts # THE OLD + + + + # 5. Average elevation (height over stream) for each bin and save to dictionary + _, elevs = loadSubcatchmentRaster(filenames['elev_above_streams'], subcatch, subcatch_crs) + elev_bins = np.zeros((hillslope['num_bins'],),'d') + for i in range(hillslope['num_bins']): + elev_bins[i] = elevs[(hillslope['bins'] == i) & (~np.isnan(elevs))].mean() + + hillslope['elevation'] = elev_bins + + # 6. Average aspect across the entire subcatchment and save to dictionary + _, aspects = loadSubcatchmentRaster(filenames['aspect'], subcatch, subcatch_crs) + hillslope['aspect'] = meanAspect(aspects) + + # 7. Get land cover using mode for each bin and save to dictionary + lc_profile, lc = loadSubcatchmentRaster(filenames['land_cover'], subcatch, subcatch_crs, False) # don't nan-it + + assert(lc_profile['nodata'] == 255) # uint8, nan is -1 == 255 + # classify land cover into veg classes + hillslope['lc'] = lc + num_missing, lc = landcover.classifyVegetation(lc) + lc_bins = np.zeros((hillslope['num_bins'],),'i') + + for i in range(hillslope['num_bins']): + bin_lc = lc[hillslope['bins'] == i] + assert(len(bin_lc) > 0) + lc_bins[i] = scipy.stats.mode(bin_lc.ravel())[0][0] + + hillslope['land_cover'] = lc_bins + hillslope['land_cover_raster'] = landcover.veg2img(lc) + + + # 8. Smooth subcatchment for easier 3D simulation + def smoothSubcatchmentShape(subcatch, subcatch_crs, smoothing_factor=100, plot=plot_smoothed): + if type(subcatch) is shapely.geometry.MultiPolygon: + subcatch_simp = shapely.ops.cascaded_union(subcatch.buffer(100)) + if type(subcatch_simp) is shapely.geometry.MultiPolygon: + # see if this is one tiny island + total_area = sum(p.area for p in subcatch_simp) + assert(total_area > 0) + polygons = [p for p in subcatch_simp if p.area/total_area > 0.01] + if len(polygons) > 1: + # give up, this will throw + return shapely.geometry.MultiPolygon(polygons) + + subcatch_simp = polygons[0] + subcatch_simp = subcatch_simp.simplify(smoothing_factor) + else: + subcatch_simp = subcatch.simplify(smoothing_factor) + + # find a buffer value that matches the original area + def optimize_func(radius): + subcatch_buff = subcatch_simp.buffer(radius).simplify(smoothing_factor) + return abs(subcatch_buff.area - subcatch.area) + + res = scipy.optimize.minimize_scalar(optimize_func) + subcatch_new = subcatch_simp.buffer(res.x).simplify(smoothing_factor) + + if plot: + print("error: ", subcatch.area - subcatch_new.area) + fig, ax = workflow.plot.get_ax(crs=subcatch_crs) + workflow.plot.shply(subcatch_new, subcatch_crs, 'r', ax=ax) + workflow.plot.shply(subcatch, subcatch_crs, 'b', ax=ax) + + return subcatch_new + + subcatch_new = smoothSubcatchmentShape(hillslope['subcatchment'], + hillslope['subcatchment_target_crs'], + smoothing_factor=100, plot=plot_smoothed) + + hillslope['subcatchment_smooth'] = subcatch_new + + return hillslope + + + +# Download DayMet +def downloadDaymet(hillslope_pars, raw_directory, save_filename, + start=datetime.date(1980,1,1), + end=datetime.date(2020,12,31)): + import daymet_to_ats + start, end = daymet_to_ats.validate_start_end(start, end) + lon, lat = hillslope_pars['centroid'] + daymet = daymet_to_ats.download_daymet(raw_directory, lat, lon, start, end) + ats = daymet_to_ats.daymet_to_ats(daymet) + attrs = daymet_to_ats.daymet_attrs(lat, lon, start, end) + + daymet_to_ats.write_ats(ats, attrs,save_filename) + + \ No newline at end of file diff --git a/new_scripts/huc_process.ipynb b/new_scripts/huc_process.ipynb new file mode 100644 index 0000000..db1e5ab --- /dev/null +++ b/new_scripts/huc_process.ipynb @@ -0,0 +1,1063 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a45b921c", + "metadata": {}, + "source": [ + "**InteRFACE Data Package Workflow**" + ] + }, + { + "cell_type": "markdown", + "id": "d8af26e3", + "metadata": {}, + "source": [ + "**Objective:** \n", + "* Download all needed data inputs for a ATS + MOSART run on a Sag River HUC.\n", + "* Process raw huc data and generate mesh.\n", + "* Generate a data package in a purely automated way." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "b92fa29a", + "metadata": {}, + "outputs": [], + "source": [ + "import sys,os\n", + "import numpy as np\n", + "import datetime\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import scipy.optimize, scipy.signal, scipy.stats\n", + "import collections\n", + "import logging\n", + "import fiona, rasterio, shapely\n", + "import rasterio.warp\n", + "\n", + "import workflow\n", + "import workflow.crs\n", + "import workflow.warp\n", + "import workflow.source_list\n", + "import workflow.utils\n", + "import workflow.ui\n", + "import workflow.conf\n", + "import workflow.mesh\n", + "\n", + "import hillslopes\n", + "import landcover\n", + "import meshing\n", + "import plot" + ] + }, + { + "cell_type": "markdown", + "id": "3ea2203a", + "metadata": {}, + "source": [ + "# Processing preparation" + ] + }, + { + "cell_type": "markdown", + "id": "ea31aacc", + "metadata": {}, + "source": [ + "## Input for this worksheet" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "b9c888d5", + "metadata": {}, + "outputs": [], + "source": [ + "# the HUC to delineate\n", + "huc = '190604020404'\n", + "\n", + "# contributing area, in pixels? [m^2]? used to provide a lower bound on pixels\n", + "# that are included in the stream network \n", + "streams_contributing_area_cutoff = -1 \n", + "\n", + "# target, in pixels, of the subcatchment size\n", + "target_subcatchment_size = 20000\n", + "\n", + "# number of horizontal grid cells in the hillslope\n", + "hillslope_bin_dx = 100\n", + "mesh_dx = 20\n", + "toeslope_min_slope = 0.0\n", + "hillslope_min_slope = 0.1\n", + "\n", + "#Don't let individual areas get too small in width -- 10% mean as a min value?\n", + "min_area_ratio = 0.1\n", + "\n", + "# Refine mesh area no larger than refine_max_area\n", + "refine_max_area = 1000\n", + "\n", + "# what fraction of the total flowpath lengths do we want to include?\n", + "#\n", + "# Effectively, some small number of pixels are often a long way away from the stream\n", + "# (whether this is real or artifact). We don't need tiny cells up at the top of the\n", + "# hillslope. Maybe keep 95% of the flowpath length?\n", + "hillslope_keep_fraction = 0.95\n", + "\n", + "# demo subcatchment\n", + "subcatch_demo_id = 35\n", + "\n", + "# The top level directory where these packages will go, one subdirectory for each HUC/data package\n", + "package_directory = '../huc'" + ] + }, + { + "cell_type": "markdown", + "id": "d26910d1", + "metadata": {}, + "source": [ + "## Set up directory" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e45012d2", + "metadata": {}, + "outputs": [], + "source": [ + "# directory for one huc\n", + "huc_dir = os.path.join(package_directory,f'{huc}')\n", + "if not os.path.isdir(huc_dir):\n", + " os.mkdir(huc_dir)\n", + "\n", + "# directory of raw hillslope data for this huc, downloaded from USGS\n", + "hillslope_raw_dir = os.path.join(huc_dir, 'data_raw-hillslope')\n", + "if not os.path.isdir(hillslope_raw_dir):\n", + " os.mkdir(hillslope_raw_dir)\n", + " \n", + "# directory of preprocessed hillslope data for meshing in this huc\n", + "pp_hillslope_dir = os.path.join(huc_dir, 'data_preprocessed-meshing')\n", + "if not os.path.isdir(pp_hillslope_dir):\n", + " os.mkdir(pp_hillslope_dir)\n", + " \n", + "# directory of generated meshes\n", + "mesh_dir = os.path.join(huc_dir, 'mesh')\n", + "if not os.path.isdir(mesh_dir):\n", + " os.mkdir(mesh_dir)\n", + " \n", + "# directory of raw daymet data for this huc\n", + "daymet_raw_dir = os.path.join(huc_dir, 'data_raw-daymet')\n", + "if not os.path.isdir(daymet_raw_dir):\n", + " os.mkdir(daymet_raw_dir) \n", + "\n", + "# directory of processed daymet data for this huc\n", + "daymet_dir = os.path.join(huc_dir, 'daymet')\n", + "if not os.path.isdir(daymet_dir):\n", + " os.mkdir(daymet_dir)" + ] + }, + { + "cell_type": "markdown", + "id": "1ab2fb57", + "metadata": {}, + "source": [ + "## Set up data source and watershed_workflow" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "62862906", + "metadata": {}, + "outputs": [], + "source": [ + "target_crs = workflow.crs.default_alaska_crs()\n", + "raster_extension = 'tif'\n", + "\n", + "# data sources\n", + "data_sources = dict()\n", + "data_sources['huc'] = workflow.source_list.huc_sources['WBD'] \n", + "data_sources['dem'] = workflow.source_list.dem_sources['NED 1 arc-second']\n", + "\n", + "# set up watershed_workflow\n", + "workflow.ui.setup_logging(1)\n", + "workflow.conf.rcParams['DEFAULT']['data_directory'] = hillslope_raw_dir" + ] + }, + { + "cell_type": "markdown", + "id": "78e88b15", + "metadata": {}, + "source": [ + "## Set up filenames needed in data package" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "9408461c", + "metadata": {}, + "outputs": [], + "source": [ + "filenames = hillslopes.get_filenames(huc, huc_dir)\n", + "# print(filenames)" + ] + }, + { + "cell_type": "markdown", + "id": "30bcee47", + "metadata": {}, + "source": [ + "# Acquire HUC shapefile and DEM from USGS" + ] + }, + { + "cell_type": "markdown", + "id": "7e427d29", + "metadata": {}, + "source": [ + "## Download and conversion" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "0883d586", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2021-10-01 10:56:12,437 - root - INFO: \n", + "2021-10-01 10:56:12,438 - root - INFO: Loading HUC 190604020404\n", + "2021-10-01 10:56:12,439 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:12,440 - root - INFO: \n", + "2021-10-01 10:56:12,441 - root - INFO: Loading level 12 HUCs in 190604020404\n", + "2021-10-01 10:56:12,441 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:12,443 - root - INFO: Using HUC file \"../huc/190604020404/data_raw-hillslope/hydrography/WBD_19_GDB/WBD_19.gdb\"\n", + "2021-10-01 10:56:14,729 - root - INFO: ... found 1 HUCs\n", + "2021-10-01 10:56:14,730 - root - INFO: -- 190604020404\n", + "2021-10-01 10:56:14,740 - root - INFO: Converting to shapely\n", + "2021-10-01 10:56:14,741 - root - INFO: ... found 1\n", + "2021-10-01 10:56:14,742 - root - INFO: \n", + "2021-10-01 10:56:14,742 - root - INFO: Loading Raster\n", + "2021-10-01 10:56:14,743 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:14,743 - root - INFO: Collecting raster\n", + "2021-10-01 10:56:14,777 - root - INFO: Collecting DEMs to tile bounds: [-149.3337383807433, 68.15835009523077, -148.79217986384458, 68.33619896061079]\n", + "2021-10-01 10:56:14,778 - root - INFO: Need:\n", + "2021-10-01 10:56:14,778 - root - INFO: ../huc/190604020404/data_raw-hillslope/dem/USGS_NED_1as_n69_w150.tif\n", + "2021-10-01 10:56:14,779 - root - INFO: ../huc/190604020404/data_raw-hillslope/dem/USGS_NED_1as_n69_w149.tif\n", + "2021-10-01 10:56:14,779 - root - INFO: source files already exist!\n", + "2021-10-01 10:56:14,871 - root - INFO: ... got raster of shape: (641, 1950)\n", + "2021-10-01 10:56:14,880 - root - INFO: Masking to shape\n", + "2021-10-01 10:56:14,912 - root - INFO: shape bounds: (-149.32373838074332, 68.16835009523078, -148.80217986384457, 68.32619896061078)\n", + "2021-10-01 10:56:14,917 - root - INFO: casting mask of dtype: float32 to: nan\n", + "2021-10-01 10:56:14,919 - root - INFO: ... got raster bounds: (-149.3337383807433, 68.33619896061079, -148.79207171405866, 68.15814340504932)\n" + ] + } + ], + "source": [ + "# download (if necessary) the HUC shapefile\n", + "huc_crs, huc_shape = workflow.get_huc(data_sources['huc'], huc)\n", + "\n", + "# download (if necessary) the DEM\n", + "dem_profile, dem_raster = workflow.get_raster_on_shape(data_sources['dem'], \n", + " huc_shape, huc_crs, \n", + " mask=True, nodata=np.nan)\n", + "# convert rasterio crs to workfkow standard crs\n", + "native_crs = workflow.crs.from_rasterio(dem_profile['crs'])\n", + "\n", + "# project the shapefile into the native CRS\n", + "huc_shape = workflow.warp.shply(huc_shape, huc_crs, native_crs)" + ] + }, + { + "cell_type": "markdown", + "id": "26f51878", + "metadata": {}, + "source": [ + "## Validate: plot DEM and HUC shape" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "98defa77", + "metadata": {}, + "outputs": [], + "source": [ + "fig_property = dict()\n", + "fig_property['figsize'] = (12,4)\n", + "fig_property['dpi'] = 80\n", + "fig, ax = workflow.plot.get_ax(native_crs, **fig_property)\n", + "workflow.plot.dem(dem_profile, dem_raster, ax=ax)\n", + "workflow.plot.shply([huc_shape,], native_crs, ax=ax)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "a95bc135", + "metadata": {}, + "source": [ + "## Save DEM and HUC shape files do disk" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8474823f", + "metadata": {}, + "outputs": [], + "source": [ + "if not os.path.isfile(filenames['huc']):\n", + " hillslopes.save_shapefile(filenames['huc'], [huc_shape,], native_crs) \n", + "\n", + "if not os.path.isfile(filenames['dem']):\n", + " hillslopes.save_demfile(filenames['dem'], dem_profile, dem_raster)" + ] + }, + { + "cell_type": "markdown", + "id": "49f89c9c", + "metadata": {}, + "source": [ + "# Delineate subcatchments and flowpaths" + ] + }, + { + "cell_type": "markdown", + "id": "4267d9da", + "metadata": {}, + "source": [ + "## Subcatchments and flowpaths source" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "54da22c4", + "metadata": {}, + "outputs": [], + "source": [ + "# NOTE: this needs to be added by Jon!\n", + "#\n", + "# At this point, we need:\n", + "assert(os.path.isfile(filenames['subcatchments'])) # subcatchment shapefile\n", + "assert(os.path.isfile(filenames['streams_raster'])) # streams raster\n", + "assert(os.path.isfile(filenames['aspect'])) # aspect raster\n", + "assert(os.path.isfile(filenames['slope'])) # slope raster\n", + "assert(os.path.isfile(filenames['flowpath_length'])) # raster of each pixel's distance to the stream\n", + "assert(os.path.isfile(filenames['elev_above_streams'])) # raster of HAND" + ] + }, + { + "cell_type": "markdown", + "id": "89ad20b5", + "metadata": {}, + "source": [ + "## Load all subcatchments and flowpaths from source" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "0c41298b", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2021-10-01 10:56:14,981 - root - INFO: \n", + "2021-10-01 10:56:14,983 - root - INFO: Loading shapes\n", + "2021-10-01 10:56:14,985 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:14,987 - root - INFO: Loading file: '../huc/190604020404/data_preprocessed-meshing/huc_190604020404_subcatchments.shp'\n", + "2021-10-01 10:56:15,001 - root - INFO: ... found 44 shapes\n", + "2021-10-01 10:56:15,002 - root - INFO: Converting to shapely\n", + "2021-10-01 10:56:15,023 - root - INFO: Converting to requested CRS\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Numbers of subcatchments in basin-190604020404: 44\n" + ] + } + ], + "source": [ + "_, subcatchments = workflow.get_shapes(filenames['subcatchments'], out_crs=native_crs)\n", + "print(f'Numbers of subcatchments in basin-{huc}: ', len(subcatchments))" + ] + }, + { + "cell_type": "markdown", + "id": "78b1a50e", + "metadata": {}, + "source": [ + "## Validate: plot subcatchments and flowpaths" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "e5a14591", + "metadata": {}, + "outputs": [], + "source": [ + "fig_property = dict()\n", + "fig_property['figsize'] = (12,4)\n", + "fig_property['dpi'] = 80\n", + "fig, ax = workflow.plot.get_ax(native_crs, **fig_property)\n", + "workflow.plot.dem(dem_profile, dem_raster, ax=ax)\n", + "workflow.plot.shply(workflow.utils.flatten(subcatchments),\n", + " native_crs, ax=ax, color='r')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c62d4e99", + "metadata": {}, + "source": [ + "# Land cover" + ] + }, + { + "cell_type": "markdown", + "id": "d3bdb7f1", + "metadata": {}, + "source": [ + "## Project NNSI land cover and save to disk" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "cf5331c8", + "metadata": {}, + "outputs": [], + "source": [ + "# if not os.path.isfile(filenames['land_cover']):\n", + "# landcover.reprojectLandCover(dem_profile, nssiImg_filename, filenames['land_cover'])\n", + "# Comment because no NSSI land cover map data" + ] + }, + { + "cell_type": "markdown", + "id": "c2c6a9c8", + "metadata": {}, + "source": [ + "## Validate: plot land cover " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "09835626", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "lc_profile, lc_raster = workflow.get_raster_on_shape(filenames['land_cover'],\n", + " huc_shape, huc_crs, \n", + " mask=True, nodata=np.nan)\n", + "fig_property = dict()\n", + "fig_property['figsize'] = (12,4)\n", + "fig_property['dpi'] = 80\n", + "\n", + "fig, ax = workflow.plot.get_ax(native_crs, **fig_property)\n", + "land_cover = workflow.plot.raster(lc_profile, lc_raster, ax=ax)\n", + "workflow.plot.shply([huc_shape,], native_crs, ax=ax)\n", + "# fig.colorbar(land_cover)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "2f8cef6b", + "metadata": {}, + "source": [ + "# Determine hillslopes geometry based on flowpaths" + ] + }, + { + "cell_type": "markdown", + "id": "4e2163b7", + "metadata": {}, + "source": [ + "**Objective:** \n", + "* Determine a single hillslope profile geometry for each subcatchment (44 in total).\n", + "\n", + "**Objects to obtain:** \n", + "* hillslope length\n", + "* an elevation profile along that length\n", + "* a width along that length\n", + "\n", + "**Main idea:**\n", + "* Regard each subcatchment as a single flowpath, with average properties.\n", + "\n", + "**Method:** \n", + "* Route the surface flow and generate a standard D8[1] flowpath direction vector for each pixel of the (smoothed and filled) DEM.\n", + "* Form rasters comprising of \"length along the flowpath to the stream network\" and the corresponding \"height above the stream network.\"\n", + "* Bin the rasters according to the flowpath length.\n", + "* Average pixels in each bin to determine:\n", + " - hillslope length = 90th % of the maximum flowpath length\n", + " - bins as a function of flowpath length\n", + " - elevation as a function of bin\n", + " - number of pixels in each bin gives an area\n", + "\n", + "\n", + "\n", + "**Notes** \n", + "[1] There are eight valid output directions relating to the eight adjacent cells into which flow could travel. This approach is commonly referred to as an [eight-direction (D8) flow model](https://pro.arcgis.com/en/pro-app/latest/tool-reference/raster-analysis/flow-direction.htm). This method models flow direction from each cell to its steepest downslope neighbor. The output is an integer raster whose values range from 1 to 255. The number denotes the flowpath." + ] + }, + { + "cell_type": "markdown", + "id": "25d16b4d", + "metadata": {}, + "source": [ + "## Get hillslope parameters for one subcatchment " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "cbceae4f", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2021-10-01 10:56:15,081 - root - INFO: \n", + "2021-10-01 10:56:15,081 - root - INFO: Loading shapes\n", + "2021-10-01 10:56:15,082 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:15,083 - root - INFO: Loading file: '../huc/190604020404/data_preprocessed-meshing/huc_190604020404_subcatchments.shp'\n", + "2021-10-01 10:56:15,087 - root - INFO: ... found 44 shapes\n", + "2021-10-01 10:56:15,088 - root - INFO: Converting to shapely\n", + "2021-10-01 10:56:15,141 - root - INFO: \n", + "2021-10-01 10:56:15,142 - root - INFO: Loading Raster\n", + "2021-10-01 10:56:15,143 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:15,143 - root - INFO: Loading file: '../huc/190604020404/data_preprocessed-meshing/huc_190604020404_flowpath_lengths.tif'\n", + "2021-10-01 10:56:15,144 - root - INFO: Collecting raster\n", + "2021-10-01 10:56:15,190 - root - INFO: ... got raster of shape: (130, 305)\n", + "2021-10-01 10:56:15,200 - root - INFO: Masking to shape\n", + "2021-10-01 10:56:15,232 - root - INFO: shape bounds: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n", + "2021-10-01 10:56:15,235 - root - INFO: casting mask of dtype: float32 to: -32768.0\n", + "2021-10-01 10:56:15,235 - root - INFO: ... got raster bounds: (-149.21512726962825, 68.20453229393975, -149.13040504740323, 68.16842118282744)\n", + "2021-10-01 10:56:15,239 - root - INFO: \n", + "2021-10-01 10:56:15,240 - root - INFO: Loading Raster\n", + "2021-10-01 10:56:15,240 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:15,241 - root - INFO: Loading file: '../huc/190604020404/data_preprocessed-meshing/huc_190604020404_elev_above_streams.tif'\n", + "2021-10-01 10:56:15,242 - root - INFO: Collecting raster\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bounds in my_crs: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2021-10-01 10:56:15,289 - root - INFO: ... got raster of shape: (130, 305)\n", + "2021-10-01 10:56:15,298 - root - INFO: Masking to shape\n", + "2021-10-01 10:56:15,328 - root - INFO: shape bounds: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n", + "2021-10-01 10:56:15,330 - root - INFO: casting mask of dtype: float32 to: -999.0\n", + "2021-10-01 10:56:15,331 - root - INFO: ... got raster bounds: (-149.21512726962825, 68.20453229393975, -149.13040504740323, 68.16842118282744)\n", + "2021-10-01 10:56:15,335 - root - INFO: \n", + "2021-10-01 10:56:15,335 - root - INFO: Loading Raster\n", + "2021-10-01 10:56:15,336 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:15,337 - root - INFO: Loading file: '../huc/190604020404/data_preprocessed-meshing/huc_190604020404_aspect.tif'\n", + "2021-10-01 10:56:15,339 - root - INFO: Collecting raster\n", + "2021-10-01 10:56:15,388 - root - INFO: ... got raster of shape: (130, 305)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bounds in my_crs: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n", + "bounds in my_crs: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2021-10-01 10:56:15,398 - root - INFO: Masking to shape\n", + "2021-10-01 10:56:15,434 - root - INFO: shape bounds: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n", + "2021-10-01 10:56:15,437 - root - INFO: casting mask of dtype: float32 to: -9999.0\n", + "2021-10-01 10:56:15,438 - root - INFO: ... got raster bounds: (-149.21512726962825, 68.20453229393975, -149.13040504740323, 68.16842118282744)\n", + "2021-10-01 10:56:15,442 - root - INFO: \n", + "2021-10-01 10:56:15,443 - root - INFO: Loading Raster\n", + "2021-10-01 10:56:15,444 - root - INFO: ------------------------------\n", + "2021-10-01 10:56:15,444 - root - INFO: Loading file: '../huc/190604020404/data_preprocessed-meshing/huc_190604020404_landcover.tif'\n", + "2021-10-01 10:56:15,445 - root - INFO: Collecting raster\n", + "2021-10-01 10:56:15,493 - root - INFO: ... got raster of shape: (130, 305)\n", + "2021-10-01 10:56:15,503 - root - INFO: Masking to shape\n", + "2021-10-01 10:56:15,536 - root - INFO: shape bounds: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n", + "2021-10-01 10:56:15,538 - root - INFO: casting mask of dtype: float32 to: 255.0\n", + "2021-10-01 10:56:15,539 - root - INFO: ... got raster bounds: (-149.21512726962825, 68.20453229393975, -149.13040504740323, 68.16842118282744)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bounds in my_crs: (-149.21512726962825, 68.16842118282744, -149.13040504740323, 68.20425451616197)\n", + "error: 0.0002331351861357689\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/3hg/opt/anaconda3/envs/watershed_workflow/lib/python3.9/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems\n", + " proj_string = self.to_proj4()\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADnCAYAAABbohl9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnfElEQVR4nO3deVxUVR/H8c8M+6LIDOKGC4q541rumbnv9qTllktappZmmZVZmba4lD6pmcujaeSehlmmWbnvKIoruSuIiAOissOc54+rZiWIzMCZGc779fIFzsC9XwZ+c7dzf0cnhEBRlAfTyw6gKLZMFYii5EAViKLkQBWIouRAFYii5MA5pyfbt28vrl+/XlBZFEWKgwcPbhJCtH/QczkWyPXr1wkLC8ufVIpiI3Q6nV92z6ldLEXJgSoQRcmBKhBFyYEqEEXJgSoQRcmBKhBFyYEqEEXJQY7XQZTC68YNiI6GpCR4/HHQ6WQnkkMViI0QAtLTwc1NdhJNzZrg5QVxcRASAp06aY9nZEB8PJhMEBAARYvKzZnfVIHYguRk3L2ccNab+bPRAMr4JoOHB3h6av/ufv6Qx66netPz7UCuxLkw5aM0uvd0AVfXB779m82QmPjXH7vJ9PfPr1+Hgwfhww9h+HDQ67XHk5PB11dbbLNmsHKlhNerAOlyuqOwQYMGQg01KQDJyRTzSsezqAsDOsbxWZ+jkJKi/TXe/Xj/5/c9JpJTSL5txnTLlW3XazAq5m0GeywnPLUaLfkDU5YvJmd/4nV+mHR+mDAQby5GQmYRvJzSMLrdwuiehME9BaNXKkbvVAxFMqlYIon+Tc9yS+9DZII/huJOGP2dKGp0Qe/tyZ5IAx2Hl6dn1zTmf5WpFamzfb7f6nS6g0KIBg98ThWIbahUPJEBhvXMuNaPnTu1d+vy5bV/8fEwdepf7+73v9vHx2vv7kajtkvUsSOMGAGff6690xuKZWEsko7RKw2DRwpGj2SMbrfxdb6FS0byw4swm+dEUjLbE2rROvY7lrgNxZTujQk/4l386dH4Ck9unSj7Jc01VSB2oE5tM2knzuEWWIp0Zy9cXbXjkX37tN2Yjz+G117TCsFg0D7e/dzDQ17uUaMgNhaMBoHR18y5PzO5tWEHC5d54NetqbxgjyCnAkEIke2/+vXrC6VgnDolBAhxrnEfIYQQly4J4eEhRIUKQnh5CTFkiOSAuXT0qBCBfjdFTa9zsqPkGhAmsqkBdR3ERgQFwROPm/nqZEvSNm6hbFk4ehR+/1073Tp/vuyEuVOzJhw87s6xpEB+n3lcdhyLqQKxEU5OMP59PV+nDKR/PzMIQaVKULEi+PjY13UIX38XBjc9xdovzsmOYjFVIDakSxfYtk3H7zfq0bXRNTZvlp0o755+qRLx182wf7/sKBZRBWJj6j/hxPw3Iil1diczZwrOn5edKG9KBLiwIrkbZ99ZIDuKRVSB2BidDv4zpSEjiy8nLjKBnj1lJ8qbp5+GenXN7Aj3hkOHZMfJM1Ugtkino8YXL7JW9ywHD8KmTbIDPTqdDp5soSe+xTMwaZLsOHmmCsRWdehAaWMaw1r9yZYtssPkTcmS8MX+ZtT9aSLnNpySHSdPVIHYKp0OPvmE2hEh2sGuHRo9Glau0lOsvA+fj7xITIzsRI9OFYgta9kSQwln4g9flJ0kT+4OaBwzxZ+wi/78d/x17G0yAVUgNs748rP8Ee5L2s002VHyrNOz7rzb6zzTFhno1092mkejCsTG1e1XkzSdO98N2yk7ikWemdOGvT7t2bwxk2nTsJstiSoQG+frC+++HM+o5Y0o4S+4aJ97W1CkCPVGt+CjyksZOxYSEmQHyh1VIHbg3VmlOdV+NI95R/Pee3DhguxEeeP8+qsMO/Mm5ctkMHiw7DS5owrEDjg5QcC0UUxKeI3LFzJZtkx2ojzy8YERIwip/QW7dsG4cbIDPZwqEHtRowZPdfKis9dW4uNlh7HAqFE03fM5cyZcY9o0OHBAdqCcqQKxJxMmcGjrTeKi7PeMFgYD+mFD6RHxAT17akNSDh+WHSp7qkDsSVAQ67I60dNs550SRo+G1atZNuUyjRppdyRaJDmZecOP0L/yHuZU+a92W7CVqAKxI/XqQUqWG3U3T4Vr12THyTs/PxgyBKZMwWgkb7uM0dEwbx507owoUZIJCwPwLO7F5HPPa3eaWYkqEDsSHQ2RkVCmX0uYPFl2HMu8+SYsW4bR/TYmUy6+XgitD9GECVC/PgQHw/btLK3wHh4ZiVxNN3LOOxhD0Qyr7rOpArEjpUrB7dtop38WL9Yqxl75+8PAgRhO7sq+QFJS4KefYOhQrUtdnz7aCzB9OrERsax7bimzwhozbJgOIbS7L4e2OgtHjlgtpn02MiqkfH1hwAA4erSUtovyyScwZ47sWHn31lsY587iXK1mgJf22JUrWlGsXw/btmlbiy5dtC3OY4/d+9Y3+8GpU1C9Ojz/vPaYiwuYilSwaoGoriZ25OJFIVxdhZg9W4jMq3FCGAxCnLOf7iEP8m2bb0XfwF1CTJggRP36Qvj6CtG7txDLlgkRH3/v606f1r7ktdeE6NNH6wDz7bd/Lcds1h4LmXtbCG9vIbKycp0B1dXEMQQEaP2xxo+H47F+Woe4ifbToO1BSr3YgZ8v16LazFfwP7ubGePiYNky6N1b22TeMX8+7NgBlSpB+/bw88/Qq9dfy7nb1MLDz0trFnbOOg0j1C6WHdHr4a23YO9e7YzWTyvG0v7rQO3IvUoV2fHypNXzfvxaCTw9i7Btm3Zi6mocPPWU9vndTpInT8LMmVrzvOyMHasV0rN16mi7WUFBFudTBWKH1qzRNh4fTPWm7etvoJ8wAZYvlx0rT3Q6bXoF0E5CxMTA1q1aN8kqVbQt5t0ukiVK5Lysbt1g+3agYW2tQJ591uJ8ahfLTo0bpw3TONR8FGzZAhERsiNZzGDQbl9/7TVITYUOHaBFC60ZXalS2hY0J35+2tb1mG9zq53qVb157VjfvrBuHazvv5qWV5ZCaKjsSFIJAa1awWMlEpm7O5jc3huQU29etQWxY0uXansRI7Y8S/uNr5O0/aDsSFLpdNqWdd4KHy6ZvKxy04kqEDv3+ecwc5aek171mfvSQVJTZSeSq3VreOIJuFKxmVWuh6gCsXPFi2t/FGPGe/DV6TasWGS9gXr2qmhROGZ4UhWI8pfXRjvTzbAT0xk7uZc1H7VrBy9t66cKRPk7ox/EX06SHUO6N97QPq7eWtziZakCcSDGkq6YrtjxzVRWotfD+LHprLrYUJuW15JlWSmTYgMMZb0wxWXJjmETmrdy5YZ7CW2UgQVUgTgQY2BR4m+oXyloV95NrqUsPg5Rr6YDMVYxYrrtJjuGTTAaIfxGRVLDjlm0HFUgDiT0cCDnU0rJjmETSpWCop4ZfPZjdYuWowrEgVy45sE4l2lw86bsKNK5uUGFsmbiYzMtWo4qEAdSvLgOva9PrscgOaqsLO0liIh0Y1h7y+awUwXiQAICYOy1Mdw8ESU7ilQzZkCdYDPBTscI+GSYRctSBeJAXnsNXPSZ/PyLHc0ZbWUnTmg3Vg0r/zNH3l1J0SqWHZOpAnEgTk7wVosDvLO2oewo0ixcCN4uqQyN+Ui7/dJCqkAczKje17h0y5cffpCdpOBlZsL06TDKfT7lJwzSRi1aSBWIg/GvXYqxJZawe7fsJAXv7u0fA9Lmw8svW2WZqkAcTfnyVL4dbt8d4PNo4f/MPOZ2AdfPPtKaZFmBatrgaPz9MWZcxRSbSWH79Rb78wAtfGLQ9/iP1ZaptiCORqfDWNIFU0zhGtWblpjKjlUxXK3c/K8mWVZQuN5iCglDgCfxUfY5t3peLX95C8uSu7P6desuV21BHJCxog+mBCfZMQpM+NZEhq9qQd8uifToYd1lqwJxQIYqxYlPcrWbqZYtYTbDd2MO85jRxMzFPlZfvioQB+QWVBZXXaY2VYKDWzT1OtMPtmDoW0UxGKy/fFUgDmhndCBJWe6kFYLj9LUzo3ip3kGGvW39rQeog3SH9POyRNoFncHPz/LmzbZMhB3kl5j6bP8m/1odqS2Io0lNZU54I9r1yof9DRuRlAQ1awo8GtYCoGo9z3xblyoQB3NmwR/cFEVp38dxC8RshrOnzSRUehyRkUlxy7v7ZEsViIM5vuQg1cvcoFo12Unyj7dHFlkZZnSffAzO+XuUoArEkcTEYIzcjU9AEdlJ8o3ZDK1qXSNDuJDWpnO+r08ViCP57ju2V3mJi5cd9yLhmaMpbDlVijOrw/Eplv83hqkCcRRCwOLF/JrViqFDZYfJH2lpUKWOB8+U3ke5bnULZJ2qQBzFgQNcSfJh22Ef2rWTHSZ/HN6WCMDKzQZrjWZ/KFUgjmLxYvY1GQ1AQwe94/b3JVF0LH4Al+qVC2ydqkAcQWoqrFzJJn17+veXHSb/GBLOUiagYNeprqQ7gh9/hLp18atQBDcH7jy640hRkgPKFOg61RbEESxeDIMGaQ2bTbLD5J9NV4PpN6Bgz9CpArF3V65ocx8/8wzFikFcnOxA+WPvhnhMZgMNOvoX6HpVgdi7kBBtqltPT957D1xdZQfKHxuXxdOwyAnKVyjYpniqQOzZnWsfDBxIUhLExMAwyzpt2qwN2zzpGmxZn928UAViz/bv1zo1N2nC1q3g4QG1a8sOZX3nzsGBqNK0aVPwLVVVgdizO1sPdDoMBggOxiHPYu3bB6WdYnm8R/kCX7c6zWuvUlNh1So4fBgAgwGHbRa3PCSDdmyCKn0KfN1qC2Kv1q2D+vWhbFkAhz3Fm5AA639xoWvFo/k+tP1BVIHYq2++0Xav7vjhB8fcgty+DW7OmXRreUvK+lWB2KPoaO0AvXv3ew+FhsKYMdIS5RujETALqFswo3f/SRWIPQoJgR49wFO7F9tshg0btMshjsbTE3Qii+Rq9aWsXxWIvbl77WPQoHsP/fGH9rFWLTmR8lP6rTRShTumUjWlrF8ViL3Zt08rkkaN7j20bp229fDykpgrn6z6MgZnMvALcJeyflUg9ua+ax93zZ4NnfP/9mwplq50olNAxN29yQKnroPYk5QU7dpHRMS9h06c0D62aiUpUz47ddGD+X0uAOoYRHmY0FB4/HFtvmfg1i1YsUI7wXPncojDCdRdRF+5krT1qy2IPbnv4FwIqFxZG4o1bpzcWPkmKwtD0mVMhjbSIqgCsRdRURAWBqGhREZqe1qxsdqV5mLFZIfLJ6dPk+ZaBFOqvLMPahfLXoSEQM+e4OHBrFmwc6d2Md0nf5qa24QbO4/xU0oratSQl0EViD0Q4t7Qktmz4auvtOuE/ziZ5XA8ThzEWZ9F8+byMqgCsQd79oBez7XAhsybB++8A4MHyw6V/9yOhuHmKrglZxgWoArEPty59jHmLR1RUdCrF+gd/TcnBISHSx+lrA7SbV1yMvGrfmN0m2OEfA8LFzrmXYP/EhUFTk4Y/Z0wmSAwUE4MVSC2LjSUvUH92LLPkw0bHPeC4L+Eh0O9ehgzdVKH8Tv6htruHZy5i04HJ9KiBXTo4LhdS/4lPBzq1pW+i6UKxJZdusRvEf7UDjYTEiI7TAG7r0AmT9aG9MugCsSWhYSw3qcvXbsVwl/ToUNQty4ffqgNPUtMlBOjEL7ydkIIXpoaxK6rQbRvLztMATOZtIqoWBF/f6hYUd5ulioQG2XeuZsVtzuzdo2gSRPZaQpYeLh2qk6vJzZW64ultiDK3+z/fDu3zV40a+7Al8qzc+f4A2D3bu2+9Hr15ERRBWKLkpNZuakYAaWy8nWKY5t1X4F8+ik0bSpvSI0qEBsUt/hn1um78/Y4x52MM0d3CmTTJm0A8yuvyIuiLhTaoAEflCNBb6R1a9lJJEhKgosXoXp1RvfWhtW0bCkvjioQG3Pr+CV+MTVk7Yp0qlaVnUaCiAioXp3Icy6cPAnLl4O7nH4NgNrFsjnhM7YC0OXZwnLJ/B/Cw8mqXY85c7TTu7LHnaktiA2JPCVosbA//2kZj7OzQXYcOQ4d4nfPbsycCd9+KzuM2oLYlG2LzlDK+RprfveVHUWe8HCWnWtI7drwwguyw6gtiE35ZeUteja8CbqCnYfPZmRkEHvCxIZLxZk+Q3YYjSoQG/F6jyhCL9Vj0xQHbNGeWydO0J1QdHodTZvKDqNRu1g24NT603y71ovv3j1O216F9NgDWPjlbfam1mHpUnk3SP2TKhDZoqMZ2tNEubKCp1+T2L7DBkxYXZ23n9rH00/LTvIXtYsl0Qdvp7FrThTb0xqxaQGUKiU7kTyHDkHUbV+GDMqyqfvtbShK4SLS0pn2ObzY6CSHDgrayGseaBPSUsw8oQ8jqLNtXR1VBSKDEPzYfg6pZjc6Ln+BuvV0Dt3fKjeMyZcx6f202UhtiCoQGcaPZ93xSnTukIWvXyEdkPgPxsuHMen8ZMf4F1UgBW3uXMSq1WzUd6TPC6o47ioWuY+bmZ5kZspO8neqQArApUvw9NNQs/xNSo/ojvvFSGJinWjYUHYy2+F05BBFPTO5cUN2kr9TBZKP0tPh6lVtirSY8yksv9WFAz9e5UaiDiG0wXgKIAQb9hq4keRqU2ewQJ3mzRdvvAELFkBqqnbM6eaUwaCbc6i14i3oVEd2PNsTE8P61Db06CEwGGzrbIXV6/VsWAIpR/609mLtwrlz2nTMixdrLULT0yH2WByXvKrz0RfejjuRoKXCw1lj7k7HjrZVHGDlLYgQEPS4L6O9VjP9RkVwLlwbqJYttd2mdu2geXPQpSRrRfH88zB0qOx4NuvIhmjiMnylTnOQHatuQSKOCABmJL1MzNI/rLlomxYers3VcemSNo3H8uVQqnimdr9o1aowaZLsiDbt9x2u1A+MJyhIdpJ/s1qBCAHLZ8bR3G0/ZXyTWDH1krUWbfNCQ7Wp0H77DSpUQHsxRo7UDkIWLHDsWW6s4NvIJ2jdWsiO8UBWK5D582HKN/483zSKce8788aJIRAdba3F2zSjUZtl9l7n9cmTtYZO339fiLpN5030iUSOpFencz/bvEnMagWydi2MMn7HiIklGDrSDYDFIw9Za/E2bf9+iIu785+QEJg3TztaL1pUai57cCj0Et76ZJo9aWPnd++wWqpff4U+Yik0boyTE0x+NYpJPwaTlSGpLXcB2rBBOw7nt99gzBjtgdKlZceyC1+HeNO6/GnZMbJltQKp4HsDv7b17s0N1u/tMpzLLM/GKUestQqbtGePdvzRwPME9OkDq1dD9eqyY9kNv7Rouj91Q3aMbFmtQIzpMZga/3Wev0yAjiFNTtDnoyrWWoXNSU6GqVOhRaNUyr3UDmbPhieflB3LriSbkon3rSQ7RrasUyAxMRjSrmIKbPC3h2etLcPNTE9O7HTM+6xDQ+GnnwTvXBkFo0fDc8/JjmRXUuJTWHOzLbXb2G6TCusUSGgoxjJuxN9y+dvD7iV86Fb2EO8Pu26V1diauXPMtCu6l/b/8dTGlyiPxP30UZzJoHkr2z3TZ5UCEWt/4Ki5xgMnORn/oRNrjz1G2AHbPM+dV+ZMM7t3C96suQm++EJ2HPsUHk4WTlIn6XwYywskIQHTnj85ftmHbt3+/XSDF4PpWWQjU8Y61lZkS9//kSWcCF72TiGYtDx/rPpej0CPh4fsJNmz+Dd7c/Um6pkPAFCu3AO+QKfj3SFxfL+1OKmplq7NRsyaxepfixJcIxNjGYmdle2cuHCR556Os+nLRRYXSNSq3VxOKZ7jFFl1x3XAiUyGvZhm6erkW7uWhE+/5mePHox6o3ANxrSqzEy8o04Sk2qbV9DvsqxAkpOZt6Ma9Wpn5vwu4OfHxhaTWbLCldhYi9Yo165d8MordCu1n/hEZ5scfWo3IiOZrBuHZxHbfpOxqEAyN/zK/zIH0m/Aw3/IJ99pQjnnK6xYbqcH65GRxHYfypvN97Mj3JvQUKhcWXYo+7X52xh2p9Zl+HDZSR5CCJHtv/r164uc3Oz1kgAhkpNz/DJNVpaY7/euACFSUnLx9bYkJkaIwECxZPA2UbeuEEuWCJGRITuUfatT/LLoUfOkSEqSnUQIIExkUwMWbUHm7KiFt0dm7gas6vW8NNobF30mixdbstYCdvs2dO5Mk8xtDFj4JI0bQ//+he5eMKvq2xcOxwXwztAEPD1lp8mZRQXiUrkCQ4IP4JTb7jWDBvGx6yQWL8yyZLUFJyODQ23f4UfDQPZcLkt0NHz1lexQ9i0rC1atEvzi1YP6z9nuEJO7LCoQQ8PKmM48wlWeUqVo3+gGSTGSZoV/FEJwod94mu37nAWuw3n1VfC33RERdmP7dsjM1NGw6Em7eEEtKhBjwyDiE514lFNTxoFdOBZt4MABS9ac/1b3+p6hG7uT5ezG+p/0zJqldqusoWRJ7ePxQPtoYGFRgfiXdubnzPac+mZPrr+nTL+WdHffyKqv4h7+xbIsWsSMdZWo80Iw+/er22WtqVo16FvzCHs9bWiOgxxYVCBPPAGNK8Uyd9EjDDZzcqJT6zTi95/511MpKRAVpTU/kGbjRk6PXcCetHr0HuIlfZZVR1Q94zCmohVkx8gViwpEp4P3PnDiy9MduXE9901V/bo05pdTFWgZdJngmmYCAsDDA3x9oWFDrRHIxYuWJMujQ4egf392Dg2hdGmoU0dChkLAcPUk8W72ccelxUNNOvX3o6j+FpNfv5rr7+ky2J+5M1J5v/RCQmJas2fAXEwXb5OSovV5CA6GHj20P9C1ay1NmEsXLkCXLjBvHj+eCHrgwEvFCmJjMWbGYkr3lp0kd7K7QCJycaHwrlXdlwoXfYZITc3DVZqICCGee04If38hPvtMiJs3xdmzQmzeLMRbbwnRqpUQs2YJcfNmHpadWyaTEFWqiMj3Q0RgoBAgxIYN+bi+wuyXX8Tvdd4QTz0lO8hfyK8LhXe1fbkCLmSwYEEevrlWLVi5Ev74AyIioFIlKi7/hNaPJzJ0KAQEaNceVq/W2k1ZXUoKdO2KuUs3vsnoh5cXJCZChw75sC4FwsMJqOvH1q2wbZvsMA9nlQLxafMEb7nOJO7crbwvpEYNWLZMO1F+6hQEBVFp6UQW//cGAwfC4MEwcaI10t4nKwv69YNy5dj01GdMngwjRqhuPfkqPJzHni7LwIGwd6/sMA9nnTt9nJ0xVi3Orz+mWP4uX7Wq1ltq9244fx6Cghib9CFL591m6VIYNAgiI62QWQgO9P0vo/f3or8uhOGv6mnRAl55xQrLVrIXHg716lG1Kg+8A9XWWO1WuLbPFWPvWX/277fSAitX1hrd7tuH7ko03d6pxmsV1hNzKZ2vvybH+09yZcYMvt5UkegG3WjV1omZM2HpUqskV7Jz8ybExECVKhiN9lEgVjlIF0IIERMjerqsFUZDlli/3qrHUJrz54V4+WXxk/fzomyRBNH7P3kbEnzokBA1yyaIkvpYAUKsXm3dmEoOtm8XomFDIYQQa9YI0b275Dx3kN8H6QCULMk3Y07Q37yEL6ekWH//skIFmDePTsensqzVQjb8kEYd/ytsWJqQ60VkZcGiSdH4X40gbJOJ5GTtdLJSQMLDoW5dALvZgli124DXp+8xbLgOn4NbGPd6Etevg9nanUfLlaPp2jf55Yc0WpT6k1kDDzLjqVDeH32LCxdy/tbl06KY/UMZ+owsTpnW1Wy6WYBDulMgFy9C167aHpfNy27TIh51F+s+Z6auEUV0N4Wba5aYMsXyTWB2Tp4Uon+PJDGqzlbxpPNO8W79TUJERT34a7fFitJOMaJXo/P5F0jJUVZwHXF6zRExdaoQRqMQ8fGyE2nIYRcrXwpECCHMoevEXK83RFCZJDF3rmU/QG6s+V+8ACGu+FQVYvhwIS5duvdcVkKi+LjkLBFkNInLl/M/i/IAqalikcvLolgxs6hZU4hPP5Ud6C9SCkQIIZI2bhefek0SQaVuic2btQvW+emxx4R4d+RtIcaOFcJgEGLoUDF5TJxw0mUKEGLaVHP+BlCyFxYmuhT5Q7z4ouwg/yatQIQQInFbuGjluk08VuKGePllixeXo+++0zbd4eFC/L72hlj1zDLR0mmr+KLGQiEyM/N35UrOFiwQOrLEli2yg/xbTgWS7y0Biz5Zh98iSvBfMYrlS9IoW1abRiM/dO2qnewaMAAmzfRhpb43lQc0pfOaQeT+vmAlP/z2YzICvd3NDFEw98hVqUKHsEnsa/EfZhZ7n+nTG1KsmI4GDR7+rY+iSBEIC/vno+o2QNkSE2Hy9sZ0bmzC398oO84jKbimsmXLUm3/Ekamf4HLiSNM+FBw/Xo+DUBUbMrqlWb2JlbjrfFusqM8soLtuuznR7U9i/is5Ey2bEqjTBmhuoQUAjOmZtDV6w+e7Ggn94Dcp+DbkhcpQvWtc7jdpQ+TK8zjq1lmvvmmwFMoBSA1FTZvhhNn3Rj9xC7ZcfJETt9+d3d0q1cxpOFRnk9axPRpmfz5p5QkSj6aP9fMS71v0cd9LVWes8+b++VNbOHsTJEls3mt+2X8Lh6idrAgy076ySkPd2DeIca/mUwnl80s3Vmeoq/0kR0pT+TO/KLTYZw1gS3v/4FP5nWCq6YxYoTURIqlLlwg/dnefDfmMFUqpPNRxDNQv77sVHkm/xyoTgfvvMNe8zKOfbGJF5f/j8BAF9q0QbXcsSe3b8PkySR+9R0VUk9yM92d+eN0+BWXHcwyNjN3WIVxfWg3uwsvpn7NljXxTJoE6emyUykPZTbDt99C1aqI8xf49o1wbqR6kJWlY/Bg2eEsJ38Lch+X3j2YbPyV8Od70OTQZpo0cXrAhT/FZuzZA6+/rn3+/fdcLt2IkeVh4UKpqawruzEowkpjsfJkzx6R6FdRuDpnir59hVi6VE4MJRuXLwvRt68QpUsLsWSJOHcmSzRrJkTZskJ4eMgO9+iQORYrTxo1ouiWdawoOpSA2INMnw5n/t2pVCloyclaa5natbVBb5GRJHTpz5of9CQkwKZN2i3njsSmdrH+pmZNngl7j8Yte7ObldSvX47ERNVIWgohYNUqGDtWa8gcFgaBgQA0ra4dnw8dqjWmdjS2WyAAgYGU3BvKtnbt8bmyl1q13GjbVscXX8gOVogcPKgdZ9y+rR2Mt2iBEPDtEvjzTzh5UmvD9NhjsoPmD9vcxbpfyZLotm3lcHB/JhhmsnKl4PvvZYcqBK5ehRdfhM6dtfsHwsKgRQtAG507ZIjWcHzxYseezNT2CwSgWDEqbl9MJ48t9HILpW9fQUICXL8uO5gDSkuDKVOgZk2t9cipU1o1ODlx5Qo0awb16kFmJrz3nlY7Ogfe87WPAgHw9MR9/Wo+b/Q9wa6n8PMTlC6tdStVrEAI+OEHqF5d62q5Zw9MmwY+PiQna83vQ0O1uVs2btS2Io5cGHfZ9jHIP7m4QEgIu9+bgFPIYj5ss5vJkwMICdF+YbNna+9uyiOKiNCOM65dg7lzoU2be08JAUFB2g2ZXl4wcKDjHm88iH0VCIBej8tnE6HJ47w+qBWVO8zG+Hxrvlms48svtd9t797qDttciYuDDz6ANWvgww+1U1F3JmL880+tj7jJpJ26jY/XJjgqbOxnF+ufunTBuG8D/Y++RadlfRn5UgopKfDmm7Bvn+xwNi4jA/77X213ytVVO84YMeJvs5SOGqVNOXH9OsybB8WKSUsrV3ZXEIXMK+mPIilJiAEDhKhRQ4jISPHCC0I4OQmxapXsYDbq55+FqFJFiHbthDh+/N7DWVlCTJggxPPPC9G6tdbv5pdfJOYsQMhs+1MgzGYh5s0Tonhxkb5yrRg/XojGjYXYuFF2MBty4oQQ7dtrzcN++kl7ze6TlKT9NSxdqr1uYWH/+hKH5fgFctf+/UKULy+iXpkkXhyYJZo2zf9mdTYvPl6IkSOF8PMTYvp0IdLS/vb0q69qvcScnITw9i6c7cNyKhD7PQZ5kMcfh7AwypzbwYenenPlciY1asgOJUlmJsyZo01IlJYGJ04gXh9NUoYrly7B4cPw++/ahb6FC7WZ6G7dUic3/sn+zmI9jJ8fbNhAuY8+4uzCCnjEX+K55/TMnQsGg+xwBePIvL3MfjcKk1Mt4iucxbTLG1Nt7UyUTqdd/zMYtI9t20LjxtoZdOXfHK9AQHsbnDgRXcOGbHuhE4N3hLBooZEmTXU0aSI7XD46cwbGjGHF1nZcDerOgLdLYvTT3SsGoxE15cOjym7fS9jjMciDnD0rZpX5TLQuESFcXMzi9m3ZgfJBYqI4M+QzEaQ/K4q6pQgQBdJR31FQaI5BHqRiRV49PYrNHWdQUXeewHKZPPOM7FDWkZGaxdlPV3Kg4nMs3l8d54rluHjVnaws7ZqfYjnHLxDQ9isWLWL31F38Ktry+6YMZszAvtsMbd/OxAqLaPRhW4b5r2GXsSvDRzpTrBjoC8dvtUA45jFINgyjXqBYk+pMbPcZE8eNoVkjV4KqOtvFEIpTp+DKFYg/dQ3TwlDizyYwX4zizYluvPNuIRg1KEmhKhAA/eP1ef10BXbXDKPDk8EUKeXN+Uu2/TKcOQP16gka+p/HcOU4xuD6GAYHM6akCy+8IDudY7Ptv4z8YjSyKroZmRM+xmPSu7RvmEC7Xr6MHi072AOYzWz5dB9OqbXY0vR9mDwZypaVnarQKLx7q3o9zhM/YNvnYbQ/MZ25n8WzfJng6lXZwe6zdy+ZjZuzfm0Gg3vchKVLVXEUsMJbIHc0ebMxg3YNoXbGQaaNuMDHH6QjhHYhWproaHjhBZ5rGYd72A7WJz5J55dLSwxUeBX6AgHwCS7PqpjmfNzwRxb8D9xchZwOHSkpMGkSBAcT5tKY79M6s3WbHiGgdWsJeRRVIPe4u9Nx4yiiZnzPTZ+yXDhvZsgQWLSoANYtBKxcqY2bOnIEwsLotmk4HTroVH9iyQrnQXoOio/qA82rsrrDMA4de5Yp29tQubIOk0kby1ShAjz9tBVXeOiQdnfSrVuYv1nCiNVPcXyAdkr311+1eRcVeXQih0kCGzRoIMIKa3Pc+HgSew2l846xmIv4YCyux7O4F1tP+BMVrSc+QUd8vHZLqskEwcFa8eTa1ataW5Cff4aJE8kaOJircU4EBGijbMuUgSpV8uuHU+6n0+kOCiEeOKWs2oJkx2DAZ+NKdoSHw5GdEBFB2uGTlL3+He6uvvi63MbgnY7RAJnuXvgFuLNhoxNZWZCQwL0tzt0CqlZNa0pIWhp8+SVMnap1QIiMBB8fxr2t3dpauzY0b65G19oKVSA50eu1yV/uTADjBlxOFbhcj0F/LEI7XoiI4Oj+FJ7YtBRfp3RuCW983NMx+mZhKOGCsZQrer2OCxcERz5ax403J2EKbIDp63Divcpi+lEroD/+0Nrejhwp90dW/k7tYlnJlfNpuF8+jc+5cJyOaYVDRASkpRFXtTmVDiwnOcuNot5CK5z77sm4+3HoUChZUvZPUvjktIulCiS/xcbC0aMkX72JW8+uOLmpjbatUccgMpUoASVK4Ck7h5In6jqIouRAFYii5EAViKLkQBWIouRAFYii5EAViKLkQBWIouRAFYii5CDHK+k6nU5dRlcKg+tCiPYPeiLHAlGUwk7tYilKDlSBKEoOVIEoSg5UgShKDlSBKEoO/g+qA0YbFEkmwgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "hillslope_demo_pars = hillslopes.parameterizeSubcatchment(\n", + " filenames, huc, subcatch_demo_id,\n", + " target_crs=target_crs,\n", + " hillslope_keep_fraction=hillslope_keep_fraction,\n", + " hillslope_bin_dx=hillslope_bin_dx,\n", + " plot_smoothed=True)\n", + "\n", + "mesh_demo_pars = meshing.parameterizeMesh(hillslope_demo_pars, mesh_dx,\n", + " toeslope_min_slope=toeslope_min_slope,\n", + " hillslope_min_slope=hillslope_min_slope,\n", + " min_area_ratio=min_area_ratio)\n", + "\n", + "# print(hillslope.keys())" + ] + }, + { + "cell_type": "markdown", + "id": "187c4e5d", + "metadata": {}, + "source": [ + "## Validate: plot the hillslope" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "8d9e358a", + "metadata": {}, + "outputs": [], + "source": [ + "# fig = plt.figure(figsize=(6,8))\n", + "# plt.tight_layout()\n", + "# gs = gridspec.GridSpec(4,1)\n", + "# axs = [fig.add_subplot(gs[2,0]), fig.add_subplot(gs[3,0])]\n", + "# plot.plot(hillslope_demo_pars, mesh_demo_pars, fig=fig, axs=axs)\n", + "\n", + "\n", + "# ax0 = workflow.plot.get_ax(native_crs, fig, axgrid=gs[0,:])\n", + "# workflow.plot.raster(hillslope_demo_pars['raster_profile'], hillslope_demo_pars['bins'],\n", + "# ax=ax0, cmap='prism')\n", + "# subcatch_demo = hillslope_demo_pars['subcatchment']\n", + "# subcatch_crs_demo = hillslope_demo_pars['subcatchment_target_crs']\n", + "# workflow.plot.shply(subcatch_demo, subcatch_crs_demo, ax=ax0, color='k')\n", + "# ax0.set_title(f'subcatchment {subcatch_demo_id}')\n", + "\n", + "# ax1 = workflow.plot.get_ax(native_crs, fig, axgrid=gs[1,:])\n", + "# vmin, vmax = min(hillslope_demo_pars['land_cover']), max(hillslope_demo_pars['land_cover'])\n", + "# cmap = plt.get_cmap(('viridis'),vmax-vmin+1)\n", + "# lc = workflow.plot.raster(hillslope_demo_pars['raster_profile'],\n", + "# hillslope_demo_pars['land_cover_raster'], \n", + "# ax=ax1, vmin=vmin-0.5, vmax=vmax+0.5, cmap=cmap)\n", + "# workflow.plot.shply(subcatch_demo, subcatch_crs_demo, ax=ax1, color='k')\n", + "# position=fig.add_axes([0.95, 0.55, 0.02,0.1])\n", + "\n", + "# fig.colorbar(lc,ticks=np.arange(vmin,vmax+1),cax=position)\n", + "\n", + "# plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "1f4cb183", + "metadata": {}, + "source": [ + "## Get parameters for all hillslopes" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "44f418ed", + "metadata": {}, + "outputs": [], + "source": [ + "# hillslope_pars = []\n", + "# mesh_pars = []\n", + "# subcatch_all_ids = 1+np.arange(len(subcatchments))\n", + "# for subcatch_id in subcatch_all_ids:\n", + "# hillslope_pars.append(hillslopes.parameterizeSubcatchment(\n", + "# filenames, huc, subcatch_id,\n", + "# target_crs=target_crs,\n", + "# hillslope_keep_fraction=hillslope_keep_fraction,\n", + "# hillslope_bin_dx=hillslope_bin_dx))\n", + " \n", + "# for i in range(len(hillslope_pars)): \n", + "# mesh_pars.append(meshing.parameterizeMesh(hillslope_pars[i], mesh_dx,\n", + "# toeslope_min_slope=toeslope_min_slope,\n", + "# hillslope_min_slope=hillslope_min_slope,\n", + "# min_area_ratio=min_area_ratio))" + ] + }, + { + "cell_type": "markdown", + "id": "d89ac989", + "metadata": {}, + "source": [ + "## Plot all hillslopes" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "789b2ed3", + "metadata": {}, + "outputs": [], + "source": [ + "# fig = plt.figure(figsize=(20,30),dpi=80)\n", + "# nx, ny = 6,8\n", + "# sep=0.02\n", + "\n", + "# axs = []\n", + "# for i in range(nx):\n", + "# for j in range(ny):\n", + "# sub_id = i*ny+j+1\n", + "# if sub_id > len(subcatchments):\n", + "# continue\n", + " \n", + "# gs = gridspec.GridSpec(4,1,bottom=j/ny+sep, left=i/nx+sep, \n", + "# top=(j+1)/ny - sep, right=(i+1)/nx - sep)\n", + "# ax1 = fig.add_subplot(gs[2,0])\n", + "# ax2 = fig.add_subplot(gs[3,0])\n", + "# axs = [ax1,ax2]\n", + "# h_par = hillslope_pars[sub_id-1]\n", + "# m_par = mesh_pars[sub_id-1] \n", + "# plot.plot(h_par, m_par, fig=fig, axs=axs)\n", + " \n", + "# ax0 = workflow.plot.get_ax(native_crs, fig, axgrid=gs[0,0])\n", + "# ax0.set_title(f'subcatchment {sub_id}')\n", + "# workflow.plot.raster(hillslope_pars[sub_id-1]['raster_profile'], \n", + "# hillslope_pars[sub_id-1]['bins'], ax=ax0, cmap='prism')\n", + "# subcatch = hillslope_pars[sub_id-1]['subcatchment']\n", + "# subcatch_crs = hillslope_pars[sub_id-1]['subcatchment_target_crs'] \n", + "# workflow.plot.shply(subcatch, subcatch_crs, ax=ax0, color='k')\n", + " \n", + "# ax1 = workflow.plot.get_ax(native_crs, fig, axgrid=gs[1,0])\n", + "# vmin = min(hillslope_pars[sub_id-1]['land_cover']) \n", + "# vmax = max(hillslope_pars[sub_id-1]['land_cover'])\n", + "# cmap = plt.get_cmap(('viridis'),vmax-vmin+1)\n", + "# lc = workflow.plot.raster(hillslope_pars[sub_id-1]['raster_profile'], \n", + "# hillslope_pars[sub_id-1]['land_cover_raster'], \n", + "# ax=ax1, vmin=vmin-0.5, vmax=vmax+0.5, cmap=cmap)\n", + "# workflow.plot.shply(subcatch, subcatch_crs, ax=ax1, color='k')\n", + "# position=fig.add_axes([i/nx+sep+0.13, j/ny+sep+0.045, 0.006,0.04])\n", + "\n", + "# fig.colorbar(lc,ticks=np.arange(vmin,vmax+1),cax=position)\n", + "\n", + "# plt.tight_layout()\n", + "# plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d60fb0f2", + "metadata": {}, + "source": [ + "# Generate mesh" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "4916162d", + "metadata": {}, + "outputs": [], + "source": [ + "layer_info = meshing.layeringStructure(organic_cells=30, organic_cell_dz=0.02, \n", + " increase2depth=9.4, increase_cells=20, largest_dz=2.0,\n", + " bottom_depth=46)\n", + "layer_types, layer_data, layer_ncells = layer_info" + ] + }, + { + "cell_type": "markdown", + "id": "366f43ba", + "metadata": {}, + "source": [ + "## Column spinup mesh" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "d75f2a06", + "metadata": {}, + "outputs": [], + "source": [ + "# if os.path.isfile(os.path.join(mesh_dir,'column.exo')):\n", + "# os.remove(os.path.join(mesh_dir,'column.exo'))\n", + "# colum_mesh = meshing.createColumnMesh(layer_info, \n", + "# os.path.join(mesh_dir,'column.exo'))" + ] + }, + { + "cell_type": "markdown", + "id": "973ad6c6", + "metadata": {}, + "source": [ + "## 2D for a given subcatchment" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "4396c703", + "metadata": {}, + "outputs": [], + "source": [ + "# m2 = meshing.createHillslopeMesh2D(mesh_demo_pars)\n", + "# m3 = meshing.createHillslopeMesh3D(m2, mesh_demo_pars, layer_info,\n", + "# os.path.join(mesh_dir,f'sag_hillslope{subcatch_demo_id}.exo'))" + ] + }, + { + "cell_type": "markdown", + "id": "889672f6", + "metadata": {}, + "source": [ + "## 3D for a given subcatchment" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "4bf13c05", + "metadata": {}, + "outputs": [], + "source": [ + "subcatch_demo = hillslope_demo_pars['subcatchment']\n", + "subcatch_demo_smooth = hillslope_demo_pars['subcatchment_smooth']\n", + "subcatch_crs_demo = hillslope_demo_pars['subcatchment_target_crs']\n", + "\n", + "nodem, m2, lc = meshing.createSubcatchmentMesh2D(filenames, subcatch_demo_smooth,\n", + " subcatch_crs_demo, mesh_demo_pars,\n", + " refine_max_area = refine_max_area, plot=True)\n", + "if nodem > 0:\n", + " _, m2, lc = meshing.createSubcatchmentMesh2D(filenames, subcatch_demo,\n", + " subcatch_crs_demo, mesh_demo_pars,\n", + " refine_max_area = refine_max_area, plot=True)\n", + "\n", + "m3 = meshing.createSubcatchmentMesh3D(m2, lc, layer_info,\n", + " os.path.join(mesh_dir, \n", + " f'sag_subcatchment{subcatch_demo_id}.exo'))" + ] + }, + { + "cell_type": "markdown", + "id": "d121163b", + "metadata": {}, + "source": [ + "## 2D for all subcatchments" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "457de366", + "metadata": {}, + "outputs": [], + "source": [ + "# subcatch_all_ids = np.arange(1, len(subcatchments)+1)\n", + "# for subcatch_id in subcatch_all_ids:\n", + "# index = subcatch_id - 1\n", + "# m2 = meshing.createHillslopeMesh2D(mesh_pars[index])\n", + "# m3 = meshing.createHillslopeMesh3D(m2, mesh_pars[index], meshing.layeringStructure(),\n", + "# os.path.join(mesh_dir,f'sag_hillslope{subcatch_id}.exo'))" + ] + }, + { + "cell_type": "markdown", + "id": "5e846eff", + "metadata": {}, + "source": [ + "## 3D for all subcatchments" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "d32ea0e5", + "metadata": {}, + "outputs": [], + "source": [ + "# give_up_id = []\n", + "# must_smooth_id = []\n", + "# no_smooth_id = []\n", + "\n", + "# subcatch_all_ids = np.arange(1, len(subcatchments)+1)\n", + "# for subcatch_id in subcatch_all_ids:\n", + "# index = subcatch_id - 1\n", + "# subcatch = hillslope_pars[index]['subcatchment']\n", + "# subcatch_smooth = hillslope_pars[index]['subcatchment_smooth']\n", + "# subcatch_crs = hillslope_pars[index]['subcatchment_target_crs']\n", + " \n", + "# if type(subcatch) == shapely.geometry.multipolygon.MultiPolygon:\n", + "# if subcatch_smooth is None:\n", + "# give_up_id.append(subcatch_id)\n", + " \n", + "# else:\n", + "# must_smooth_id.append(subcatch_id)\n", + "# _, m2, lc = meshing.createSubcatchmentMesh2D(filenames, subcatch_smooth, subcatch_crs,\n", + "# mesh_pars[index], plot=False)\n", + "# m3 = meshing.createSubcatchmentMesh3D(m2, lc, layer_info, \n", + "# os.path.join(mesh_dir, f'sag_subcatchment{subcatch_id}.exo'))\n", + " \n", + "# else:\n", + "# nodem, m2, lc = meshing.createSubcatchmentMesh2D(filenames, subcatch_smooth, \n", + "# subcatch_crs, mesh_pars[index],\n", + "# plot=False)\n", + "# if nodem > 0:\n", + "# no_smooth_id.append(subcatch_id)\n", + "# _, m2, lc = meshing.createSubcatchmentMesh2D(filenames, subcatch, subcatch_crs,\n", + "# mesh_pars[index], plot=False)\n", + " \n", + " \n", + "# m3 = meshing.createSubcatchmentMesh3D(m2, lc, layer_info, \n", + "# os.path.join(mesh_dir, f'sag_subcatchment{subcatch_id}.exo'))\n", + " \n", + "# print(give_up_id)\n", + "# print(must_smooth_id)" + ] + }, + { + "cell_type": "markdown", + "id": "bd6e1541", + "metadata": {}, + "source": [ + "# DayMet" + ] + }, + { + "cell_type": "markdown", + "id": "42957ce2", + "metadata": {}, + "source": [ + "## For a given subcatchment" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "f40ce9ae", + "metadata": {}, + "outputs": [], + "source": [ + "# start, end = datetime.date(1980,1,1), datetime.date(2020,12,31)\n", + "# hillslopes.downloadDaymet(hillslope_demo_pars, daymet_raw_dir, \n", + "# os.path.join(daymet_dir, f'huc_{huc}_subcatchment{subcatch_demo_id}_'\n", + "# +str(start.year)+'_'+str(end.year)+'.h5'),\n", + "# start, end)" + ] + }, + { + "cell_type": "markdown", + "id": "ddf3aa46", + "metadata": {}, + "source": [ + "## For all subcatchment" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "b2e67606", + "metadata": {}, + "outputs": [], + "source": [ + "# start, end = datetime.date(2011,1,1), datetime.date(2020,12,31)\n", + "# subcatch_all_ids = np.arange(1, len(subcatchments)+1)\n", + "# for subcatch_id in subcatch_all_ids:\n", + "# index = subcatch_id - 1\n", + "# hillslopes.downloadDaymet(hillslope_pars[index], daymet_raw_dir, \n", + "# os.path.join(daymet_dir, f'huc_{huc}_subcatchment{subcatch_id}_'\n", + "# +str(start.year)+'_'+str(end.year)+'.h5'),\n", + "# start, end)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:watershed_workflow]", + "language": "python", + "name": "conda-env-watershed_workflow-py" + }, + "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.9.4" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "586px", + "left": "455px", + "top": "691.125px", + "width": "288px" + }, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/new_scripts/landcover.py b/new_scripts/landcover.py new file mode 100644 index 0000000..059a232 --- /dev/null +++ b/new_scripts/landcover.py @@ -0,0 +1,118 @@ +import numpy as np +import rasterio +import attr + +_lc_ids = [1,2,3,5,6,7,8,9,10,11,14,16,18,19,20,21,22,23,50,51,91,92,95,96,99] +_lc_names = ['Bare Ground','Sparsely Vegetated','Open Water','FWM: Arctophila Fulva', + 'FWM: Carex Aquatillis','Wet Sedge','Wet Sedge - Sphagnum','Mesic Herbaceous','Tussock Tundra', + 'Tussock Shrub Tundra','Mesic Sedge-Dwarf Shrub Tundra','Dwarf Shrub - Dryas','Dwarf Shrub - Other', + 'Birch Ericaceous Low Shrub','Low-Tall Willow','Alder','Marine Beach/Beach Meadow','Coastal Marsh', + 'Ice / Snow','Burned Area','Open Needleleaf','Woodland Needleleaf','Open Mixed Needleleaf/Deciduous', + 'Deciduous','Unclassified (cloud, terrain shadow, etc.)'] + +name_to_id = dict(zip(_lc_names, _lc_ids)) +id_to_name = dict(zip(_lc_ids, _lc_names)) + +class_names, class_ids = dict(), dict() +class_names['shrub'] = ['Birch Ericaceous Low Shrub','Low-Tall Willow','Alder'] +class_names['sedge'] = ['FWM: Carex Aquatillis', 'Mesic Sedge-Dwarf Shrub Tundra', 'Wet Sedge','Wet Sedge - Sphagnum'] +class_names['tussock'] = ['Tussock Tundra','Tussock Shrub Tundra', 'Dwarf Shrub - Dryas'] +class_names['sparse_veg'] = ['Bare Ground','Sparsely Vegetated','Dwarf Shrub - Other','Open Water','Ice / Snow'] +for group, names in class_names.items(): + class_ids[group] = [name_to_id[name] for name in names] + + +def vegClasses(): + return list(class_ids.keys()) + + +# Remaps from NSSI IDs to lumped ATS IDs around shrub,sedge,tussock,and sparse_veg +def classifyVegetation(lc): + lc_remap = np.zeros(lc.shape, dtype=np.uint8) + lc_remap[lc == 255] = 255 + for i, (group,ids) in enumerate(class_ids.items()): + ats_id = 100+i + for veg_id in ids: + lc_remap[lc == veg_id] = ats_id + missing_id = set(lc[[lc_remap == 0]]) + for mid in missing_id: + print(f'Missing vegetation id {mid} type: ', id_to_name[mid]) + + return len(missing_id), lc_remap + + +# Return a float version with NaNs to improve plotting +def veg2img(lc): + new_lc = np.nan * np.ones(lc.shape, 'd') + for i in set(lc.ravel()): + new_lc[lc == i] = i + new_lc[lc == 255] = np.nan + return new_lc + + +# Reprojects land cover from the NSSI default CRS into (ugh) Lat/Long the CRS of the bands +def reprojectLandCover(dem_profile, nssiImg_filename, lc_filename): + rio_profile = dem_profile.copy() + rio_profile.pop('blockxsize') + rio_profile.pop('blockysize') + rio_profile.pop('tiled') + rio_profile['nodata'] = 255 + rio_profile['driver'] = 'GTiff' + + with rasterio.open(nnsiImg_filename, 'r') as fin: + with rasterio.open(lc_filename, 'w', **rio_profile) as fout: + rasterio.warp.reproject( + source=rasterio.band(fin, 1), + destination=rasterio.band(fout, 1), + src_transform=fin.transform, + src_crs=fin.crs, + dst_transform=rio_profile['transform'], + dst_crs=rio_profile['crs'], + resampling=rasterio.enums.Resampling.nearest) + + + + +# Determine soil structure according to land cover +@attr.s +class SoilHorizons: + acrotelm = attr.ib() + catotelm = attr.ib() + +soil_horizons = dict() +soil_horizons['riparian sparse_veg'] = SoilHorizons(acrotelm=0, catotelm=0) +soil_horizons['hillslope sparse_veg'] = SoilHorizons(acrotelm=0, catotelm=0) + +soil_horizons['riparian sedge'] = SoilHorizons(acrotelm=0.10, catotelm=0.30) +soil_horizons['hillslope sedge'] = SoilHorizons(acrotelm=0.08, catotelm=0.20) + +soil_horizons['riparian shrub'] = SoilHorizons(acrotelm=0.12, catotelm=0.24) +soil_horizons['hillslope shrub'] = SoilHorizons(acrotelm=0.14, catotelm=0.14) + +soil_horizons['riparian tussock'] = SoilHorizons(acrotelm=0.10, catotelm=0.14) +soil_horizons['hillslope tussock'] = SoilHorizons(acrotelm=0.10, catotelm=0.14) + +soil_horizons['average'] = SoilHorizons(acrotelm=0.1, catotelm=0.2) + + + +# Given a soil cell thickness dz, determine the material label of each cell +def soilStructure(dz, land_cover): + ''' + 1000 = mineral + 1001 = acrotelm + 1002 = catotelm + ''' + if land_cover == -999: + lc_name = 'average' + elif land_cover < 110: + lc_name = 'hillslope ' + vegClasses()[land_cover - 100] + else: + lc_name = 'riparian ' + vegClasses()[land_cover - 110] + + horizons = soil_horizons[lc_name] + z_depth = np.cumsum(dz) + soil_type = np.where(z_depth < horizons.acrotelm, 1001, + np.where(z_depth < horizons.acrotelm + horizons.catotelm, 1002, 1000)) + + return soil_type \ No newline at end of file diff --git a/new_scripts/meshing.py b/new_scripts/meshing.py new file mode 100644 index 0000000..3ed68b2 --- /dev/null +++ b/new_scripts/meshing.py @@ -0,0 +1,287 @@ +import numpy as np +import scipy +import hillslopes +import landcover +import workflow +import os +import matplotlib.pyplot as plt + +# Given a hillslope's parameters, generate the mesh parameters +def parameterizeMesh(hillslope, dx, + toeslope_min_slope=0.0, + hillslope_min_slope=0.1, + min_area_ratio=0.1): + mesh = dict() + + # 1. Determine x coordinate + mesh['huc'] = hillslope['huc'] + mesh['subcatchment_id'] = hillslope['subcatchment_id'] + mesh['dx'] = dx + mesh['num_cells'] = int(np.round(hillslope['total_length'] / dx)) + mesh['length'] = mesh['dx'] * mesh['num_cells'] + mesh['x'] = np.linspace(0, mesh['length'], mesh['num_cells']) + + # 2. Determine z coordinate + # 2.1 Interpolate from (x,z)=(0,0), and (x,z)=(bin_centroid,bin_average) + x_bin = np.concatenate([np.array([0.,]), (np.arange(0,hillslope['num_bins']) + 0.5)*hillslope['bin_dx']]) + z_bin = np.concatenate([np.array([0.,]), hillslope['elevation']]) + z_native = np.interp(mesh['x'], x_bin, z_bin) + mesh['z_native'] = z_native + z = scipy.signal.savgol_filter(mesh['z_native'], window_length=11, polyorder=3) + + # 2.2 Determine riparian area and hillslope area according to slope + for i in range(1,int(np.round(len(z)/2))): + if z[i] < toeslope_min_slope * (mesh['x'][i] - mesh['x'][i-1]) + z[i-1]: + z[i] = toeslope_min_slope * (mesh['x'][i] - mesh['x'][i-1]) + z[i-1] + for i in range(int(np.round(len(z)/2)),len(z)): + if z[i] < hillslope_min_slope * (mesh['x'][i] - mesh['x'][i-1]) + z[i-1]: + z[i] = hillslope_min_slope * (mesh['x'][i] - mesh['x'][i-1]) + z[i-1] + mesh['z'] = scipy.signal.savgol_filter(z, 5, 3) + + slope = (mesh['z'][1:] - mesh['z'][0:-1]) / (mesh['x'][1:] - mesh['x'][0:-1]) + riparian = np.zeros(slope.shape, 'i') + i = 0 + while slope[i] < 0.1: + riparian[i] = 1 + i += 1 + mesh['riparian'] = riparian + if i == 1: + mesh['riparian_width'] = 0 + else: + mesh['riparian_width'] = (mesh['x'][i-1] + mesh['x'][i])/2. + + + # smooth z once more to deal with discontinuities + z = scipy.signal.savgol_filter(z, window_length=5, polyorder=3) + z = z - z[0] + mesh['z'] = z + + # 3. Determine y coordinate + # 3.1 Interpolate bins to create a bin-consistent profile + y = np.interp(mesh['x'], x_bin[1:], hillslope['bin_counts']) + + # 3.2 Smooth y + y = scipy.signal.savgol_filter(y, window_length=51, polyorder=3, mode='nearest') + + # 3.3 Don't let individual areas get too small -- 10% mean as a min value? + min_y = min_area_ratio * y.mean() + y = np.maximum(y, min_y) + + # 3.4 Scale by area ratios to ensure that the final mesh has the identical surface area as the subcatchment it represents + y_factor = hillslope['total_area'] / np.trapz(y, mesh['x']) + y = y_factor * y + mesh['y'] = y + + # 4. Resample land cover onto mesh + land_cover_mesh = np.zeros(len(mesh['x'])-1, 'i') + def lc_index(lc_type, is_riparian): + if is_riparian: + return lc_type + 10 + else: + return lc_type + + for i in range(len(land_cover_mesh)): + x = (mesh['x'][i] + mesh['x'][i+1])/2 + j_bin = np.where(x_bin>x,x_bin,np.inf).argmin()-1 +# j_bin = int(np.round(x / hillslope['bin_dx'] - 0.5)) # THE OLD + land_cover_mesh[i] = lc_index(hillslope['land_cover'][j_bin], riparian[i]) + mesh['land_cover'] = land_cover_mesh + + # 5. Add aspect to mesh + mesh['aspect'] = hillslope['aspect'] + + return mesh + + + +# Take mesh parameters and turn those into a 2D surface transect mesh +def createHillslopeMesh2D(mesh_pars): + labeled_sets = list() + for i,vtype in zip(range(100, 104), landcover.vegClasses()): + labeled_sets.append(workflow.mesh.LabeledSet(f'hillslope {vtype}', i, 'CELL', + [int(c) for c in np.where(mesh_pars['land_cover'] == i)[0]])) + for i,vtype in zip(range(110, 114), landcover.vegClasses()): + labeled_sets.append(workflow.mesh.LabeledSet(f'riparian {vtype}', i, 'CELL', + [int(c) for c in np.where(mesh_pars['land_cover'] == i)[0]])) + assert(min(mesh_pars['y']) > 0) + m2 = workflow.mesh.Mesh2D.from_Transect(mesh_pars['x'], mesh_pars['z'], mesh_pars['y'], + labeled_sets=labeled_sets) + rotation = 90 + mesh_pars['aspect'] # 90 makes the aspect due north + m2.transform(mat=workflow.mesh.transform_rotation(np.radians(rotation))) + return m2 + + +# Preparing layer extrusion data +# +# Meshes are extruded in the vertical by "layer", where a layer may consist of multiple cells in the z direction. +# These layers are logical unit to make construction easier, +# and may or may not correspond to material type (organic/mineral soil). +# +# The extrusion process is then given four lists, each of length num_layers. +def layeringStructure(organic_cells=30, organic_cell_dz=0.02, + increase2depth=9.4, increase_cells=20, largest_dz=2.0, + bottom_depth=46): + layer_types = [] # a list of strings that tell the extruding code how to do the layers. + # See meshing_ats documentation for more, but here we will use only "constant", + # which means that dz within the layer is constant. + + layer_data = [] # this data depends upon the layer type, but for constant is the thickness of the layer + + layer_ncells = [] # number of cells (in the vertical) in the layer. + # The dz of each cell is the layer thickness / number of cells. + + # 30 layers at 2cm makes 60cm of cells, covering the deepest organic layer and getting into mineral soil + n_top = organic_cells + dz = organic_cell_dz + current_depth = 0 + for i in range(n_top): + layer_types.append('constant') + layer_data.append(dz) + layer_ncells.append(1) + + # telescope from 2cm to 2m + dzs, res = workflow.mesh.optimize_dzs(organic_cell_dz, largest_dz, increase2depth, increase_cells) + for dz in dzs: + layer_types.append('constant') + layer_data.append(dz) + layer_ncells.append(1) + + num_at_2m = int(np.round((bottom_depth - sum(layer_data)) / largest_dz)) + for i in range(num_at_2m): + layer_types.append('constant') + layer_data.append(2) + layer_ncells.append(1) + + return layer_types, layer_data, layer_ncells + + + +# Exstrude the 2D hillslope surface mesh downward to construct 3D hillslope mesh +def createHillslopeMesh3D(m2, mesh_pars, layer_info, mesh_fname): + if os.path.isfile(mesh_fname): + os.remove(mesh_fname) + + layer_types, layer_data, layer_ncells = layer_info + layer_mat_ids_near_surface = np.array([landcover.soilStructure(layer_data, mesh_pars['land_cover'][c]) + for c in range(m2.num_cells())]).transpose() + layer_mat_ids = list(layer_mat_ids_near_surface) + m3 = workflow.mesh.Mesh3D.extruded_Mesh2D(m2, layer_types, layer_data, layer_ncells, layer_mat_ids) + m3.write_exodus(mesh_fname) + + + +# Create column mesh for spinup +def createColumnMesh(layer_info, filename): + m2 = workflow.mesh.Mesh2D.from_Transect(np.array([-0.5, 0.5]), np.array([0., 0.]), 1.0) + + layer_types, layer_data, layer_ncells = layer_info + layer_mat_ids = list(landcover.soilStructure(layer_data, -999)) + + m3 = workflow.mesh.Mesh3D.extruded_Mesh2D(m2, layer_types, layer_data, layer_ncells, layer_mat_ids) + m3.write_exodus(filename) + + +# Smooth DEM data +def smoothDEM(dem): + dem1 = dem.copy() + nodata = 0 + for i in range(len(dem1)): + if len(list(set(dem1[i]))) == 1 and list(set(dem1[i]))[0] == -9999: + nodata += 1 + + else: + if dem1[i][0] == -9999: + index = (dem1[i]!=-9999).argmax() + dem1[i][:index] = dem1[i][index]+0.5 + + if dem1[i][-1] == -9999: + index = (dem1[i][::-1]!=-9999).argmax() + dem1[i][-index:] = dem1[i][-index-1]-0.5 + + if -9999 in dem1[i]: + index1 = np.where(dem1[i] == -9999)[0][0] + index2 = np.where(dem1[i] == -9999)[0][-1] + d = (dem1[i][index2+1] - dem1[i][index1-1])/(index2-index1+3) + for j in range(index1, index2+1): + dem1[i][j] = dem1[i][index1-1]+(j-index1+1)*d + + return nodata, dem1 + +# Take mesh parameters and create a 2D subcatchment surface mesh +def createSubcatchmentMesh2D(filenames, subcatch, subcatch_crs, mesh_pars, + refine_max_area = 200, tol = 1, plot=False): + # 1. triangulate the subcatchment + verts, tris, areas, dist = workflow.triangulate([subcatch,], list(), tol = tol, + refine_max_area = refine_max_area) + centroids = np.array([verts[t].mean(0) for t in tris]) + + # 2. elevate the triangulation + dem_profile, dem0 = workflow.get_raster_on_shape(filenames['dem'], subcatch, + subcatch_crs, mask=False) + + nodem, dem = smoothDEM(dem0) + + dem_crs = workflow.crs.from_rasterio(dem_profile['crs']) + verts3 = workflow.elevate(verts, subcatch_crs, dem, dem_profile) + + # 3. get a land cover raster + lc_profile, lc_raster = workflow.get_raster_on_shape(filenames['land_cover'],subcatch, + subcatch_crs,mask=False) + + lc = workflow.values_from_raster(centroids, subcatch_crs, lc_raster, lc_profile) + lc = lc.astype(int) + + # 4. renumber from NSSI ids to ATS IDs + _, lc = landcover.classifyVegetation(lc) + + # 5. riparian vs hillslope + # 5.1 load raster of flowpath length + fpl_profile, fpl_raster = workflow.get_raster_on_shape(filenames['flowpath_length'],subcatch, + subcatch_crs,mask=False) + + + fpl = workflow.values_from_raster(centroids, subcatch_crs, fpl_raster, fpl_profile) + riparian = np.where(fpl > mesh_pars['riparian_width'], 1, 0) + + # 5.2 riparians are incremented by 10 for unique indices + lc = np.where(riparian, 10+lc, lc) + + # 6. label land cover + labeled_sets = list() + for i,vtype in zip(range(100, 104), landcover.vegClasses()): + labeled_sets.append(workflow.mesh.LabeledSet(f'hillslope {vtype}', i, 'CELL', + [int(c) for c in np.where(lc == i)[0]])) + + for i,vtype in zip(range(110, 114), landcover.vegClasses()): + labeled_sets.append(workflow.mesh.LabeledSet(f'riparian {vtype}', i, 'CELL', + [int(c) for c in np.where(lc == i)[0]])) + + # 7. create surface mesh + m2 = workflow.mesh.Mesh2D(verts3, tris, labeled_sets=labeled_sets) + if plot: + fig = plt.figure(figsize=(6,6)) + ax = workflow.plot.get_ax('3d', fig) + cax = fig.add_axes([1.1,0.3,0.03,0.5]) + mp = ax.plot_trisurf(verts3[:,0], verts3[:,1], verts3[:,2], + triangles=tris, cmap='viridis', + edgecolor=(0,0,0,.2), linewidth=0.5) + + cb = fig.colorbar(mp, orientation="vertical", cax=cax) + t = cax.set_title('elevation [m]') + + return nodem, m2, lc + + + +# Extrude 2D subcatchment suface mesh to create 3D subcatchment mesh +def createSubcatchmentMesh3D(m2, lc, layer_info, mesh_fname): + if os.path.isfile(mesh_fname): + os.remove(mesh_fname) + + layer_types, layer_data, layer_ncells = layer_info + assert(len(lc) == m2.num_cells()) + layer_mat = np.array([landcover.soilStructure(layer_data, lc[c]) for c in range(m2.num_cells())]).transpose() + layer_mat_ids = list(layer_mat) + m3 = workflow.mesh.Mesh3D.extruded_Mesh2D(m2, layer_types, layer_data, layer_ncells, layer_mat_ids) + m3.write_exodus(mesh_fname) + diff --git a/new_scripts/plot.py b/new_scripts/plot.py new file mode 100644 index 0000000..ad315ae --- /dev/null +++ b/new_scripts/plot.py @@ -0,0 +1,20 @@ +import matplotlib.pyplot as plt + + +def plot(hillslope_pars, mesh_pars, fig=None, axs=None, color='k'): + """plotting routine""" + if fig is None: + fig = plt.figure() + if axs is None: + axs = fig.subplots(2,1, sharex=True) + fig.subplots_adjust(hspace=0) + + axs[0].plot(mesh_pars['x'], mesh_pars['z']) + axs[0].set_xticklabels([]) + axs[0].set_ylabel('elev [m]') + axs[0].yaxis.set_label_position("right") + + axs[1].plot(mesh_pars['x'], mesh_pars['y']) + axs[1].set_xlabel('hillslope length [m]') + axs[1].set_ylabel('width') + axs[1].yaxis.set_label_position("right") \ No newline at end of file