Author

Corey T. White

Modified

April 2, 2026

[WIP] SSURGO Soil Data in GRASS

This notebook demonstrates how to import SSURGO soil data into GRASS using the r.in.ssurgo extension. The r.in.ssurgo tool can import both local SSURGO data and data from the Soil Data Access (SDA) web service. The imported data includes soil areas, hydrologic groups, ksat maps, and mukey maps, which can be used for hydrological modeling and analysis in GRASS.

The flowchart below (@fig_ssurgo_flowchart) illustrates the workflow for using SSURGO data in GRASS, including the tools and outputs that can be generated from the imported soil data. We will cover the following steps in this notebook:

%% label: fig_ssurgo_flowchart
%% echo: false
%% eval: true
%% fig-cap: "(Draft) Flowchart of SSURGO SIMWE workflows"
flowchart TD
    A[r.in.ssurgo] --> B("Soils map vector")
    A[r.in.ssurgo] --> HG("Hydrologic group map")
    A[r.in.ssurgo] --> KSAT("Ksat maps rasters")
    A[r.in.ssurgo] --> E("MUKEY map raster")
    B --> F("Soil attributes in vector attribute table")
    G[r.curvenumber] --> H("Curve number map")
    HG --> G
    LC(NLCD land cover map) --> G
    ELEV[Elevation raster] --> W[r.watershed]
    W --> I("Flow direction")
    W --> J("Streams")
    W --> K("Basins")
    TC[r.timeofconcentration] --> TCR("Time of concentration")
    ELEV --> TC
    I --> TC
    J --> TC
    H --> R[r.runoff]
    I --> R
    TCR --> R
    R --> RD("Runoff depth")
    R --> RV("Runoff volume")
    R --> TP("Time to peak")
    R --> PD("Peak discharge raster")
    R --> URD("Upstream runoff depth")
    R --> URV("Upstream runoff volume")
    R --> UA("Upstream area")
    LC --> MC["r.manning"]
    MC --> MN("Manning's n map")
    ELEV --> RSA["r.slope.aspect"]
    RSA --> DX("dx")
    RSA --> DY("dy")
    ELEV --> SIMWE["r.sim.water"]
    DX --> SIMWE
    DY --> SIMWE
    MN --> SIMWE
    HG --> SIMWE
    KSAT --> SIMWE
    SIMWE --> DEPTH("Water depth raster")
    SIMWE --> DISCHARGE("Discharge raster")

Objectives

  1. Import SSURGO soil data into GRASS using r.in.ssurgo.
  2. Visualize the imported soil maps, including soil areas, hydrologic groups, ksat maps, and mukey maps.
  3. Calculate curve number using the imported hydrologic group map and a land cover map.
  4. Calculate runoff depth and peak discharge using the curve number map, flow direction, and time of concentration.
  5. Run a basic SIMWE simulation using the imported soil data and visualize the results.

We will examine the Coweeta site for this demonstration, which has a mapset called ssurgo_demo that we will use to import the SSURGO data and run the SIMWE model.

Import SSURGO Data

Install the r.in.ssurgo extension if it is not already installed. This extension is available in the GRASS Addons repository and can be installed using g.extension. The source code for this extension is also available in the tools/r.soildb/r.in.ssurgo directory of this repository.

g.extension r.in.ssurgo
tools.g_extension(
    extension="r.in.ssurgo"
)

Set the computational region to the elevation raster and create a relief raster for hillshading.

Code 
region_text = tools.g_region(raster="elevation", flags="p").text
print(region_text)
projection: 1 (UTM)
zone:       17
datum:      nad83
ellipsoid:  grs80
north:      3884960
south:      3876990
west:       272970
east:       280150
nsres:      10
ewres:      10
rows:       797
cols:       718
cells:      572246

Let’s also create a relief raster for hillshading or soil maps.

Code 
tools.r_relief(input="elevation", output="relief")

Import SSURGO data using r.in.ssurgo. Note that this will import the soil areas map, hydrologic group map, and ksat maps for low, regular, and high values. The mukey column is also imported as a raster with a categorical color scheme applied. The hydrologic group map also has a categorical color scheme applied.

Code 
tools.r_in_ssurgo(
    ssurgo_path="../data/gSSURGO_NC.zip",
    soils="soil_areas",
    hydgrp="hydgrp",
    ksat_l="ksat_l",
    ksat_r="ksat_r",
    ksat_h="ksat_h",
    mukey="mukey",
    hzdept_r=0,
    hzdepb_r=100,
    desgnmaster="A",
    nprocs=30
)
Code 
tools.r_in_ssurgo(
    soils="soil_areas_sda",
    hydgrp="hydgrp_sda",
    ksat_l="ksat_l_sda",
    ksat_r="ksat_r_sda",
    ksat_h="ksat_h_sda",
    mukey="mukey_sda",
    hzdept_r=0,
    hzdepb_r=100,
    desgnmaster="A",
    nprocs=30
)

Visualize Imported Data

The tool imports the soil polygons as a vector map (soil_areas_sda), which contains the soil polygons with attributes from the SSURGO database.

::: {#cell-SDA soil attributes .cell execution_count=8}

cat AREASYMBOL SPATIALVER MUSYM MUKEY Shape_Length Shape_Area
0 1 NC113 5 PwE 545853 562.108989 16443.825
1 2 NC113 5 PwD 545852 781.282178 27038.250
2 3 NC113 5 EvC 545835 1195.591245 28988.135
3 4 NC113 5 DrB 545822 782.617111 20451.195
4 5 NC113 5 PwE 545853 470.492194 9361.455

:::

MUKEY Map

Ksat Maps

Low

Low

Regular

Regular

High

High

Ksat maps for low, regular, and high values

Hydrologic Group

We can visualize the hydrologic group maps to compare the local and SDA imports. The hydrologic group map is a categorical raster with values of A, B, C, D, and A/D. The color scheme is applied using r.colors with the hydgrp color table.

Hydrologic groups

Hydrologic groups

Calculate Curve Number

Install the r.curvenumber extension if it is not already installed. This extension is available in the GRASS Addons repository and can be installed using g.extension.

# | label: install_curvenumber_text
# | eval: false
tools.g_extension(
    extension="r.curvenumber",
)

To calculate the curve number, we need both a land cover map and a hydrologic soil group map. The land cover map is from the National Land Cover Database (NLCD) 2024 dataset, which was imported as a COG raster. The landcover_source parameter specifies that the land cover map is from the NLCD dataset, which allows the tool to use the appropriate lookup table for curve number values.

Code 
tools.r_import(
    output="nlcd_2024",
    input="<path_to_nlcd_cog>",
    extent="region",
    resolution="region",
)

NLCD 2024 land cover map:

Code 
with gs.RegionManager(raster="nlcd_2024", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="nlcd_2024")
    m.d_legend(
        raster="nlcd_2024",
        title="NLCD 2024",
        flags="",
        at="5,80,2,45",
        fontsize=12,
    )
    m.show()

Calculate the hydrolic curve number using the r.curvenumber extension, which requires both a land cover map and a hydrologic soil group map. The land cover map is from the National Land Cover Database (NLCD) 2024 dataset, which was imported as a COG raster. The landcover_source parameter specifies that the land cover map is from the NLCD dataset, which allows the tool to use the appropriate lookup table for curve number values.

Code 
tools.r_curvenumber(
    landcover="nlcd_2024",
    soil="hydgrp_sda",
    landcover_source="nlcd",
    output="curvenumber",
)
Code 
with gs.RegionManager(raster="curvenumber", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="curvenumber")
    m.d_legend(raster="curvenumber", title="Curve Number", flags="s", at="5, 30, 2, 10", fontsize=12)
    m.show()

Calculate Runoff and time of concentration

Let’s also run r.runoff to calculate runoff depth using the curve number method. This tool requires a precipitation raster, which we can create using r.mapcalc with a constant value for precipitation.

Code 
# tools.r_noaa_atlas14(
#     mode="grid",
#     statistic="intensity",
#     unit="metric",
#     series="pds",
#     bound="expected",
#     resample="bilinear",
#     flags="i",
#     # region="se",
#     base_gis_url="https://hdsc.nws.noaa.gov/pub/hdsc/data",
# )
# | Latitude | 35.0667° |
# | Longitude | -83.4333° |
outjson = tools.r_noaa_atlas14(
    mode="point",
    statistic="intensity",
    vector_output="noaa_atlas14_station",
    coordinates=[-83.4333, 35.0667],
    unit="metric",
    series="pds",
    bound="expected",
    resample="bilinear",
    # flags="i",
    # region="se",
    # base_gis_url="https://hdsc.nws.noaa.gov/pub/hdsc/data",
).json

pprint(outjson)
{'bound': 'expected',
 'metadata': {'data_type': 'Precipitation intensity',
              'elevation_(usgs)': 'None None',
              'latitude': '35.0667 Degree',
              'location_name_(esri_maps)': 'None',
              'longitude': '-83.4333 Degree',
              'project_area': 'Ohio River Basin',
              'station_name': 'None',
              'time_series_type': 'Partial duration'},
 'request': {'lat': 35.0667,
             'lon': -83.4333,
             'series': 'pds',
             'statistic': 'intensity',
             'units': 'metric',
             'url': 'https://hdsc.nws.noaa.gov/cgi-bin/new/fe_text.csv?lat=35.0667&lon=-83.4333&data=intensity&units=metric&series=pds'},
 'table': {'return_periods_years': [1, 2, 5, 10, 25, 50, 100, 200, 500, 1000],
           'rows': [{'duration': '5-min',
                     'values': {'1': 123.0,
                                '10': 192.0,
                                '100': 269.0,
                                '1000': 353.0,
                                '2': 144.0,
                                '200': 293.0,
                                '25': 220.0,
                                '5': 166.0,
                                '50': 246.0,
                                '500': 323.0}},
                    {'duration': '10-min',
                     'values': {'1': 98.0,
                                '10': 154.0,
                                '100': 214.0,
                                '1000': 278.0,
                                '2': 115.0,
                                '200': 232.0,
                                '25': 175.0,
                                '5': 133.0,
                                '50': 196.0,
                                '500': 256.0}},
                    {'duration': '15-min',
                     'values': {'1': 82.0,
                                '10': 130.0,
                                '100': 180.0,
                                '1000': 232.0,
                                '2': 96.0,
                                '200': 195.0,
                                '25': 148.0,
                                '5': 112.0,
                                '50': 165.0,
                                '500': 214.0}},
                    {'duration': '30-min',
                     'values': {'1': 56.0,
                                '10': 94.0,
                                '100': 138.0,
                                '1000': 188.0,
                                '2': 66.0,
                                '200': 152.0,
                                '25': 110.0,
                                '5': 80.0,
                                '50': 124.0,
                                '500': 171.0}},
                    {'duration': '60-min',
                     'values': {'1': 35.0,
                                '10': 61.0,
                                '100': 95.0,
                                '1000': 137.0,
                                '2': 42.0,
                                '200': 107.0,
                                '25': 73.0,
                                '5': 51.0,
                                '50': 84.0,
                                '500': 122.0}},
                    {'duration': '2-hr',
                     'values': {'1': 21.0,
                                '10': 36.0,
                                '100': 56.0,
                                '1000': 81.0,
                                '2': 25.0,
                                '200': 63.0,
                                '25': 43.0,
                                '5': 30.0,
                                '50': 49.0,
                                '500': 72.0}},
                    {'duration': '3-hr',
                     'values': {'1': 15.0,
                                '10': 26.0,
                                '100': 41.0,
                                '1000': 61.0,
                                '2': 18.0,
                                '200': 46.0,
                                '25': 31.0,
                                '5': 22.0,
                                '50': 36.0,
                                '500': 54.0}},
                    {'duration': '6-hr',
                     'values': {'1': 10.0,
                                '10': 17.0,
                                '100': 26.0,
                                '1000': 39.0,
                                '2': 12.0,
                                '200': 30.0,
                                '25': 20.0,
                                '5': 14.0,
                                '50': 23.0,
                                '500': 35.0}},
                    {'duration': '12-hr',
                     'values': {'1': 7.0,
                                '10': 11.0,
                                '100': 16.0,
                                '1000': 23.0,
                                '2': 8.0,
                                '200': 18.0,
                                '25': 13.0,
                                '5': 9.0,
                                '50': 14.0,
                                '500': 20.0}},
                    {'duration': '24-hr',
                     'values': {'1': 4.0,
                                '10': 7.0,
                                '100': 11.0,
                                '1000': 15.0,
                                '2': 5.0,
                                '200': 12.0,
                                '25': 8.0,
                                '5': 6.0,
                                '50': 10.0,
                                '500': 13.0}},
                    {'duration': '2-day',
                     'values': {'1': 3.0,
                                '10': 4.0,
                                '100': 6.0,
                                '1000': 8.0,
                                '2': 3.0,
                                '200': 7.0,
                                '25': 5.0,
                                '5': 4.0,
                                '50': 6.0,
                                '500': 8.0}},
                    {'duration': '3-day',
                     'values': {'1': 2.0,
                                '10': 3.0,
                                '100': 4.0,
                                '1000': 6.0,
                                '2': 2.0,
                                '200': 5.0,
                                '25': 3.0,
                                '5': 3.0,
                                '50': 4.0,
                                '500': 5.0}},
                    {'duration': '4-day',
                     'values': {'1': 1.0,
                                '10': 2.0,
                                '100': 3.0,
                                '1000': 4.0,
                                '2': 2.0,
                                '200': 4.0,
                                '25': 3.0,
                                '5': 2.0,
                                '50': 3.0,
                                '500': 4.0}},
                    {'duration': '7-day',
                     'values': {'1': 1.0,
                                '10': 2.0,
                                '100': 2.0,
                                '1000': 3.0,
                                '2': 1.0,
                                '200': 2.0,
                                '25': 2.0,
                                '5': 1.0,
                                '50': 2.0,
                                '500': 3.0}},
                    {'duration': '10-day',
                     'values': {'1': 1.0,
                                '10': 1.0,
                                '100': 2.0,
                                '1000': 2.0,
                                '2': 1.0,
                                '200': 2.0,
                                '25': 1.0,
                                '5': 1.0,
                                '50': 2.0,
                                '500': 2.0}},
                    {'duration': '20-day',
                     'values': {'1': 1.0,
                                '10': 1.0,
                                '100': 1.0,
                                '1000': 1.0,
                                '2': 1.0,
                                '200': 1.0,
                                '25': 1.0,
                                '5': 1.0,
                                '50': 1.0,
                                '500': 1.0}},
                    {'duration': '30-day',
                     'values': {'1': 0.0,
                                '10': 1.0,
                                '100': 1.0,
                                '1000': 1.0,
                                '2': 1.0,
                                '200': 1.0,
                                '25': 1.0,
                                '5': 1.0,
                                '50': 1.0,
                                '500': 1.0}},
                    {'duration': '45-day',
                     'values': {'1': 0.0,
                                '10': 1.0,
                                '100': 1.0,
                                '1000': 1.0,
                                '2': 0.0,
                                '200': 1.0,
                                '25': 1.0,
                                '5': 0.0,
                                '50': 1.0,
                                '500': 1.0}},
                    {'duration': '60-day',
                     'values': {'1': 0.0,
                                '10': 0.0,
                                '100': 1.0,
                                '1000': 1.0,
                                '2': 0.0,
                                '200': 1.0,
                                '25': 1.0,
                                '5': 0.0,
                                '50': 1.0,
                                '500': 1.0}}]}}

We need the flow direction and stream maps to run r.runoff with time of concentration, which can be calculated using the r.watershed tool.

Code 
tools.r_watershed(
    elevation="elevation",
    drainage="flow_direction",
    accumulation="accumulation",
    stream="stream",
    basin="basins",
    threshold=10,
)
Code 
with gs.RegionManager(raster="accumulation", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="accumulation")
    m.d_legend(
        raster="accumulation",
        title="Accumulation",
        flags="s",
        at="5, 30, 2, 10",
        fontsize=12,
    )
    m.show()

    m2 = gj.Map(use_region=True)
    m2.d_shade(shade="relief", color="basins")
    m2.d_legend(
        raster="basins",
        title="Basins 50K threshold",
        flags="s",
        at="5, 30, 2, 10",
        fontsize=12,
    )
    m2.show()

Let calculate the time of concentration using the r.timeofconcentration tool, which requires an elevation raster, flow direction raster, and stream raster as inputs. The output is a raster of time of concentration values in hours.

Code 
tools.r_timeofconcentration(
    elevation="elevation",
    direction="flow_direction",
    stream="stream",
    length_min=100,  # minimum length of flow path
    tc="time_concentration",
)
Code 
with gs.RegionManager(raster="elevation", s="s-2100", flags="a"):
    tools.r_colors(map="time_concentration", color="magma", flags="e")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="time_concentration")

    m.d_legend(
        raster="time_concentration",
        at="8, 12, 20, 80",
        title="Time of Concentration (hours)",
        font="Fira Sans Condensed Light",
        fontsize=14,
        flags="s",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Let’s check the univarite statistics for the time of concentration raster to see the range of values.

Table 1: Time of Concentration — Summary Statistics
Statistic Value
n 57,632.000
null_cells 514,614.000
cells 572,246.000
min 0.011
max 0.783
range 0.772
mean 0.055
mean_of_abs 0.055
stddev 0.077
variance 0.006
coeff_var 141.459
sum 3,147.735
first_quartile 0.022
median 0.031
third_quartile 0.053
percentile_percentile 90.000
percentile_value 0.100

We can use the the time_concentration raster as an input to the r.runoff tool, which calculates runoff depth and peak discharge using the curve number method. The r.runoff tool also requires a precipitation raster, which we can create using r.mapcalc with a constant value for precipitation (\(50\) mm) over a 24 hour period.

Code 
df_tc_min = df_tc[["min", "mean", "max"]].apply(
    lambda x: round(x * 60, 2)
)  # convert to minutes
max_tc = df_tc_min["max"].max()
print(max_tc)

df_rf_intensity = pd.read_csv("data/r_intensity_pds.csv")
display(
    df_rf_intensity.style.hide(axis="index")
    .set_caption("Rainfall Intensity — NOAA Atlas 14, Coweeta, NC")
    .format(
        {
            "Value": lambda v: f"{v:,.3f}"
            if isinstance(v, (int, float, np.number))
            else v
        }
    )
)

# | label: rainfall_intensity_figure
# | echo: false
# | fig-cap: "NOAA Atlas 14 rainfall intensity by storm duration for Coweeta, NC, with the maximum time of concentration marked."

rf = df_rf_intensity.copy()

duration_col = next(
    (c for c in rf.columns if "duration" in c.lower()),
    rf.columns[0],
)

def _duration_to_minutes(value):
    text = str(value).strip().lower()
    match = re.search(r"(\d+(?:\.\d+)?)", text)
    if not match:
        return np.nan
    number = float(match.group(1))
    if "day" in text:
        return number * 24 * 60
    if "hr" in text or "hour" in text:
        return number * 60
    return number

import re

rf["duration_min"] = rf[duration_col].apply(_duration_to_minutes)

value_cols = [c for c in rf.columns if c not in {duration_col, "duration_min"}]
rf_long = rf.melt(
    id_vars=["duration_min"],
    value_vars=value_cols,
    var_name="return_period",
    value_name="intensity_mm_hr",
)
rf_long["intensity_mm_hr"] = pd.to_numeric(rf_long["intensity_mm_hr"], errors="coerce")
rf_long = rf_long.dropna(subset=["duration_min", "intensity_mm_hr"]).sort_values(
    ["return_period", "duration_min"]
)

sns.set_theme(style="whitegrid", context="talk")
fig, ax = plt.subplots(figsize=(10, 6))

sns.lineplot(
    data=rf_long,
    x="duration_min",
    y="intensity_mm_hr",
    hue="return_period",
    marker="o",
    linewidth=2.5,
    ax=ax,
)

ax.axvline(
    max_tc,
    color="black",
    linestyle="--",
    linewidth=1.5,
    label=f"Max Tc = {max_tc:.1f} min",
)

ax.annotate(
    f"Max Tc\n{max_tc:.1f} min",
    xy=(max_tc, rf_long["intensity_mm_hr"].max()),
    xytext=(8, -10),
    textcoords="offset points",
    ha="left",
    va="top",
    fontsize=11,
    bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="0.6", alpha=0.9),
)

ax.set_xscale("log")
ax.set_xlabel("Storm duration (minutes)")
ax.set_ylabel("Rainfall intensity (mm/hr)")
ax.set_title("Rainfall intensity–duration curves")
ax.legend(title="Return period", frameon=True, ncol=2)
sns.despine()
plt.tight_layout()
plt.show()


return_periods = rf_long["return_period"].dropna().unique()
rp_100 = next(
    (
        rp
        for rp in return_periods
        if re.search(r"100", str(rp)) and re.search(r"(yr|year)", str(rp).lower())
    ),
    next((rp for rp in return_periods if re.search(r"100", str(rp))), None),
)

if rp_100 is None:
    raise ValueError("Could not identify the 100-year storm column in rainfall intensity data.")

rf_100 = (
    rf_long.loc[rf_long["return_period"] == rp_100, ["duration_min", "intensity_mm_hr"]]
    .dropna()
    .sort_values("duration_min")
)

storm_intensity_100yr = float(
    np.interp(
        np.log10(max_tc),
        np.log10(rf_100["duration_min"].to_numpy()),
        rf_100["intensity_mm_hr"].to_numpy(),
    )
)

rain_mm_per_hour = storm_intensity_100yr

display(
    pd.DataFrame(
        {
            "return_period": [rp_100],
            "max_tc_min": [round(max_tc, 2)],
            "intensity_mm_hr": [round(storm_intensity_100yr, 2)],
        }
    ).style.hide(axis="index").set_caption(
        "Interpolated rainfall intensity at maximum time of concentration"
    )
)

print(f"100-year storm intensity at max Tc ({max_tc:.2f} min): {storm_intensity_100yr:.2f} in/hr")
46.97
Table 2: Rainfall Intensity — NOAA Atlas 14, Coweeta, NC
duration 1 2 5 10 25 50 100 200 500 1000
5-min 4.850000 5.660000 6.550000 7.570000 8.660000 9.670000 10.600000 11.500000 12.700000 13.900000
10-min 3.880000 4.520000 5.240000 6.050000 6.910000 7.700000 8.420000 9.150000 10.100000 10.900000
15-min 3.230000 3.790000 4.420000 5.100000 5.840000 6.500000 7.100000 7.700000 8.440000 9.150000
30-min 2.210000 2.620000 3.140000 3.700000 4.320000 4.890000 5.430000 5.990000 6.720000 7.410000
60-min 1.380000 1.640000 2.010000 2.410000 2.880000 3.320000 3.740000 4.200000 4.820000 5.410000
2-hr 0.816000 0.969000 1.180000 1.410000 1.680000 1.940000 2.190000 2.460000 2.830000 3.190000
3-hr 0.598000 0.707000 0.855000 1.020000 1.220000 1.420000 1.620000 1.830000 2.130000 2.410000
6-hr 0.395000 0.464000 0.552000 0.653000 0.781000 0.903000 1.030000 1.170000 1.360000 1.550000
12-hr 0.263000 0.308000 0.365000 0.427000 0.502000 0.571000 0.638000 0.710000 0.805000 0.896000
24-hr 0.166000 0.198000 0.244000 0.281000 0.332000 0.374000 0.418000 0.464000 0.528000 0.581000
2-day 0.100000 0.120000 0.146000 0.166000 0.195000 0.219000 0.243000 0.268000 0.303000 0.331000
3-day 0.072000 0.086000 0.104000 0.118000 0.137000 0.153000 0.168000 0.185000 0.207000 0.225000
4-day 0.058000 0.069000 0.083000 0.093000 0.108000 0.119000 0.131000 0.143000 0.159000 0.172000
7-day 0.039000 0.047000 0.056000 0.063000 0.073000 0.081000 0.090000 0.098000 0.109000 0.118000
10-day 0.032000 0.038000 0.045000 0.050000 0.058000 0.064000 0.070000 0.076000 0.084000 0.091000
20-day 0.021000 0.025000 0.029000 0.032000 0.036000 0.039000 0.042000 0.045000 0.049000 0.052000
30-day 0.017000 0.020000 0.023000 0.025000 0.028000 0.030000 0.032000 0.034000 0.036000 0.038000
45-day 0.015000 0.017000 0.019000 0.021000 0.023000 0.024000 0.025000 0.027000 0.028000 0.029000
60-day 0.013000 0.015000 0.017000 0.019000 0.020000 0.021000 0.022000 0.023000 0.024000 0.025000

Table 3: Interpolated rainfall intensity at maximum time of concentration
return_period max_tc_min intensity_mm_hr
100 46.970000 4.340000 
100-year storm intensity at max Tc (46.97 min): 4.34 in/hr
Code 
# storm_intensity_100yr 46.97 min = 0.7828 hr
duration_min = 46.97
intensity_in_hr = 4.34
intensity_mm_hr = intensity_in_hr * 25.4

duration_hr = duration_min / 60.0
rain_mm = intensity_mm_hr * duration_hr

print(f"Storm duration: {duration_hr:.4f} hr")
print(f"Rainfall intensity: {intensity_mm_hr:.3f} mm/hr")
print(f"Rain depth: {rain_mm:.3f} mm")
tools.r_mapcalc(
    expression=f"rain = {rain_mm}",  # mm of precipitation
)
Storm duration: 0.7828 hr
Rainfall intensity: 110.236 mm/hr
Rain depth: 86.296 mm
Code 
tools.r_runoff(
    rainfall="rain",
    direction="flow_direction",
    duration=duration_hr,  # 0.7828 hrs
    curvenumber="curvenumber",
    time_concentration="time_concentration",
    runoff_depth="runoff_depth",
    runoff_volume="runoff_volume",
    upstream_area="upstream_area",
    upstream_runoff_depth="upstream_runoff_depth",
    upstream_runoff_volume="upstream_runoff_volume",
    time_to_peak="time_to_peak",
    peak_discharge="peak_discharge",
)

Let’s visualize the runoff depth

Code 
with gs.RegionManager(raster="runoff_depth", w="w-1600", flags="a"):
    tools.r_colors(map="runoff_depth", color="water")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="runoff_depth")
    m.d_legend(
        raster="runoff_depth", at="5, 30, 2, 10", title="Runoff Depth (mm)", flags="b"
    )
    m.show()

Peak discharge and time to peak discharge maps

Code 
with gs.RegionManager(raster="peak_discharge", w="w-1600", flags="a"):
    tools.r_colors(map="peak_discharge", color="magma", flags="")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="peak_discharge")

    m.d_legend(
        raster="peak_discharge",
        at="5, 30, 2, 10",
        title="Peak Discharge (m³/s)",
        flags="d",
    )
    m.show()

    tools.r_colors(map="time_to_peak", color="magma", flags="e")
    m1 = gj.Map(use_region=True)
    m1.d_shade(shade="relief", color="time_to_peak")

    m1.d_legend(
        raster="time_to_peak",
        at="5, 30, 2, 10",
        title="Time to Peak (hours)",
        flags="d",
    )
    m1.show()

Manning’s n Map

We will calculate the Manning’s n map using the r.manning extension, which uses the National Land Cover Dataset (2024) map to assign Manning’s n values based on the lookup table from Kalyanapu et al. (2009). The output is a raster of Manning’s n values that can be used as an input to the SIMWE model.

We first must install the r.manning extension if it is not already installed. This extension is available in the GRASS Addons repository and can be installed using g.extension.

Code 
tools.g_extension(
    extension="r.manning"
)

Now we can calculate the roughness map using the r.manning tool.

Code 
tools.r_manning(
    input="nlcd_2024",
    output="mannings_n",
    landcover="nlcd",
    source="kalyanapu",
    seed=42,
)


with gs.RegionManager(raster="mannings_n", s="s-2100", flags="a"):
    tools.r_colors(map="mannings_n", color="plasma", flags="")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="mannings_n")

    m.d_legend(
        raster="mannings_n",
        at="8, 12, 20, 80",
        title="Manning's n",
        font="Fira Sans Condensed Light",
        fontsize=14,
        flags="",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Run SIMWE model

Calcuate the slope and aspect from the elevation raster, which are required inputs for the r.sim.water tool.

Precipitation and Rainfall Excess

We will use a constant rainfall value for the simulation, which is based on the 60 min 100 year storm from NOAA Atlas 14 for Coweeta, NC. The rainfall excess is calculated as the runoff depth divided by the duration of the storm (24 hours) to get a rainfall excess in mm/hr.

Location Information Value
Name Otto, North Carolina, USA*
Station name COWEETA EXP STN
Site ID 31-2102
Latitude 35.0667°
Longitude -83.4333°
Elevation 2700 ft**

Coweeta Hydrologic Lab

Coweeta Hydrologic Lab

Rainfall excess map

Code 
# Time of Concentration: 46.8 min
tools.r_mapcalc(
    expression=f"rainfall_excess = runoff_depth / 24.0",
)


with gs.RegionManager(raster="rainfall_excess", s="s-2100", flags="a"):
    tools.r_colors(map="rainfall_excess", color="magma", flags="e")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="rainfall_excess")

    m.d_legend(
        raster="rainfall_excess",
        title="Rainfall Excess (mm/hr)",
        at="8, 12, 20, 80",
        font="Fira Sans Condensed Light",
        fontsize=14,
        flags="",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Calcuate the overland flow for a 47 min 100 year storm using the r.sim.water tool, which simulates overland flow using a particle-based approach. The tool requires the elevation raster, rainfall excess raster, and Manning’s n raster as inputs, along with the number of walkers to simulate and the number of iterations to run. The output is a series of depth and discharge rasters that can be visualized as maps or animations.

Code 
# duration_min = 46.97 round to 47 min for simulation
tools.r_slope_aspect(elevation="elevation", dx="dx", dy="dy")
tools.r_sim_water(
    elevation="elevation",
    depth="tc_100yr_depth",
    discharge="tc_100yr_discharge",
    walkers_output="tc_100yr_walkers",
    rain="rainfall_excess",  # Infiltration / Veg Intercept
    man="mannings_n",  # Mixed forest Kalyanapu et al. (2009)
    nwalkers=100000,
    # infil="ksat_r_sda", # Infiltration of flowing water
    # niterations=1440 / 2,
    observation="sample_outlets",
    logfile="coweeta_100yr_simulation.log",
    duration=48, # Minutes
    # mintimestep=0,
    output_step=2,
    nprocs=30,
    flags="t",
)

# Register the output maps into a space time dataset
gs.run_command(
    "t.create",
    output="tc_100yr_depth_sum",
    type="strds",
    temporaltype="absolute",
    title="Runoff Depth",
    description="Runoff Depth in [m]",
    overwrite=True,
)

# Get the list of depth maps
depth_list = gs.read_command(
    "g.list", type="raster", pattern="tc_100yr_depth.*", separator="comma"
).strip()

# Register the maps
gs.run_command(
    "t.register",
    input="tc_100yr_depth_sum",
    type="raster",
    start="2024-01-01",
    increment="2 minutes",
    maps=depth_list,
    flags="i",
    overwrite=True,
)
Code 
with gs.RegionManager(raster="accumulation", s="s-2100", flags="a"):
    depth_list = gs.read_command(
        "g.list", type="raster", pattern="tc_100yr_depth.*", separator="comma"
    ).strip()

    max_depth = depth_list.split(",")[-1]
    stats = tools.r_info(map=max_depth, format="json").json
    m = gj.Map(use_region=True, width=600)
    m.d_shade(shade="relief", color="naip_2022_rgb@naip")
    m.d_rast(map=max_depth, values=f"0.005-{stats['max']}")
    m.d_vect(
        map="sample_outlets",
        type="point",
        color="black",
        fill_color="orange",
        icon="basic/circle",
        size=10,
    )
    m.d_legend(
        raster=max_depth,
        title="Deqpth [m]",
        at="8, 12, 15, 85",
        font="Fira Sans Condensed Light",
        range=[0.005, stats["max"]],
        fontsize=14,
        flags="",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Animation of depth maps

Code 
depth_list = gs.read_command(
    "g.list", type="raster", pattern="tc_100yr_depth.*", separator="comma"
).strip()

max_depth = depth_list.split(",")[-1]
stats = tools.r_info(map=max_depth, format="json").json
m = gj.TimeSeriesMap(use_region=True)
m.d_rast(map="relief")
m.add_raster_series("tc_100yr_depth_sum", values=f"0.025-{stats['max']}")
m.d_legend(
    raster=max_depth,
    title="Depth [m]",
    flags="bt",
    border_color="none",
    range=[0.005, stats["max"]],
    at=(1, 50, 0, 5),
)

m.save("output/tc_100yr_simwe_depth_animation.gif")
m.show()

Overland Flow Simulation with Soil Data

Overland Flow Simulation with Soil Data

Overland Flow Simulation (100yr - 12hr storm) with Soil Data

Code 
m = gj.TimeSeriesMap(use_region=True)
m.d_rast(map="relief")
m.add_raster_series("depth_sum", values="0.025-1.25")
m.d_legend(
    raster="depth.720",
    title="Depth [m]",
    flags="bt",
    border_color="none",
    at=(1, 50, 0, 5),
)
# m.d_shade(shade="relief", color="depth")
m.save("output/simwe_depth_animation.gif")
m.show()

Overland Flow Simulation (100yr - 12hr storm) with Soil Data

Overland Flow Simulation (100yr - 12hr storm) with Soil Data

No soils

Code 
# 12 hr 100 year storm is 16.2052 mm/hr
tools.r_sim_water(
    elevation="elevation",
    depth="no_soil_depth",
    discharge="no_soil_discharge",
    walkers_output="no_soil_walkers",
    rain_value=rain_mm_per_hour,
    # rain="rainfall_excess",  # Infiltration / Veg Intercept
    man="mannings_n",  # Mixed forest Kalyanapu et al. (2009)
    nwalkers=100000,
    # infil="ksat_r_sda", # Infiltration of flowing water
    niterations=1440 / 2,
    # mintimestep=0,
    output_step=2,
    nprocs=30,
    flags="t",
)

gs.run_command(
    "t.create",
    output="no_soil_depth_sum",
    type="strds",
    temporaltype="absolute",
    title="Runoff Depth",
    description="Runoff Depth in [m]",
    overwrite=True,
)

# Get the list of depth maps
depth_list = gs.read_command(
    "g.list", type="raster", pattern="no_soil_depth.*", separator="comma"
).strip()

# Register the maps
gs.run_command(
    "t.register",
    input="no_soil_depth_sum",
    type="raster",
    start="2024-01-01",
    increment="2 minutes",
    maps=depth_list,
    flags="i",
    overwrite=True,
)

with gs.RegionManager(raster="accumulation", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="naip_2022_rgb@naip")
    m.d_rast(map="no_soil_depth.720", values="0.025-1.25")
    m.d_legend(
        raster="no_soil_depth.720",
        title="Depth [m]",
        flags="btl",
        border_color="none",
        at=(5, 50, 2, 5),
    )
    m.show()

no soil animation

Code 
m = gj.TimeSeriesMap(use_region=True)
m.d_rast(map="relief")
m.add_raster_series("no_soil_depth_sum", values="0.025-1.25")
m.d_legend(
    raster="no_soil_depth.720",
    title="Depth [m]",
    flags="bt",
    border_color="none",
    at=(1, 50, 0, 5),
)
m.save("output/no_soil_simwe_depth_animation.gif")
m.show()

Overland Flow Simulation without Soil Data

Overland Flow Simulation without Soil Data