Merge pull request dailyerosion#276 from akrherz/240903

Dynamic Tillage Omnibus

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.6.3"
+    rev: "v0.6.4"
     hooks:
       - id: ruff
         args: [--fix, --exit-non-zero-on-fix]

scripts/dynamic_tillage/compute_smstate.py

@@ -34,12 +34,16 @@ def job(dates, tmpdir, huc12) -> int:
         huc12df = pd.read_sql(
             text(
                 """
-                select o.ofe, p.fpath, o.fbndid, g.plastic_limit,
-                p.fpath || '_' || o.ofe as combo, p.huc_12 as huc12,
-                substr(o.landuse, :charat, 1) as crop
-                from flowpaths p, flowpath_ofes o, gssurgo g
-                WHERE o.flowpath = p.fid and p.huc_12 = :huc12
-                and p.scenario = 0 and o.gssurgo_id = g.id
+    select o.ofe, p.fpath, o.fbndid,
+    case when g.plastic_limit < 40 then
+        g.plastic_limit else
+        g.wepp_min_sw + (g.wepp_max_sw - g.wepp_min_sw) * 0.42
+    end as plastic_limit,
+    p.fpath || '_' || o.ofe as combo, p.huc_12 as huc12,
+    substr(o.landuse, :charat, 1) as crop
+    from flowpaths p, flowpath_ofes o, gssurgo g
+    WHERE o.flowpath = p.fid and p.huc_12 = :huc12
+    and p.scenario = 0 and o.gssurgo_id = g.id
             """
             ),
             conn,

scripts/dynamic_tillage/huc12_timeseries.py

@@ -19,10 +19,9 @@ def get_plastic_limit(huc12: str, year: int) -> pd.DataFrame:
             text(
                 """
                 select o.ofe, p.fpath, o.fbndid,
-                coalesce(
-                    (g.wepp_max_sw + g.wepp_min_sw) / 2.0,
-                    g.plastic_limit
-                ) as plastic_limit,
+                g.wepp_min_sw + (g.wepp_max_sw - g.wepp_min_sw) * 0.58
+                as plastic_limit58,
+                g.plastic_limit,
                 p.fpath || '_' || o.ofe as combo,
                 substr(landuse, :charat, 1) as crop
                 from flowpaths p, flowpath_ofes o, gssurgo g
@@ -72,10 +71,10 @@ def main(huc12: str, year: int):
     huc12sm["below_pl"] = huc12sm["sw1"] < huc12sm["pl0.8"]
 
     fig = figure(
-        title=f"DEP: {huc12} Brush Creek Dynamic Tillage Diagnostic",
+        title=f"DEP: {huc12} Dynamic Tillage Diagnostic",
         subtitle=f"11 April - 15 June {year}",
         logo="dep",
-        figsize=(8, 8),
+        figsize=(8, 10),
     )
     yaxsize = 0.15
     # Scatter plot of Soil Moisture
@@ -147,7 +146,7 @@ def main(huc12: str, year: int):
                 color="red" if row[f"limited_by_{label}"] else "green",
             )
 
-    fig.savefig(f"plotsv2/sm_{huc12}_{year}.png")
+    fig.savefig(f"plotsv3/sm_{huc12}_{year}.png")
 
 
 if __name__ == "__main__":

scripts/dynamic_tillage/nass_comparison.py

@@ -187,7 +187,8 @@ def get_labels(year):
     default="corn",
     help="Crop to Process",
 )
-def main(year, district, state, crop):
+@click.option("--plotv2", is_flag=True, help="Plot v2")
+def main(year, district, state, crop, plotv2: bool):
     """Go Main Go."""
     nass = get_nass(year, state, district, crop)
     fields = get_fields(year, district, state, crop)
@@ -265,7 +266,7 @@ def main(year, district, state, crop):
     # axes showing the days suitable as little boxes for each seven day period
     # and the value of the days suitable for that period in the middle
     ax3 = fig.add_axes([0.1, 0.25, 0.8, 0.07])
-    for valid, row in nass.iterrows():
+    for valid, row in nass[nass["days suitable"].notna()].iterrows():
         ax3.add_patch(
             Rectangle(
                 (valid - pd.Timedelta(days=7), 0),
@@ -320,8 +321,18 @@ def main(year, district, state, crop):
     ax2.plot(
         accum.index,
         accum["acres"] / fields["acres"].sum() * 100.0,
-        label="DEP DynTillv2",
+        label="DEP DynTillv3",
     )
+    if plotv2:
+        v2df = pd.read_csv(
+            f"plotsv2/{crop}_{year}_{state}.csv", parse_dates=["valid"]
+        )
+        v2df["valid"] = v2df["valid"].dt.date
+        ax2.plot(
+            v2df["valid"].values,
+            v2df[f"dep_{crop}_planted"],
+            label="DEP DynTillv2",
+        )
     ax2.scatter(
         nass.index.values,
         nass[f"{crop} planted"],
@@ -334,7 +345,7 @@ def main(year, district, state, crop):
 
     # create a table in the lower right corner with the values from
     txt = ["Weekly Progress", "Date   NASS DEP"]
-    for valid, row in nass.iterrows():
+    for valid, row in nass[nass["days suitable"].notna()].iterrows():
         if valid not in accum.index:
             continue
         dep = accum.at[valid, "percent"]
@@ -360,14 +371,14 @@ def main(year, district, state, crop):
     ax3.set_xlim(*ax2.get_xlim())
 
     nass.to_csv(
-        f"plotsv2/{crop}_{year}_"
+        f"plotsv3/{crop}_{year}_"
         f"{district if district is not None else state}.csv",
         float_format="%.2f",
     )
 
     ss = f"state_{state}"
     fig.savefig(
-        f"plotsv2/{crop}_progress_{year}_"
+        f"plotsv3/{crop}_progress_{year}_"
         f"{district if district is not None else ss}.png"
     )
 

scripts/dynamic_tillage/workflow.py

@@ -24,7 +24,7 @@
 import numpy as np
 import pandas as pd
 from pyiem.database import get_dbconnc, get_sqlalchemy_conn
-from pyiem.iemre import SOUTH, WEST, daily_offset, get_daily_mrms_ncname
+from pyiem.iemre import daily_offset, get_daily_mrms_ncname
 from pyiem.util import archive_fetch, logger, ncopen
 from sqlalchemy import text
 from tqdm import tqdm
@@ -123,11 +123,14 @@ def do_huc12(dt, huc12) -> Tuple[int, int]:
     """
     dbcolidx = dt.year - 2007 + 1
     crops = ["B", "C", "L", "W"]
+    tillagerate = 0.10
     maxrate = 0.06
     # If before April 20, we only plant corn
     if f"{dt:%m%d}" < "0420":
         crops = ["C"]
         maxrate = 0.03  # Life choice
+    elif f"{dt:%m%d}" > "0510":
+        maxrate = 0.10
     with get_sqlalchemy_conn("idep") as conn:
         # build up the cross reference of everyhing we need to know
         df = pd.read_sql(
@@ -208,7 +211,7 @@ def do_huc12(dt, huc12) -> Tuple[int, int]:
 
     total_acres = fields["acres"].sum()
     # NB: Crude assessment of NASS peak daily planting rate, was 10%
-    limit = (total_acres * maxrate) if not mud_it_in else total_acres + 1
+    limit = (total_acres * tillagerate) if not mud_it_in else total_acres + 1
 
     # Work on tillage first, so to avoid planting on tilled fields
     for fbndid, row in fields[fields["till_needed"]].iterrows():
@@ -221,6 +224,8 @@ def do_huc12(dt, huc12) -> Tuple[int, int]:
         if acres_tilled > limit:
             break
 
+    # Redine limit for planting
+    limit = (total_acres * maxrate) if not mud_it_in else total_acres + 1
     # Now we need to plant
     for fbndid, row in fields[fields["plant_needed"]].iterrows():
         # We can't plant fields that were tilled or need tillage GH251
@@ -306,9 +311,11 @@ def estimate_rainfall(huc12df: pd.DataFrame, dt: date):
     LOG.info("Using %s[tidx:%s] for precip", ncfn, tidx)
     with ncopen(ncfn) as nc:
         p01d = nc.variables["p01d"][tidx].filled(0)
+        south = nc.variables["lat"][0]
+        west = nc.variables["lon"][0]
     for idx, row in huc12df.iterrows():
-        y = int((row["lat"] - SOUTH) * 100.0)
-        x = int((row["lon"] - WEST) * 100.0)
+        y = int((row["lat"] - south) * 100.0)
+        x = int((row["lon"] - west) * 100.0)
         huc12df.at[idx, "precip_mm"] = p01d[y, x]
 
 
@@ -379,7 +386,7 @@ def main(scenario, dt, huc12, edr, run_prj2wepp):
                 continue
             # Require that 50% of modelled OFEs are 0.8 of plastic limit
             toowet = gdf["sw1"].gt(gdf["plastic_limit"] * 0.8).sum()
-            if toowet > len(gdf.index) / 2:
+            if toowet > (len(gdf.index) * 0.9):
                 huc12df.at[huc12, "limited_by_soilmoisture"] = True
 
     # restrict to HUC12s that are not limited