r.soils.rosetta: new addon

src/raster/r.in.ssurgo/r.in.ssurgo.html

@@ -43,6 +43,13 @@ <h2>NOTES</h2>
 rainfall-excess models such as the SCS Curve Number method (see
 <em>r.curvenumber</em>).
 
+<p>
+The soil texture separates (<b>sand</b>, <b>silt</b>, <b>clay</b>, in percent)
+and the bulk density at 1/3 bar (<b>bulk_density</b>, g/cm3) can also be
+written as rasters, aggregated the same way as Ksat. These match the inputs of
+the ROSETTA pedotransfer model, so they can be passed directly to
+<em>r.soils.rosetta</em> to estimate van Genuchten soil hydraulic parameters.
+
 <p>
 The <b>soils</b> output is the source SSURGO Map Unit polygon vector with all
 aggregated attributes attached. Use it for reference or to derive custom
@@ -94,8 +101,9 @@ <h3>Aggregation methods</h3>
 <h3>r3 (3D raster) output for depth profiles</h3>
 
 When <b>depths</b> is set (a comma-separated, strictly-increasing list of
-cm boundaries with at least 2 values), the depth-weighted Ksat outputs
-(<b>ksat_l</b>, <b>ksat_r</b>, <b>ksat_h</b>) are produced as <b>3D
+cm boundaries with at least 2 values), the depth-weighted outputs
+(<b>ksat_l</b>, <b>ksat_r</b>, <b>ksat_h</b>, <b>sand</b>, <b>silt</b>,
+<b>clay</b>, <b>bulk_density</b>) are produced as <b>3D
 rasters</b> with one z-slice per depth bin instead of a single 2D
 average. The number of slices is <code>len(depths) - 1</code>.
 <b>hzdept_r</b> and <b>hzdepb_r</b> are ignored when <b>depths</b> is
@@ -197,6 +205,26 @@ <h3>Import from a local SSURGO geodatabase</h3>
 )
 </pre></div>
 
+<h3>Import texture and estimate hydraulic parameters</h3>
+
+Import soil texture and bulk density for the top 25.4 cm and estimate van
+Genuchten hydraulic parameters with <em>r.soils.rosetta</em> (see its manual
+for the full SSURGO-to-SIMWE workflow):
+
+<div class="code"><pre>
+r.in.ssurgo \
+    ssurgo_path="gSSURGO_NC.zip" \
+    soils="soils" \
+    sand="sand" \
+    silt="silt" \
+    clay="clay" \
+    bulk_density="bd" \
+    hzdept_r=0 \
+    hzdepb_r=25.4
+r.soils.rosetta sand="sand" silt="silt" clay="clay" bulk_density="bd" \
+    ksat="ksat_vg" version=3
+</pre></div>
+
 <h3>Force the SQLite/OGR backend (skip DuckDB)</h3>
 
 <div class="code"><pre>
@@ -254,7 +282,8 @@ <h2>SEE ALSO</h2>
 <em>
     <a href="https://grass.osgeo.org/grass-stable/manuals/addons/r.curvenumber.html">r.curvenumber</a>,
     <a href="https://grass.osgeo.org/grass-stable/manuals/addons/r.runoff.html">r.runoff</a>,
-    <a href="https://grass.osgeo.org/grass-stable/manuals/r.sim.water.html">r.sim.water</a>
+    <a href="https://grass.osgeo.org/grass-stable/manuals/r.sim.water.html">r.sim.water</a>,
+    <a href="https://grass.osgeo.org/grass-stable/manuals/addons/r.soils.rosetta.html">r.soils.rosetta</a>
 </em>
 
 <h2>AUTHORS</h2>

src/raster/r.in.ssurgo/r.in.ssurgo.md

@@ -45,6 +45,12 @@ permeability. The HGS raster (**hydgrp**) are used in rainfall excess models
 such as the SCS Curve Number method to estimate runoff from rainfall
 events (see *r.curvenumber*).
 
+The soil texture separates (**sand**, **silt**, **clay**, in percent) and the
+bulk density at 1/3 bar (**bulk_density**, g/cm3) can also be written as
+rasters, aggregated the same way as Ksat. These match the inputs of the ROSETTA
+pedotransfer model, so they can be passed directly to *r.soils.rosetta* to
+estimate van Genuchten soil hydraulic parameters.
+
 The **soils** output is a vector layer containing the source SSURGO Map Unit polygons
 and attribute data. This can be used for reference or to create custom rasters for
 Curve Number or other applications based on the SSURGO attributes.
@@ -99,8 +105,9 @@ master horizon.
 ### r3 (3D raster) output for depth profiles
 
 When `depths` is set (a comma-separated, strictly-increasing list of cm
-boundaries with at least 2 values), the depth-weighted Ksat outputs
-(**ksat_l**, **ksat_r**, **ksat_h**) are produced as **3D rasters** with
+boundaries with at least 2 values), the depth-weighted outputs
+(**ksat_l**, **ksat_r**, **ksat_h**, **sand**, **silt**, **clay**,
+**bulk_density**) are produced as **3D rasters** with
 one z-slice per depth bin instead of a single 2D average. The number of
 slices is `len(depths) - 1`. `hzdept_r` and `hzdepb_r` are ignored when
 `depths` is set; the aggregation runs once per slice.
@@ -226,6 +233,24 @@ tools.r_in_ssurgo(
 )
 ```
 
+Import soil texture and bulk density for the top 25.4 cm and estimate van
+Genuchten hydraulic parameters with *r.soils.rosetta* (see its manual for the
+full SSURGO-to-SIMWE workflow):
+
+```sh
+r.in.ssurgo \
+    ssurgo_path="gSSURGO_NC.zip" \
+    soils="soils" \
+    sand="sand" \
+    silt="silt" \
+    clay="clay" \
+    bulk_density="bd" \
+    hzdept_r=0 \
+    hzdepb_r=25.4
+r.soils.rosetta sand="sand" silt="silt" clay="clay" bulk_density="bd" \
+    ksat="ksat_vg" version=3
+```
+
 Force the SQLite/OGR backend (skip DuckDB) for a local import:
 
 ```sh
@@ -317,9 +342,10 @@ Accessed [04/23/2025].
 
 ## SEE ALSO
 
-[r.curvenumber](https://grass.osgeo.org/grass-stable/manuals/addons/r.curvenumber.html),
+[r.curvenumber](r.curvenumber.md),
 [r.runoff](https://grass.osgeo.org/grass-stable/manuals/addons/r.runoff.html),
-[r.sim.water](https://grass.osgeo.org/grass-stable/manuals/addons/r.sim.water.html)
+[r.sim.water](https://grass.osgeo.org/grass-stable/manuals/r.sim.water.html),
+[r.soils.rosetta](https://grass.osgeo.org/grass-stable/manuals/addons/r.soils.rosetta.html)
 
 ## AUTHORS
 

src/raster/r.in.ssurgo/r.in.ssurgo.py

@@ -70,6 +70,34 @@
 # % required: no
 # %end
 
+# %option G_OPT_R_OUTPUT
+# % key: sand
+# % description: Sand content raster map [percent]
+# % guisection: Outputs
+# % required: no
+# %end
+
+# %option G_OPT_R_OUTPUT
+# % key: silt
+# % description: Silt content raster map [percent]
+# % guisection: Outputs
+# % required: no
+# %end
+
+# %option G_OPT_R_OUTPUT
+# % key: clay
+# % description: Clay content raster map [percent]
+# % guisection: Outputs
+# % required: no
+# %end
+
+# %option G_OPT_R_OUTPUT
+# % key: bulk_density
+# % description: Bulk density (one-third bar) raster map [g/cm3]
+# % guisection: Outputs
+# % required: no
+# %end
+
 # %option
 # % key: desgnmaster
 # % type: string
@@ -1201,7 +1229,17 @@ def _run_ogr2ogr(args, label):
 
 
 def _rasterize_and_style(
-    ssurgo_vector, hydgrp, ksat_h, ksat_r, ksat_l, mukey, slices=None
+    ssurgo_vector,
+    hydgrp,
+    ksat_h,
+    ksat_r,
+    ksat_l,
+    mukey,
+    sand="",
+    silt="",
+    clay="",
+    bulk_density="",
+    slices=None,
 ):
     """Convert imported SSURGO vector attributes to raster maps and apply color schemes.
 
@@ -1211,8 +1249,13 @@ def _rasterize_and_style(
     :param str ksat_r: Output name for Ksat regular raster (or empty to skip).
     :param str ksat_l: Output name for Ksat low raster (or empty to skip).
     :param str mukey: Output name for map unit key raster (or empty to skip).
-    :param slices: Optional list of ``(top, bottom)`` cm slices. When set,
-        ``ksat_l/r/h`` are produced as r3 (3D) rasters with one slice per
+    :param str sand: Output name for sand percent raster (or empty to skip).
+    :param str silt: Output name for silt percent raster (or empty to skip).
+    :param str clay: Output name for clay percent raster (or empty to skip).
+    :param str bulk_density: Output name for bulk density raster (or empty to skip).
+    :param slices: Optional list of ``(top, bottom)`` cm slices. When set, the
+        depth-weighted rasters (``ksat_l/r/h``, ``sand``, ``silt``, ``clay``,
+        ``bulk_density``) are produced as r3 (3D) rasters with one slice per
         depth bin; ``hydgrp`` and ``mukey`` are still produced as 2D rasters
         because they're profile-level / identity values.
     """
@@ -1223,19 +1266,31 @@ def _rasterize_and_style(
         ("hsg", hydgrp, "hydgrp"),
         ("mukey_int", mukey, None),
     ]
+    # Depth-weighted attributes are 2D in the single-interval (non-sliced) path
+    # and 3D when depth slices are requested. Each entry maps a chorizon column
+    # to its requested output name.
+    depth_weighted = [
+        ("ksat_h", ksat_h),
+        ("ksat_r", ksat_r),
+        ("ksat_l", ksat_l),
+        ("sandtotal_r", sand),
+        ("silttotal_r", silt),
+        ("claytotal_r", clay),
+        ("dbthirdbar_r", bulk_density),
+    ]
     if not slices:
-        # No slicing — depth-weighted Ksat rasters are 2D too (legacy path).
-        _2d_maps.extend(
-            [
-                ("ksat_h", ksat_h, None),
-                ("ksat_r", ksat_r, None),
-                ("ksat_l", ksat_l, None),
-            ]
-        )
+        _2d_maps.extend((col, name, None) for col, name in depth_weighted)
 
     for col, map_name, label_column in _2d_maps:
         if not map_name:
             continue
+        # SDA results arrive as schemaless GeoJSON: a column that is NULL for
+        # every feature imports as TEXT and v.to.rast rejects it, so recast
+        # to double precision before rasterizing (no-op for numeric columns).
+        _ensure_numeric_column(ssurgo_vector, col)
+        # v.to.rast replaces NULL attributes with 0, which would fake
+        # meaningful values (e.g. Ksat 0 mm/hr, sand 0 percent) for map units
+        # without data; filter them out so their cells stay NULL.
         gs.run_command(
             "v.to.rast",
             input=ssurgo_vector,
@@ -1244,6 +1299,7 @@ def _rasterize_and_style(
             attribute_column=col,
             output=map_name,
             label_column=label_column if label_column else "",
+            where=f"{col} IS NOT NULL",
         )
 
         if col == "mukey_int":
@@ -1253,10 +1309,19 @@ def _rasterize_and_style(
             hydrologic_soil_group_color_scheme(map_name)
 
     if slices:
-        for base, name in (("ksat_l", ksat_l), ("ksat_r", ksat_r), ("ksat_h", ksat_h)):
+        for base, name in depth_weighted:
             if not name:
                 continue
-            _rasterize_3d(ssurgo_vector, base_field=base, output=name, slices=slices)
+            # Ksat gets its dedicated 3D color ramp; texture and bulk density
+            # keep the default ramp.
+            color_func = ksat_color_scheme_3d if base.startswith("ksat_") else None
+            _rasterize_3d(
+                ssurgo_vector,
+                base_field=base,
+                output=name,
+                slices=slices,
+                color_func=color_func,
+            )
 
     else:
         ksat_map_names = [name for name in (ksat_l, ksat_r, ksat_h) if name]
@@ -1293,14 +1358,24 @@ def _ensure_numeric_column(vector_map, col):
     gs.run_command("v.db.renamecolumn", map=vector_map, column=f"{tmp_col},{col}")
 
 
-def _rasterize_3d(ssurgo_vector, *, base_field: str, output: str, slices):
+def _rasterize_3d(
+    ssurgo_vector,
+    *,
+    base_field: str,
+    output: str,
+    slices,
+    color_func=ksat_color_scheme_3d,
+):
     """Build a 3D raster ``output`` from per-slice attribute columns.
 
     For each slice ``i`` in ``slices``, rasterizes
     ``<ssurgo_vector>.<base_field>__sN`` to a temporary 2D raster, then stacks
     all of them into a 3D raster with ``r.to.rast3``. The 3D region is set to
     ``b=0 t=max_depth`` so the z-axis represents depth from the surface.
 
+    ``color_func`` is applied to the resulting 3D raster (defaults to the Ksat
+    ramp); pass ``None`` to leave the default color table.
+
     Temporary 2D rasters are deleted on completion.
     """
     gs.message(_("Building 3D raster %s from %d slices...") % (output, len(slices)))
@@ -1320,6 +1395,8 @@ def _rasterize_3d(ssurgo_vector, *, base_field: str, output: str, slices):
             col = f"{base_field}{_slice_suffix(i)}"
             _ensure_numeric_column(ssurgo_vector, col)
             tmp = f"_tmp_{output}_s{i}"
+            # See _rasterize_and_style: keep NULL-attribute map units NULL
+            # instead of v.to.rast's 0 replacement.
             gs.run_command(
                 "v.to.rast",
                 input=ssurgo_vector,
@@ -1328,6 +1405,7 @@ def _rasterize_3d(ssurgo_vector, *, base_field: str, output: str, slices):
                 attribute_column=col,
                 output=tmp,
                 overwrite=True,
+                where=f"{col} IS NOT NULL",
             )
             temp_2d.append(tmp)
         gs.run_command(
@@ -1336,7 +1414,8 @@ def _rasterize_3d(ssurgo_vector, *, base_field: str, output: str, slices):
             output=output,
             overwrite=True,
         )
-        ksat_color_scheme_3d(output)
+        if color_func:
+            color_func(output)
     finally:
         if temp_2d:
             gs.run_command(
@@ -1878,6 +1957,10 @@ def main():
     ksat_r = options["ksat_r"]
     ksat_l = options["ksat_l"]
     mukey = options["mukey"]
+    sand = options["sand"]
+    silt = options["silt"]
+    clay = options["clay"]
+    bulk_density = options["bulk_density"]
 
     # Vector outputs
     ssurgo_vector = options["soils"]
@@ -1995,7 +2078,17 @@ def main():
 
     if ssurgo_vector:
         _rasterize_and_style(
-            ssurgo_vector, hydgrp, ksat_h, ksat_r, ksat_l, mukey, slices=slices
+            ssurgo_vector,
+            hydgrp,
+            ksat_h,
+            ksat_r,
+            ksat_l,
+            mukey,
+            sand=sand,
+            silt=silt,
+            clay=clay,
+            bulk_density=bulk_density,
+            slices=slices,
         )
 
 

src/raster/r.soils.rosetta/CMakeLists.txt

@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 3.22)
+project(r.soils.rosetta)
+
+include(build_addon)
+
+build_addon(
+    ${PROJECT_NAME}
+    SOURCES
+        r.soils.rosetta.py
+    DOCFILES
+        r.soils.rosetta.html
+        r.soils.rosetta.md)

src/raster/r.soils.rosetta/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.soils.rosetta
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script