compute dominant_tillage in import workflow

refs dailyerosion/depmap#13

scripts/import/README.md

@@ -9,6 +9,7 @@ is one file per HUC12.
 1. go to ../util and run `python make_dirs.py <scenario>`
 1. cd to ../import and run `python flowpath2prj.py <scenario>`
 1. `python prj2wepp.py <scenario>`
+1. `python assign_dominant_tillage.py <scenario>`
 1. `python dbset_ofe.py <scenario>`
 1. `python package_myhucs.py <scenario>`
 1. `python check_huc12_zero_flowpaths.py <scenario>`

scripts/import/assign_dominant_tillage.py

@@ -0,0 +1,40 @@
+"""Compute the dominant tillage code for each HUC12."""
+# stdlib
+import sys
+
+# third party
+from pyiem.util import get_dbconn
+
+
+def main(argv):
+    """Do great things."""
+    scenario = int(argv[1])
+    pgconn = get_dbconn("idep")
+    cursor = pgconn.cursor()
+    # Very small percentage have ties, so take max tillage code in that case
+    cursor.execute(
+        """
+        with data as (
+            select huc_12, regexp_split_to_table(management, '') as ch
+            from flowpaths p JOIN flowpath_points t on (p.fid = t.flowpath)
+            where p.scenario = %s
+            ORDER by management desc, flowpath, fpath),
+        agg as (select ch, huc_12, count(*) from data
+            WHERE ch != '0' GROUP by ch, huc_12),
+        agg2 as (select huc_12, ch::int, rank() OVER (PARTITION by huc_12
+            ORDER by count DESC, ch DESC) from agg)
+
+        UPDATE huc12 h SET dominant_tillage = a.ch FROM agg2 a
+        WHERE h.scenario = %s and h.huc_12 = a.huc_12 and
+        a.rank = 1 and
+        (h.dominant_tillage is null or h.dominant_tillage != a.ch)
+        """,
+        (scenario, scenario),
+    )
+    print(f"Updated dominant_tillage for {cursor.rowcount} HUC12s")
+    pgconn.commit()
+    pgconn.close()
+
+
+if __name__ == "__main__":
+    main(sys.argv)

scripts/plots/plot_huc12.py

@@ -14,21 +14,9 @@ def main(argv):
     with get_sqlalchemy_conn("idep") as conn:
         df = gpd.read_postgis(
             """
-            with data as (
-                select huc_12, regexp_split_to_table(management, '') as ch
-                from flowpaths p JOIN flowpath_points t on (p.fid = t.flowpath)
-                where p.scenario = 0
-                ORDER by management desc, flowpath, fpath),
-            agg as (select ch, huc_12, count(*) from data 
-                WHERE ch != '0' GROUP by ch, huc_12),
-            agg2 as (select huc_12, ch, rank() OVER (PARTITION by huc_12
-                ORDER by count DESC) from agg),
-            agg3 as (select huc_12, ch from agg2 where rank = 1)
-            
-            SELECT ST_Transform(simple_geom, 4326) as geom, a.huc_12,
-            a.ch::int as ch
-            from agg3 a JOIN huc12 h on (a.huc_12 = h.huc_12)
-            WHERE h.scenario = 0
+            SELECT ST_Transform(simple_geom, 4326) as geom, dominant_tillage,
+            huc_12
+            from huc12 WHERE scenario = 0
         """,
             conn,
             geom_col="geom",
@@ -57,7 +45,7 @@ def main(argv):
     cmap = get_cmap("plasma")
     bins = list(range(1, 8))
     norm = mpcolors.BoundaryNorm(bins, cmap.N)
-    df["color"] = [c for c in cmap(norm(df["ch"].values))]
+    df["color"] = [c for c in cmap(norm(df["dominant_tillage"].values))]
 
     df.to_crs(mp.panels[0].crs).plot(
         ax=mp.panels[0].ax,