Merge pull request dailyerosion#221 from akrherz/240408

Omnibus

pyproject.toml

@@ -9,3 +9,8 @@ lint.select = [
 ]
 target-version = "py39"
 
+# per memory issues microsoft/pylance-release/issues/5739
+[tool.pyright]
+autoSearchPaths = false
+exclude = ["data/", "prj2wepp/"]
+ignore = ["data/", "prj2wepp/"]

scripts/dynamic_tillage/accum_acres_plot.py

@@ -222,5 +222,5 @@ def main():
 
 
 if __name__ == "__main__":
-    for _i, _dt in enumerate(pd.date_range("2023-04-15", "2023-06-04")):
+    for _i, _dt in enumerate(pd.date_range("2008-04-15", "2008-06-05")):
         plot_map_progress(_i, _dt)

scripts/dynamic_tillage/fourpanel_limiter.py

@@ -0,0 +1,68 @@
+"""A four panel plot diagnostic showing the limiting factor."""
+
+import click
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+from pyiem.database import get_sqlalchemy_conn
+from pyiem.plot import figure
+from pyiem.util import load_geodf
+
+
+@click.command()
+@click.option("--date", "dt", type=click.DateTime(), help="Date to plot.")
+def main(dt):
+    """Go Main Go."""
+    with get_sqlalchemy_conn("idep") as conn:
+        huc12df = gpd.read_postgis(
+            "SELECT huc_12, st_transform(simple_geom, 4326) as geom "
+            "from huc12 where scenario = 0",
+            conn,
+            geom_col="geom",
+            index_col="huc_12",
+        )
+
+    huc12data = pd.read_feather(
+        f"/mnt/idep2/data/huc12status/{dt:%Y}/{dt:%Y%m%d}.feather"
+    )
+    huc12df = huc12df.join(huc12data, how="inner")
+
+    states = load_geodf("us_states", 4326)
+
+    fig = figure(
+        title=f"{dt:%d %b %Y} Dynamic Tillage Limiting Factor",
+        subtitle="Red indicates limiting factor",
+        logo="dep",
+        figsize=(10.24, 7.68),
+    )
+    opts = [
+        [0.03, 0.03, "Plastic Limit", "limited_by_soilmoisture"],
+        [0.5, 0.03, "Combined Limit", "limited"],
+        [0.03, 0.48, "Soil Temp Limit", "limited_by_soiltemp"],
+        [0.5, 0.48, "Rainfall Limit", "limited_by_precip"],
+    ]
+    for opt in opts:
+        ax = fig.add_axes([opt[0], opt[1], 0.42, 0.38])
+        ax.set_xticklabels([])
+        ax.set_yticklabels([])
+        title = (
+            f"{opt[2]} {huc12df[opt[3]].sum()}/{len(huc12df)} "
+            f"{huc12df[opt[3]].sum() / len(huc12df):.1%}"
+        )
+        ax.set_title(title, fontsize=16)
+        # plot in red when limited is true, else blue
+        huc12df.plot(
+            ax=ax,
+            aspect=None,
+            color=np.where(huc12df[opt[3]], "r", "b"),
+        )
+        xlim = ax.get_xlim()
+        ylim = ax.get_ylim()
+        states.plot(ax=ax, color="None", edgecolor="k", lw=0.5, aspect=None)
+        ax.set_xlim(xlim)
+        ax.set_ylim(ylim)
+    fig.savefig(f"huc12limiter_{dt:%Y%m%d}.png")
+
+
+if __name__ == "__main__":
+    main()

scripts/dynamic_tillage/nass_comparison.py

@@ -0,0 +1,131 @@
+"""Plot our planting progress vs that of NASS."""
+
+import click
+import geopandas as gpd
+import pandas as pd
+from pyiem.database import get_sqlalchemy_conn
+from pyiem.plot import figure_axes
+from sqlalchemy import text
+
+XREF = {
+    "nw": "IAC001",
+    "nc": "IAC002",
+    "ne": "IAC003",
+    "wc": "IAC004",
+    "c": "IAC005",
+    "ec": "IAC006",
+    "sw": "IAC007",
+    "sc": "IAC008",
+    "se": "IAC009",
+}
+
+
+@click.command()
+@click.option("--year", type=int, help="Year to plot")
+@click.option("--district", type=str, help="NASS District")
+def main(year, district):
+    """Go Main Go."""
+    charidx = year - 2007 + 1
+    with get_sqlalchemy_conn("coop") as conn:
+        climdiv = gpd.read_postgis(
+            text(
+                """
+                select t.iemid, r.geom from climodat_regions r JOIN stations t
+                on (r.iemid = t.iemid) where t.id = :sid
+                """
+            ),
+            conn,
+            params={"sid": XREF[district]},
+            geom_col="geom",
+            index_col="iemid",
+        )
+
+    # Get the NASS data
+    with get_sqlalchemy_conn("coop") as conn:
+        nass = pd.read_sql(
+            text(
+                """
+            select * from nass_iowa where metric = 'corn planted'
+            and extract(year from valid) = :year ORDER by valid ASc
+            """
+            ),
+            conn,
+            params={"year": year},
+        )
+
+    # Figure out the planting dates
+    with get_sqlalchemy_conn("idep") as conn:
+        fields = gpd.read_postgis(
+            text("""
+            select st_transform(geom, 4326) as geom, plant, huc12, fbndid,
+            acres from fields f LEFT JOIN field_operations o
+            on (f.field_id = o.field_id and o.year = :year)
+            where scenario = 0 and
+            substr(landuse, :charidx, 1) = 'C'
+            """),
+            conn,
+            params={"year": year, "charidx": charidx},
+            geom_col="geom",
+            parse_dates="plant",
+        )
+
+    # Spatially filter fields that are inside the climdiv region
+    fields = gpd.sjoin(fields, climdiv, predicate="intersects")
+
+    # accumulate acres planted by date
+    accum = fields[["plant", "acres"]].groupby("plant").sum().cumsum()
+    accum["percent"] = accum["acres"] / fields["acres"].sum() * 100.0
+    accum = accum.reindex(
+        pd.date_range(f"{year}-04-15", f"{year}-06-23")
+    ).ffill()
+
+    fig, ax = figure_axes(
+        logo="dep",
+        title=(
+            f"{year} Corn Planting Progress vs "
+            f"NASS Iowa Crop District {district.upper()}"
+        ),
+        figsize=(10.24, 7.68),
+    )
+    ax.set_ylim(0, 100)
+    ax.set_ylabel("Percent of Acres Planted")
+    ax.set_xlabel("Date")
+    ax.grid(True)
+    ax.plot(
+        accum.index,
+        accum["acres"] / fields["acres"].sum() * 100.0,
+        label="DEP",
+    )
+
+    ax.scatter(
+        nass["valid"],
+        nass[district],
+        marker="*",
+        s=30,
+        color="r",
+        label="NASS",
+    )
+    ax.legend(loc=2)
+
+    # create a table in the lower right corner with the values from
+    txt = ["Date       | NASS | DEP"]
+    for _, row in nass.iterrows():
+        dep = accum.at[pd.Timestamp(row["valid"]), "percent"]
+        txt.append(f"{row['valid']} | {row['se']:.0f} | {dep:.1f}")
+    print("\n".join(txt))
+    ax.text(
+        0.99,
+        0.01,
+        "\n".join(txt),
+        transform=ax.transAxes,
+        fontsize=10,
+        va="bottom",
+        ha="right",
+        bbox=dict(facecolor="white", edgecolor="black", pad=5),
+    )
+
+    fig.savefig(f"corn_progress_{year}_{district}.png")
+
+
+if __name__ == "__main__":
+    main()

scripts/dynamic_tillage/workflow.py

@@ -31,6 +31,7 @@
 pd.set_option("future.no_silent_downcasting", True)
 LOG = logger()
 CPU_COUNT = min([4, cpu_count() / 2])
+HUC12STATUSDIR = "/mnt/idep2/data/huc12status"
 PROJDIR = "/opt/dep/prj2wepp"
 EXE = f"{PROJDIR}/prj2wepp"
 
@@ -112,6 +113,7 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]:
     Returns the number of acres planted and tilled.  A negative value is a
     sentinel that the HUC12 failed due to soil moisture constraint.
     """
+    dbcolidx = dt.year - 2007 + 1
     with get_sqlalchemy_conn("idep") as conn:
         # build up the cross reference of everyhing we need to know
         df = pd.read_sql(
@@ -126,7 +128,9 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]:
                 SELECT ofe, fpath, f.fbndid, plastic_limit, f.landuse,
                 f.management, f.acres, f.field_id
                 from fields f LEFT JOIN myofes m on (f.fbndid = m.fbndid)
-                WHERE f.huc12 = :huc12 and f.isag ORDER by fpath, ofe)
+                WHERE f.huc12 = :huc12 and
+                substr(landuse, :dbcolidx, 1) = ANY(:crops)
+                ORDER by fpath, ofe)
             select a.*, o.till1, o.till2, o.till3, o.plant, o.year,
             plant >= :dt as plant_needed,
             coalesce(
@@ -138,7 +142,13 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]:
             """
             ),
             conn,
-            params={"huc12": huc12, "year": dt.year, "dt": dt},
+            params={
+                "huc12": huc12,
+                "year": dt.year,
+                "dt": dt,
+                "crops": ["B", "C", "L", "W"],
+                "dbcolidx": dbcolidx,
+            },
             index_col=None,
         )
     if df.empty:
@@ -171,15 +181,16 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]:
         # bit generous here, are 50% of **model** fields above plastic_limit?
         fields["sw1"] = pd.Series(fieldsw)
         modelled = (fields["sw1"] > 0).sum()
-        hits = (fields["sw1"] > fields["plastic_limit"]).sum()
-        if hits > modelled / 2:
-            LOG.debug(
-                "HUC12: %s has %s/%s fields above plastic limit",
-                huc12,
-                hits,
-                modelled,
-            )
-            return -1, -1
+        if modelled > 1:
+            hits = (fields["sw1"] > fields["plastic_limit"]).sum()
+            if hits > modelled / 2:
+                LOG.debug(
+                    "HUC12: %s has %s/%s fields above plastic limit",
+                    huc12,
+                    hits,
+                    modelled,
+                )
+                return -1, -1
 
     # Compute things for year
     char_at = (dt.year - 2007) + 1
@@ -244,7 +255,7 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]:
 
 def estimate_soiltemp(huc12df, dt):
     """Write mean GFS soil temperature C into huc12df."""
-    ppath = f"{dt:%Y/%m/%d}/model/gfs/gfs_{dt:%Y%m%d}00_iemre.nc"
+    ppath = f"{dt:%Y/%m/%d}/model/gfs/gfs_{dt:%Y%m%d}12_iemre.nc"
     with archive_fetch(ppath) as ncfn, ncopen(ncfn) as nc:
         tsoil = np.mean(nc.variables["tsoil"][1:7], axis=0) - 273.15
         y = np.digitize(huc12df["lat"].values, nc.variables["lat"][:])
@@ -349,6 +360,7 @@ def main(scenario, dt, huc12):
             progress.update()
             if planted < 0:
                 huc12df.at[huc12, "limited_by_soilmoisture"] = True
+                huc12df.at[huc12, "limited"] = True
             else:
                 total_planted += planted
                 total_tilled += tilled
@@ -365,6 +377,13 @@ def main(scenario, dt, huc12):
         huc12df["limited_by_soiltemp"].sum(),
     )
 
+    # save a debug output file
+    mydir = f"{HUC12STATUSDIR}/{dt:%Y}"
+    os.makedirs(mydir, exist_ok=True)
+    if len(huc12df.index) == 1:  # Don't save when run for single huc12
+        return
+    huc12df.drop(columns="geom").to_feather(f"{mydir}/{dt:%Y%m%d}.feather")
+
 
 if __name__ == "__main__":
     main()