r.curvenumber: Add dual HSG support and null-soil fallback (#1636)

  Resolve dual hydrologic soil groups (A/D, B/D, C/D, D/D encoded as
  11-14) to single groups based on drainage condition. Add -d flag to
  select drained behavior (uses first group); default is undrained
  (uses second group). Fall back to HSG D when soil is null but
  landcover is valid. Add tests for both behaviors.

Author: petrasovaa

src/raster/r.curvenumber/global.h

@@ -57,4 +57,9 @@ int lookup_lut_cn_ii(const struct cn_table *tbl, int lc, int hsg, int hc_idx,
    returns cn_ii unchanged if not found in table */
 int convert_arc_from_table(const struct arc_table *tbl, int cn_ii, int arc_idx);
 
+/* resolve dual HSG codes (11-14) to single codes (1-4)
+   drained=1: use first group (A/D->A), drained=0: use second group (A/D->D)
+   returns -1 for invalid HSG values */
+int resolve_dual_hsg(int hsg, int drained);
+
 #endif

src/raster/r.curvenumber/lut.c

@@ -359,6 +359,28 @@ int load_builtin_arc(struct arc_table *tbl)
     return 0;
 }
 
+/* resolve dual HSG codes to single codes */
+int resolve_dual_hsg(int hsg, int drained)
+{
+    switch (hsg) {
+    case 1:
+    case 2:
+    case 3:
+    case 4:
+        return hsg;
+    case 11: /* A/D */
+        return drained ? 1 : 4;
+    case 12: /* B/D */
+        return drained ? 2 : 4;
+    case 13: /* C/D */
+        return drained ? 3 : 4;
+    case 14: /* D/D */
+        return 4;
+    default:
+        return -1;
+    }
+}
+
 /* lookups */
 int lookup_lut_cn_ii(const struct cn_table *tbl, int lc, int hsg, int hc_idx,
                      int *cn_ii_out)

src/raster/r.curvenumber/main.c

@@ -29,6 +29,7 @@ int main(int argc, char **argv)
     struct GModule *module;
     struct Option *opt_land, *opt_soil, *opt_source, *opt_lookup;
     struct Option *opt_hc, *opt_arc, *opt_out;
+    struct Flag *flag_drained;
     struct Cell_head region;
 
     const char *land_name, *soil_name, *out_name;
@@ -44,7 +45,8 @@ int main(int argc, char **argv)
 
     struct cn_table cn_tbl;
     struct arc_table arc_tbl;
-    int hc_idx, arc_idx;
+    int hc_idx, arc_idx, drained;
+    int null_soil_warned = 0;
 
     G_gisinit(argv[0]);
 
@@ -100,6 +102,11 @@ int main(int argc, char **argv)
     opt_arc->description = _("Antecedent Runoff Condition (the degree of "
                              "wetness of a watershed before a rainfall event)");
 
+    flag_drained = G_define_flag();
+    flag_drained->key = 'd';
+    flag_drained->description =
+        _("Use drained condition for dual hydrologic soil groups (e.g. A/D)");
+
     opt_out = G_define_standard_option(G_OPT_R_OUTPUT);
     opt_out->description = _("Name for output curve number raster");
 
@@ -116,6 +123,7 @@ int main(int argc, char **argv)
     lookup = opt_lookup->answer;
     hc_str = opt_hc->answer ? opt_hc->answer : "fair";
     arc_str = opt_arc->answer ? opt_arc->answer : "ii";
+    drained = flag_drained->answer;
 
     /* map hydrologic condition */
     if (G_strcasecmp(hc_str, "poor") == 0)
@@ -199,11 +207,31 @@ int main(int argc, char **argv)
             int cn_ii, cn;
             int ok;
 
-            if (Rast_is_c_null_value(&lc_val) ||
-                Rast_is_c_null_value(&soil_val)) {
+            if (Rast_is_c_null_value(&lc_val)) {
+                Rast_set_c_null_value(&out_row[col], 1);
+                continue;
+            }
+
+            /* null soil: fall back to HSG D (most conservative) */
+            if (Rast_is_c_null_value(&soil_val)) {
+                if (!null_soil_warned) {
+                    G_verbose_message(
+                        _("Null soil value(s) found, falling back to HSG D"));
+                    null_soil_warned = 1;
+                }
+                soil_val = 4;
+            }
+
+            /* resolve dual HSG codes (11-14) to single codes (1-4) */
+            int resolved = resolve_dual_hsg((int)soil_val, drained);
+
+            if (resolved < 0) {
+                G_warning(_("Invalid HSG value %d at row %d col %d"),
+                          (int)soil_val, row, col);
                 Rast_set_c_null_value(&out_row[col], 1);
                 continue;
             }
+            soil_val = (CELL)resolved;
 
             ok = lookup_lut_cn_ii(&cn_tbl, (int)lc_val, (int)soil_val, hc_idx,
                                   &cn_ii);

src/raster/r.curvenumber/r.curvenumber.html

@@ -20,6 +20,13 @@ <h3>Hydrologic Soil Groups (HSG)</h3>
   <li>11–14 – Dual groups: A/D, B/D, C/D, D/D (for drained/undrained soils)</li>
 </ul>
 <p>Curve Number increases from HSG A to D due to reduced infiltration associated with increasing clay content and compaction.</p>
+<p>Dual hydrologic soil groups (A/D, B/D, C/D) represent soils that behave
+differently depending on drainage. By default, the undrained (second) group
+is used (e.g., A/D &rarr; D). Use the <b>-d</b> flag to use the drained
+(first) group instead (e.g., A/D &rarr; A).</p>
+<p>When the soil raster contains null values but a valid landcover value exists,
+the module falls back to HSG D (the most conservative group with highest
+runoff potential).</p>
 <p>If you are using SSURGO soil data, Hydrologic Soil Group (HSG) values are typically provided directly. For other soil datasets that do not include HSG classifications, see Reference 1 (Chapter 7) for guidance on assigning HSGs based on soil texture, hydraulic conductivity, and water table depth. This information can be used to develop a custom HSG raster.</p>
 
 <h3>Hydrologic Condition (HC)</h3>

src/raster/r.curvenumber/r.curvenumber.md

@@ -31,6 +31,15 @@ Curve Number increases from HSG A to D
 due to reduced infiltration associated
 with increasing clay content and compaction.
 
+Dual hydrologic soil groups (A/D, B/D, C/D) represent soils that behave
+differently depending on drainage. By default, the undrained (second) group
+is used (e.g., A/D → D). Use the **-d** flag to use the drained (first)
+group instead (e.g., A/D → A).
+
+When the soil raster contains null values but a valid landcover value exists,
+the module falls back to HSG D (the most conservative group with highest
+runoff potential).
+
 If you are using SSURGO soil data,
 Hydrologic Soil Group (HSG) values are
 typically provided directly. For other

src/raster/r.curvenumber/testsuite/test_r_curvenumber.py

@@ -81,5 +81,199 @@ def test_esa_against_expected(self):
         )
 
 
+class TestDualHSG(TestCase):
+    """Test dual HSG resolution and null-soil fallback."""
+
+    lc_single = "lc_dual_test"
+    hsg_dual = "hsg_dual_test"
+    hsg_single_a = "hsg_single_a_test"
+    hsg_single_d = "hsg_single_d_test"
+    lc_null = "lc_null_test"
+    hsg_null = "hsg_null_test"
+    hsg_d_only = "hsg_d_only_test"
+    hsg_invalid = "hsg_invalid_test"
+    out_drained = "cn_drained"
+    out_undrained = "cn_undrained"
+    out_expected_a = "cn_expected_a"
+    out_expected_d = "cn_expected_d"
+    out_null_soil = "cn_null_soil"
+    out_null_expected = "cn_null_expected"
+    out_null_lc = "cn_null_lc"
+    out_invalid_hsg = "cn_invalid_hsg"
+    diff = "cn_dual_diff"
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule("g.region", n=3, s=0, e=3, w=0, rows=1, cols=1, res=3)
+
+        # landcover = ESA 30 (grassland), single pixel
+        cls.runModule("r.mapcalc", expression=f"{cls.lc_single} = 30")
+        # dual HSG A/D = 11
+        cls.runModule("r.mapcalc", expression=f"{cls.hsg_dual} = 11")
+        # single HSG A = 1 (expected result for drained)
+        cls.runModule("r.mapcalc", expression=f"{cls.hsg_single_a} = 1")
+        # single HSG D = 4 (expected result for undrained)
+        cls.runModule("r.mapcalc", expression=f"{cls.hsg_single_d} = 4")
+        # null landcover raster
+        cls.runModule("r.mapcalc", expression=f"{cls.lc_null} = null()")
+        # null soil raster
+        cls.runModule("r.mapcalc", expression=f"{cls.hsg_null} = null()")
+        # HSG D only (expected result for null soil fallback)
+        cls.runModule("r.mapcalc", expression=f"{cls.hsg_d_only} = 4")
+        # invalid HSG value
+        cls.runModule("r.mapcalc", expression=f"{cls.hsg_invalid} = 99")
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        cls.runModule(
+            "g.remove",
+            flags="f",
+            type="raster",
+            name=(
+                cls.lc_single,
+                cls.lc_null,
+                cls.hsg_dual,
+                cls.hsg_single_a,
+                cls.hsg_single_d,
+                cls.hsg_null,
+                cls.hsg_d_only,
+                cls.hsg_invalid,
+            ),
+        )
+
+    def tearDown(self):
+        self.runModule(
+            "g.remove",
+            flags="f",
+            type="raster",
+            name=(
+                self.out_drained,
+                self.out_undrained,
+                self.out_expected_a,
+                self.out_expected_d,
+                self.out_null_soil,
+                self.out_null_expected,
+                self.out_null_lc,
+                self.out_invalid_hsg,
+                self.diff,
+            ),
+        )
+
+    def test_dual_hsg_drained(self):
+        """Dual HSG A/D with -d flag should match single HSG A."""
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_dual,
+            landcover_source="esa",
+            flags="d",
+            output=self.out_drained,
+        )
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_single_a,
+            landcover_source="esa",
+            output=self.out_expected_a,
+        )
+        self.runModule(
+            "r.mapcalc",
+            expression=f"{self.diff} = {self.out_drained} - {self.out_expected_a}",
+        )
+        stats = gs.parse_command("r.univar", flags="g", map=self.diff)
+        self.assertAlmostEqual(
+            float(stats["min"]), 0.0, msg="Drained dual HSG A/D != single HSG A"
+        )
+        self.assertAlmostEqual(
+            float(stats["max"]), 0.0, msg="Drained dual HSG A/D != single HSG A"
+        )
+
+    def test_dual_hsg_undrained(self):
+        """Dual HSG A/D without -d flag should match single HSG D."""
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_dual,
+            landcover_source="esa",
+            output=self.out_undrained,
+        )
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_single_d,
+            landcover_source="esa",
+            output=self.out_expected_d,
+        )
+        self.runModule(
+            "r.mapcalc",
+            expression=f"{self.diff} = {self.out_undrained} - {self.out_expected_d}",
+        )
+        stats = gs.parse_command("r.univar", flags="g", map=self.diff)
+        self.assertAlmostEqual(
+            float(stats["min"]), 0.0, msg="Undrained dual HSG A/D != single HSG D"
+        )
+        self.assertAlmostEqual(
+            float(stats["max"]), 0.0, msg="Undrained dual HSG A/D != single HSG D"
+        )
+
+    def test_invalid_hsg_produces_null(self):
+        """Invalid HSG value should produce null output."""
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_invalid,
+            landcover_source="esa",
+            output=self.out_invalid_hsg,
+        )
+        stats = gs.parse_command("r.univar", flags="g", map=self.out_invalid_hsg)
+        self.assertEqual(
+            int(stats["n"]), 0, msg="Invalid HSG should produce all null cells"
+        )
+
+    def test_null_landcover_produces_null(self):
+        """Null landcover should produce null output."""
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_null,
+            soil=self.hsg_single_d,
+            landcover_source="esa",
+            output=self.out_null_lc,
+        )
+        stats = gs.parse_command("r.univar", flags="g", map=self.out_null_lc)
+        self.assertEqual(
+            int(stats["n"]), 0, msg="Null landcover should produce all null cells"
+        )
+
+    def test_null_soil_fallback(self):
+        """Null soil with valid landcover should fall back to HSG D."""
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_null,
+            landcover_source="esa",
+            output=self.out_null_soil,
+        )
+        self.assertModule(
+            "r.curvenumber",
+            landcover=self.lc_single,
+            soil=self.hsg_d_only,
+            landcover_source="esa",
+            output=self.out_null_expected,
+        )
+        self.runModule(
+            "r.mapcalc",
+            expression=f"{self.diff} = {self.out_null_soil} - {self.out_null_expected}",
+        )
+        stats = gs.parse_command("r.univar", flags="g", map=self.diff)
+        self.assertAlmostEqual(
+            float(stats["min"]), 0.0, msg="Null soil fallback != HSG D"
+        )
+        self.assertAlmostEqual(
+            float(stats["max"]), 0.0, msg="Null soil fallback != HSG D"
+        )
+
+
 if __name__ == "__main__":
     test()