diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402.webp deleted file mode 100644 index e9e07c5..0000000 Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402.webp and /dev/null differ diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_masked.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_masked.webp deleted file mode 100644 index 2327880..0000000 Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_masked.webp and /dev/null differ diff --git a/content/tutorials/parallelization/SIMWE_parallelization.qmd b/content/tutorials/parallelization/SIMWE_parallelization.qmd index 6747649..8bc504c 100644 --- a/content/tutorials/parallelization/SIMWE_parallelization.qmd +++ b/content/tutorials/parallelization/SIMWE_parallelization.qmd @@ -9,6 +9,7 @@ description: > and erosion simulation (SIMWE) split by subwatersheds in parallel. image: SIMWE_images/thumbnail.webp +bibliography: references.bib format: ipynb: default html: @@ -41,9 +42,9 @@ subwatersheds, each approximately 100 km² in size) and run simulations independently for each unit. This approach demonstrates how to leverage new Python API features introduced -in GRASS 8.5—specifically, the [Tools API](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.tools.html#) -and the [MaskManager](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.MaskManager) -and [RegionManager](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.RegionManager) +in GRASS 8.5—specifically, the [Tools API](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.tools.html#) +and the [MaskManager](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.MaskManager) +and [RegionManager](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.RegionManager) context managers—which simplify region and mask handling and are especially useful in parallel workflows. @@ -136,7 +137,7 @@ gs.create_project(nlcd_project, filename=nlcd_filename) session = gj.init(nlcd_project) ``` -Link the NLCD raster with [r.external](https://grass.osgeo.org/grass-devel/manuals/r.external.html). This command creates a virtual raster without importing the full dataset, allowing us to work with just the small sections needed from the larger, nationwide NLCD file. +Link the NLCD raster with [r.external](https://grass.osgeo.org/grass-stable/manuals/r.external.html). This command creates a virtual raster without importing the full dataset, allowing us to work with just the small sections needed from the larger, nationwide NLCD file. ```{python} tools = Tools() @@ -219,7 +220,7 @@ tools.v_in_ogr( ``` Next, we want to extract the subwatersheds along the river. -If we simply overlap (with [v.select](https://grass.osgeo.org/grass-devel/manuals/v.select.html)) +If we simply overlap (with [v.select](https://grass.osgeo.org/grass-stable/manuals/v.select.html)) the river and the subwatersheds, we will miss some of them because the river data don't always overlap or touch the subwatersheds. @@ -239,7 +240,7 @@ basin_map.show() ![](SIMWE_images/basins.webp) -So instead we will [buffer](https://grass.osgeo.org/grass-devel/manuals/v.buffer.html) the river first: +So instead we will [buffer](https://grass.osgeo.org/grass-stable/manuals/v.buffer.html) the river first: ```{python} tools.v_buffer(input="river", output="river_buffer", distance=0.01) @@ -265,12 +266,12 @@ The rest of the workflow will be done in a CRS used in North Carolina (EPSG 3358 Since we want our project to use a different CRS more suitable for our study area (EPSG:3358 for North Carolina), we will create it now: ```{python} -nc_project = path / "NC" +nc_project = Path(path) / "NC" gs.create_project(nc_project, epsg="3358") session = gj.init(nc_project) ``` -We will reproject the river subwatersheds vector with [v.proj]((https://grass.osgeo.org/grass-devel/manuals/v.proj.html)): +We will reproject the river subwatersheds vector with [v.proj]((https://grass.osgeo.org/grass-stable/manuals/v.proj.html)): ```{python} tools.v_proj(input="river_basins", project="hydro") @@ -278,14 +279,15 @@ tools.v_proj(input="river_basins", project="hydro") ## Downloading NED elevation data -We will use GRASS addon [r.in.usgs](https://grass.osgeo.org/grass-devel/manuals/addons/r.in.usgs.html) to download and import NED layer using [USGS TNM Access API](https://apps.nationalmap.gov/tnmaccess/). +We will use GRASS addon [r.in.usgs](https://grass.osgeo.org/grass-stable/manuals/addons/r.in.usgs.html) to download and import NED layer using [USGS TNM Access API](https://apps.nationalmap.gov/tnmaccess/). ```{python} tools.g_extension(extension="r.in.usgs") ``` -First, we will set the computational region to match the river subwatersheds and then list available datasets -for that region (the output at the time of creating this tutorial, it may vary later on): +First, we will set the computational region to match the river subwatersheds +and then list the latest available datasets for that region +(the output at the time of creating this tutorial, it may vary later on): ```{python} tools.g_region(vector="river_basins") @@ -293,6 +295,7 @@ print( tools.r_in_usgs( product="ned", ned_dataset="ned13sec", + ned_release="current", flags="i", output_directory=os.getcwd(), ).text @@ -302,41 +305,33 @@ print( ```text USGS file(s) to download: ------------------------- -Total download size: 5.32 GB -Tile count: 17 +Total download size: 2.90 GB +Tile count: 6 USGS SRS: wgs84 USGS tile titles: -USGS 1/3 Arc Second n36w080 20240611 -USGS 1/3 Arc Second n36w081 20240611 -USGS 1/3 Arc Second n36w082 20220504 -USGS 1/3 Arc Second n36w082 20220512 -USGS 1/3 Arc Second n36w082 20240611 -USGS 1/3 Arc Second n37w080 20210305 -USGS 1/3 Arc Second n37w081 20210305 -USGS 1/3 Arc Second n37w081 20240611 -USGS 1/3 Arc Second n37w082 20210305 -USGS 1/3 Arc Second n37w082 20220512 -USGS 1/3 Arc Second n37w082 20240611 -USGS 13 arc-second n36w080 1 x 1 degree -USGS 13 arc-second n36w081 1 x 1 degree -USGS 13 arc-second n36w082 1 x 1 degree -USGS 13 arc-second n37w080 1 x 1 degree -USGS 13 arc-second n37w081 1 x 1 degree -USGS 13 arc-second n37w082 1 x 1 degree +USGS 1/3 Arc Second n36w080 20100929 +USGS 1/3 Arc Second n36w081 20100929 +USGS 1/3 Arc Second n36w082 20100929 +USGS 1/3 Arc Second n37w080 20100929 +USGS 1/3 Arc Second n37w081 20100929 +USGS 1/3 Arc Second n37w082 20100929 ------------------------- ``` -We will select the 2024 data and use the `title_filter` option to filter them. -We will download and import the datasets. +Now we will download and import the datasets. ```{python} tools.r_in_usgs( product="ned", ned_dataset="ned13sec", + ned_release="current", output_directory=os.getcwd(), - title_filter="20240611", output_name="ned", ) + +elevation_map = gj.Map() +elevation_map.d_rast(map="ned") +elevation_map.show() ``` Let's now compute the erosion/deposition for one of the subwatersheds. @@ -345,8 +340,8 @@ Let's now compute the erosion/deposition for one of the subwatersheds. We begin our analysis by selecting a single HUC12 subwatershed from the set of intersected basins. This subset will be used to test the erosion and deposition simulation workflow. For all outputs, we will use unique names with the subwatershed code as a prefix. That will help us with parallelization later on. -Use [v.db.select](https://grass.osgeo.org/grass-devel/manuals/v.db.select.html) to retrieve the list of HUC12 codes from the attribute table that has "areasqkm" and "huc12" attributes. -Extract the smallest HUC12 polygon using [v.extract](https://grass.osgeo.org/grass-devel/manuals/v.extract.html). +Use [v.db.select](https://grass.osgeo.org/grass-stable/manuals/v.db.select.html) to retrieve the list of HUC12 codes from the attribute table that has "areasqkm" and "huc12" attributes. +Extract the smallest HUC12 polygon using [v.extract](https://grass.osgeo.org/grass-stable/manuals/v.extract.html). ```{python} # Select for example the smallest huc12 @@ -360,122 +355,147 @@ tools.v_extract( ) ``` -Set the computational region with [g.region](https://grass.osgeo.org/grass-devel/manuals/g.region.html) to match the bounds of the selected subwatershed. +Set the computational region with [g.region](https://grass.osgeo.org/grass-stable/manuals/g.region.html) to match the extent of the basin and align resolution with the elevation raster. ```{python} -tools.g_region(vector=f"basin_{huc12}") +tools.g_region(raster="ned", vector=f"basin_{huc12}") ``` -Display the NED data for that region: +Convert the selected vector to a raster using [v.to.rast](https://grass.osgeo.org/grass-stable/manuals/v.to.rast.html). ```{python} -basin_map = gj.Map(current_region=True) -basin_map.d_rast(map="ned") -basin_map.d_vect(map=f"basin_{huc12}", fill_color="none") -basin_map.show() +tools.v_to_rast(input=f"basin_{huc12}", output=f"basin_{huc12}", use="val") +``` + +Apply a raster mask with [MaskManager](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.MaskManager) to restrict all subsequent raster operations to this subwatershed. +MaskManager is new in GRASS version 8.5. Alternatively, you can use [r.mask](https://grass.osgeo.org/grass-stable/manuals/r.mask.html) to set a mask. + +```{python} +mask = gs.MaskManager(mask_name=f"basin_{huc12}") +mask.activate() ``` -![](SIMWE_images/basin_030401010402.webp) -Set the region to match the extent of the basin and align resolution with the elevation raster. -Convert the selected vector to a raster using [v.to.rast](https://grass.osgeo.org/grass-devel/manuals/v.to.rast.html). +Reproject NLCD from a "nlcd" project to this project (see [NLCD legend](https://www.usgs.gov/media/images/annual-nlcd-land-cover-change-legend)): ```{python} -tools.g_region(raster="ned", vector=f"basin_{huc12}") -tools.v_to_rast(input=f"basin_{huc12}", output=f"basin_{huc12}", use="val") +tools.r_proj(project="nlcd", mapset="PERMANENT", input="nlcd", output=f"nlcd_{huc12}") + +nlcd_map = gj.Map() +nlcd_map.d_rast(map=f"nlcd_{huc12}") +nlcd_map.show() ``` -Apply a raster mask with [MaskManager](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.MaskManager) to restrict all subsequent raster operations to this subwatershed. -MaskManager is new in GRASS version 8.5. Alternatively, you can use [r.mask](https://grass.osgeo.org/grass-devel/manuals/r.mask.html) to set a mask. +![](SIMWE_images/basin_030401010402_nlcd.webp) + +Next, we will estimate Manning's roughness coefficient from the NLCD landcover +using [r.manning](https://grass.osgeo.org/grass-stable/manuals/addons/r.manning.html) addon tool. +Specifically, we will use [@kalyanapu2009] lookup table that is more suitable for shallow flow. ```{python} -mask = gs.MaskManager(mask_name=f"basin_{huc12}") -mask.activate() +tools.g_extension(extension="r.manning") +tools.r_manning(input=f"nlcd_{huc12}", output=f"mannings_{huc12}", landcover="nlcd", source="kalyanapu") -basin_map = gj.Map(current_region=True) -basin_map.d_rast(map="ned") -basin_map.d_legend( - raster="ned", +manning_map = gj.Map() +manning_map.d_rast(map=f"mannings_{huc12}") +manning_map.d_legend( + raster=f"mannings_{huc12}", flags="t", at=[10, 15, 40, 95], - title="Elevation [m]", + title="Roughness", fontsize=12, ) -basin_map.show() +manning_map.show() ``` -![](SIMWE_images/basin_030401010402_masked.webp) +Then we will estimate rainfall excess and to do that we will first process +soil data for our area using [r.in.ssurgo](https://grass.osgeo.org/grass-stable/manuals/addons/r.in.ssurgo.html) addon tool. +Since we need a fairly large area, we will first download the soil data itself for North Carolina +from [Gridded Soil Survey Geographic (gSSURGO) Database](https://www.nrcs.usda.gov/resources/data-and-reports/gridded-soil-survey-geographic-gssurgo-database) +and then have r.in.ssurgo extract Hydrologic Soil Group for our area. -Reproject NLCD from a "nlcd" project to this project (see [NLCD legend](https://www.usgs.gov/media/images/annual-nlcd-land-cover-change-legend)): +::: {.callout-note title="What is Hydrologic Soil Group?"} +Hydrologic Soil Group (HSG) classifies soils based on their minimum +infiltration rate and potential for rainfall runoff. +For smaller areas, r.in.ssurgo can download the data directly using the +[Soil Data Access (SDA) API](https://sdmdataaccess.nrcs.usda.gov/). +::: ```{python} -tools.r_proj(project="nlcd", mapset="PERMANENT", input="nlcd", output=f"nlcd_{huc12}") +tools.g_extension(extension="r.in.ssurgo") +tools.r_in_ssurgo(ssurgo_path="gSSURGO_NC.zip", soils=f"soils_{huc12}", hydgrp=f"hydgrp_{huc12}", nprocs=1) -basin_map = gj.Map() -basin_map.d_rast(map=f"nlcd_{huc12}") -basin_map.show() +soil_map = gj.Map() +soil_map.d_rast(map=f"hydgrp_{huc12}") +soil_map.d_legend( + raster=f"hydgrp_{huc12}", + at=[0, 30, 0, 10], + title="HSG", + fontsize=12, + flags="cn", +) +soil_map.show() ``` -![](SIMWE_images/basin_030401010402_nlcd.webp) +With the soil and landcover data available, we will generate Curve Number, +which is an empirical parameter to estimate surface runoff. -We will create a file `mannings.txt` to reclassify NLCD to manning's coefficients. The values are suggestions from the [r.sim.water](https://grass.osgeo.org/grass-devel/manuals/r.sim.water.html) documentation. The format is described in [r.recode](https://grass.osgeo.org/grass-devel/manuals/r.recode.html) documentation. +```{python} +tools.g_extension(extension="r.curvenumber") +tools.r_curvenumber(landcover=f"nlcd_{huc12}", landcover_source="nlcd", soil=f"hydgrp_{huc12}", + output=f"curvenumber_{huc12}") -```text -%%writefile mannings.txt -11:11:0.001 -21:21:0.0404 -22:23:0.0678 -24:24:0.0404 -31:31:0.0113 -41:41:0.36 -42:42:0.32 -43:52:0.4 -52:52:0.4 -71:71:0.368 -81:82:0.325 -90:90:0.086 -95:95:0.1825 +cn_map = gj.Map() +cn_map.d_rast(map=f"curvenumber_{huc12}") +cn_map.d_legend( + raster=f"curvenumber_{huc12}", + at=[10, 15, 40, 95], + flags="t", + title="CN", + fontsize=12, +) +cn_map.show() ``` +We will model a 50-year storm event that lasts 3 hours, using data from +[NOAA Atlas 14](https://hdsc.nws.noaa.gov/pfds/) Precipitation Frequency Data Server. + ```{python} -tools.r_recode(input=f"nlcd_{huc12}", output=f"mannings_{huc12}", rules="mannings.txt") +tools.g_extension(extension="r.noaa.atlas14") +duration = 3 +depths = tools.r_noaa_atlas14(mode="point", statistic="depth", units="metric", format="json") +depth = [row["values"] for row in depths["table"]["rows"] if row["duration"] == f"{duration}-hr"][0]["50"] +print(depth) ``` -Similarly, we will reclassify rainfall excess (rainfall minus infiltration). -These rainfall excess values were derived using the SCS Curve Number method for 50 mm/h of precipitation, assigning typical Curve Numbers based on common hydrologic soil groups for each NLCD land cover class; -they provide a rough estimate and should be refined with local soil and land management data when available. +We will then create an input rainfall depth raster: -```text -%%writefile runoff.txt -11:11:0.0 -21:21:3.3 # Developed Open Space (CN ~68, HSG B) -22:22:7.9 # Developed Low Intensity (CN ~81, HSG C) -23:23:19.8 # Developed Medium Intensity (CN ~90, HSG C) -24:24:27.9 # Developed High Intensity (CN ~94, HSG D) -31:31:6.1 # Barren (CN ~75, HSG B) -41:41:0.0 # Forest (CN ~55, HSG A) — P < 0.2S → Q = 0 -42:42:0.0 -43:43:0.0 -52:52:3.3 # Shrub/Scrub (CN ~68, HSG B) -71:71:2.0 # Grassland (CN ~61, HSG A) -81:81:5.1 # Pasture (CN ~74, HSG B) -82:82:14.7 # Cropland (CN ~82, HSG C) -90:90:0.0 # Woody Wetlands (CN ~70, HSG D) — often saturated, still low Q here -95:95:0.0 # Herbaceous Wetlands ``` - -```{python} -tools.r_recode(input=f"nlcd_{huc12}", output=f"runoff_{huc12}", rules="runoff.txt") +tools.r_mapcalc(expression=f"rainfall_{huc12} = {depth}") ``` -Then, compute slope components (dx and dy) from the elevation raster to support hydrologic modeling. +With that we estimate runoff (rainfall excess): ```{python} -tools.r_slope_aspect(elevation="ned", dx=f"dx_{huc12}", dy=f"dy_{huc12}") +tools.g_extension(extension="r.runoff") +tools.r_runoff(rainfall=f"rainfall_{huc12}", duration=duration, curve_number=f"curvenumber_{huc12}", + runoff_depth=f"runoff_depth_{huc12}") +tools.r_mapcalc(expression=f"runoff_{huc12} = runoff_depth_{huc12} / {duration}") + +runoff_map = gj.Map() +runoff_map.d_rast(map=f"runoff_{huc12}") +runoff_map.d_legend( + raster=f"runoff_{huc12}", + at=[10, 15, 40, 95], + flags="t", + title="Runoff [mm/h]", + fontsize=12, +) +runoff_map.show() ``` Run hydrologic and sediment simulations for the selected subwatershed. -First, simulate overland flow using [r.sim.water](https://grass.osgeo.org/grass-devel/manuals/r.sim.water.html) +First, simulate overland flow using [r.sim.water](https://grass.osgeo.org/grass-stable/manuals/r.sim.water.html) with inputs for topography, Manning’s coefficients, and rainfall intensity. Running the simulation may take a while. @@ -483,8 +503,6 @@ Running the simulation may take a while. simulation_time = 180 tools.r_sim_water( elevation="ned", - dx=f"dx_{huc12}", - dy=f"dy_{huc12}", depth=f"depth_{huc12}", niterations=simulation_time, man=f"mannings_{huc12}", @@ -505,7 +523,7 @@ basin_map.show() ![](SIMWE_images/basin_030401010402_simwe.webp) -Then, define parameters for sediment transport and run [r.sim.sediment](https://grass.osgeo.org/grass-devel/manuals/r.sim.sediment.html) +Then, define parameters for sediment transport and run [r.sim.sediment](https://grass.osgeo.org/grass-stable/manuals/r.sim.sediment.html) to compute erosion and deposition patterns. We will trim the resulting raster's edges to avoid extreme values at the edge. The additional parameters for sediment erosion modeling are based on [WEPP](https://www.ars.usda.gov/midwest-area/west-lafayette-in/national-soil-erosion-research/docs/wepp/research/) model, here we use just a single-value estimate. @@ -516,8 +534,6 @@ tools.r_mapcalc(expression=f"shear_stress_{huc12} = 0.01") tools.r_sim_sediment( elevation="ned", water_depth=f"depth_{huc12}", - dx=f"dx_{huc12}", - dy=f"dy_{huc12}", detachment_coeff=f"detachment_coef_{huc12}", transport_coeff=f"transport_capacity_{huc12}", shear_stress=f"shear_stress_{huc12}", @@ -544,9 +560,9 @@ basin_map.show() ![](SIMWE_images/basin_030401010402_erdep.webp) Calculate total erosion for the subwatershed. -Negative values from [r.sim.sediment](https://grass.osgeo.org/grass-devel/manuals/r.sim.sediment.html) +Negative values from [r.sim.sediment](https://grass.osgeo.org/grass-stable/manuals/r.sim.sediment.html) output are treated as erosion and converted to sediment mass. -Summary statistics are computed using [r.univar](https://grass.osgeo.org/grass-devel/manuals/r.univar.html). +Summary statistics are computed using [r.univar](https://grass.osgeo.org/grass-stable/manuals/r.univar.html). ```{python} # Erosion: extract negative values and convert to positive mass [kg/m^2] @@ -585,7 +601,7 @@ Each subwatershed needs to set different computational region and mask. However those setting are usually global for each mapset. So, to use different regions and masks for each parallel process, we will use the region and mask context managers. -* Computational region is handled using `RegionManager`, a [context manager for setting and managing computational region](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.RegionManager), +* Computational region is handled using `RegionManager`, a [context manager for setting and managing computational region](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.RegionManager), making it possible to have custom region for the current process. This feature is available only since GRASS 8.5. ```python @@ -594,7 +610,7 @@ making it possible to have custom region for the current process. This feature i tools.r_sim_water(...) ``` -* Masking is handled using `MaskManager`, a [context manager for setting and managing raster mask](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.MaskManager), +* Masking is handled using `MaskManager`, a [context manager for setting and managing raster mask](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.MaskManager), making it possible to have custom mask for the current process. This feature is available only since GRASS 8.5. ```python @@ -665,15 +681,8 @@ def compute(huc12): output=f"runoff_{huc12}", rules="runoff.txt", ) - tools.r_slope_aspect( - elevation="ned", - dx=f"dx_{huc12}", - dy=f"dy_{huc12}", - ) tools.r_sim_water( elevation="ned", - dx=f"dx_{huc12}", - dy=f"dy_{huc12}", depth=f"depth_{huc12}", niterations=simulation_time, man=f"mannings_{huc12}", @@ -686,8 +695,6 @@ def compute(huc12): tools.r_sim_sediment( elevation="ned", water_depth=f"depth_{huc12}", - dx=f"dx_{huc12}", - dy=f"dy_{huc12}", detachment_coeff=f"detachment_coef_{huc12}", transport_coeff=f"transport_capacity_{huc12}", shear_stress=f"shear_stress_{huc12}", diff --git a/content/tutorials/parallelization/references.bib b/content/tutorials/parallelization/references.bib new file mode 100644 index 0000000..e69de29