Merge pull request dailyerosion#267 from akrherz/240823

Omnibus

scripts/dynamic_tillage/assemble_csvs.py

@@ -8,19 +8,12 @@
 def main():
     """Go main Go."""
     dfs = []
-    for csvfn in glob.glob("*_??.csv"):
+    for csvfn in glob.glob("plotsv2/*_??.csv"):
         crop, _year, state = csvfn[:-4].split("_")
         if state not in ["MN", "IA"]:
             continue
         progress = pd.read_csv(csvfn, parse_dates=["valid"])
         progress["state"] = state
-        # Need to fill out NASS with 100s
-        firstidx = progress[f"{crop} planted"].first_valid_index()
-        if firstidx > 0:
-            progress.loc[:firstidx, f"{crop} planted"] = 0
-        lastidx = progress[f"{crop} planted"].last_valid_index()
-        if lastidx is not None:
-            progress.loc[(lastidx + 1) :, f"{crop} planted"] = 100
         progress = progress.rename(
             columns={"dep_planted": f"{crop}_dep_planted"}
         ).set_index(["state", "valid"])
@@ -29,24 +22,17 @@ def main():
     rectified = jumbo[
         jumbo["corn planted"].notna() | jumbo["dep_corn_planted"].notna()
     ].copy()
-    soybeans = jumbo[
-        jumbo["soybeans planted"].notna()
-        | jumbo["dep_soybeans_planted"].notna()
-    ].copy()
-    rectified["soybeans planted"] = soybeans["soybeans planted"]
-    rectified["dep_soybeans_planted"] = soybeans["dep_soybeans_planted"]
 
     # Drop things with all nulls
     rectified = rectified[
         rectified["dep_corn_planted"].notna()
-        | rectified["dep_soybeans_planted"].notna()
         | rectified["dep_days_suitable"].notna()
     ]
 
     (
         rectified.reset_index()
         .sort_values(["state", "valid"])
-        .to_csv("IA_MN_dep_vs_nass_240715.csv", index=False)
+        .to_csv("IA_MN_dep_vs_nass_240828.csv", index=False)
     )
 
 

scripts/dynamic_tillage/fourpanel_limiter.py

@@ -10,7 +10,9 @@
 
 
 @click.command()
-@click.option("--date", "dt", type=click.DateTime(), help="Date to plot.")
+@click.option(
+    "--date", "dt", required=True, type=click.DateTime(), help="Date to plot."
+)
 def main(dt):
     """Go Main Go."""
     with get_sqlalchemy_conn("idep") as conn:
@@ -61,7 +63,7 @@ def main(dt):
         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")
+    fig.savefig(f"plotsv2/huc12limiter_{dt:%Y%m%d}.png")
 
 
 if __name__ == "__main__":

scripts/dynamic_tillage/huc12_timeseries.py

@@ -10,15 +10,21 @@
 from sqlalchemy import text
 
 
-def get_plastic_limit(huc12: str) -> pd.DataFrame:
+def get_plastic_limit(huc12: str, year: int) -> pd.DataFrame:
     """Figure out what the plastic limit is."""
+    charat = year - 2007
     with get_sqlalchemy_conn("idep") as conn:
         # build up the cross reference of everyhing we need to know
         return pd.read_sql(
             text(
                 """
-                select o.ofe, p.fpath, o.fbndid, g.plastic_limit,
-                p.fpath || '_' || o.ofe as combo
+                select o.ofe, p.fpath, o.fbndid,
+                coalesce(
+                    (g.wepp_max_sw + g.wepp_min_sw) / 2.0,
+                    g.plastic_limit
+                ) as plastic_limit,
+                p.fpath || '_' || o.ofe as combo,
+                substr(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
@@ -27,20 +33,21 @@ def get_plastic_limit(huc12: str) -> pd.DataFrame:
             conn,
             params={
                 "huc12": huc12,
+                "charat": charat,
             },
             index_col=None,
         )
 
 
 @click.command()
-@click.option("--huc12", help="HUC12 To Plot")
-@click.option("--year", type=int, help="Year to Plot")
-def main(huc12, year):
+@click.option("--huc12", required=True, help="HUC12 To Plot")
+@click.option("--year", required=True, type=int, help="Year to Plot")
+def main(huc12: str, year: int):
     """Go Main Go."""
-    pldf = get_plastic_limit(huc12)
+    pldf = get_plastic_limit(huc12, year)
     smdfs = []
     hucstatusdfs = []
-    for dt in pd.date_range(f"{year}/04/14", f"{year}/05/31"):
+    for dt in pd.date_range(f"{year}/04/11", f"{year}/06/13"):
         fn = f"/mnt/idep2/data/smstate/{dt:%Y}/smstate{dt:%Y%m%d}.feather"
         if not os.path.isfile(fn):
             continue
@@ -54,24 +61,35 @@ def main(huc12, year):
         statusdf["date"] = dt
         hucstatusdfs.append(statusdf.loc[[huc12]])
     hucstatusdf = pd.concat(hucstatusdfs)
-    print(hucstatusdf)
     huc12sm = pd.concat(smdfs)
     huc12sm["combo"] = (
         huc12sm["fpath"].astype(str) + "_" + huc12sm["ofe"].astype(str)
     )
-    huc12sm["plastic_limit"] = (
+    huc12sm["pl0.8"] = (
         huc12sm["combo"].map(pldf.set_index("combo")["plastic_limit"]) * 0.8
     )
-    huc12sm["below_pl"] = huc12sm["sw1"] < huc12sm["plastic_limit"]
+    huc12sm["crop"] = huc12sm["combo"].map(pldf.set_index("combo")["crop"])
+    huc12sm["below_pl"] = huc12sm["sw1"] < huc12sm["pl0.8"]
 
     fig = figure(
-        title=f"DEP: {huc12} Walnut Creek Dynamic Tillage Diagnostic",
-        subtitle="15 April - 23 April 2024",
+        title=f"DEP: {huc12} Brush Creek Dynamic Tillage Diagnostic",
+        subtitle=f"11 April - 15 June {year}",
         logo="dep",
-        figsize=(8, 6),
+        figsize=(8, 8),
     )
-    yaxsize = 0.2
-    ax = fig.add_axes([0.1, 0.68, 0.8, yaxsize])
+    yaxsize = 0.15
+    # Scatter plot of Soil Moisture
+    ax0 = fig.add_axes([0.1, 0.75, 0.8, yaxsize])
+    past = huc12sm[huc12sm["crop"] == "P"]
+    ax0.scatter(past["date"], past["sw1"] - past["pl0.8"], c="r", s=2)
+    nonp = huc12sm[huc12sm["crop"] != "P"]
+    ax0.scatter(nonp["date"], nonp["sw1"] - nonp["pl0.8"], c="b", s=2)
+    ax0.xaxis.set_major_formatter(DateFormatter("%-d %b"))
+    ax0.xaxis.set_major_formatter(DateFormatter("%-d %b"))
+    ax0.grid(True)
+
+    # Percent of OFEs below 0.8*PL
+    ax = fig.add_axes([0.1, 0.51, 0.8, yaxsize], sharex=ax0)
     ax.text(
         0, 1, "* SM Evaluation using previous date", transform=ax.transAxes
     )
@@ -89,9 +107,8 @@ def main(huc12, year):
     ax.set_yticks([0, 10, 25, 50, 75, 90, 100])
     ax.grid(True)
     ax.set_ylabel("Percent of OFEs\nBelow 0.8*PL")
-    ax.xaxis.set_major_formatter(DateFormatter("%-d %b"))
 
-    ax2 = fig.add_axes([0.1, 0.4, 0.8, yaxsize], sharex=ax)
+    ax2 = fig.add_axes([0.1, 0.3, 0.8, yaxsize], sharex=ax0)
     ax2.bar(
         hucstatusdf["date"],
         hucstatusdf["precip_mm"],
@@ -103,7 +120,7 @@ def main(huc12, year):
     ax2.axhline(10, lw=2, color="k")
     ax2.grid(True)
 
-    ax3 = fig.add_axes([0.1, 0.1, 0.8, yaxsize], sharex=ax)
+    ax3 = fig.add_axes([0.1, 0.05, 0.8, yaxsize], sharex=ax0)
     ax3.bar(
         hucstatusdf["date"],
         hucstatusdf["tsoil_avg"],
@@ -130,7 +147,7 @@ def main(huc12, year):
                 color="red" if row[f"limited_by_{label}"] else "green",
             )
 
-    fig.savefig(f"sm_{huc12}_{year}.png")
+    fig.savefig(f"plotsv2/sm_{huc12}_{year}.png")
 
 
 if __name__ == "__main__":

scripts/dynamic_tillage/nass_comparison.py

@@ -5,6 +5,7 @@
 
 import click
 import geopandas as gpd
+import matplotlib.colors as mpcolors
 import numpy as np
 import pandas as pd
 from matplotlib.patches import Rectangle
@@ -123,6 +124,16 @@ def get_nass(year, state, district, crop) -> pd.DataFrame:
             )
     nass["valid"] = pd.to_datetime(nass["valid"])
     nass = nass.pivot(index="valid", columns="metric", values="num_value")
+    # NASS does not always include 0% and 100% in their data, so we need to
+    # fill in the blanks
+    col = f"{crop} planted"
+    if nass[col].min() > 0:
+        firstrow = nass[nass[col] > 0].index[0]
+        nass.loc[firstrow - pd.Timedelta(days=7), col] = 0
+    # TODO, this may not always be right for in-progress years
+    if nass[col].max() < 100:
+        lastrow = nass[nass[col] > 0].index[-1]
+        nass.loc[lastrow + pd.Timedelta(days=7), col] = 100
     # Create space for DEP to store things
     nass[f"dep_{crop}_planted"] = np.nan
     nass["dep_days_suitable"] = np.nan
@@ -155,6 +166,17 @@ def get_fields(year, district, state, crop) -> gpd.GeoDataFrame:
     return fields
 
 
+def get_labels(year):
+    """Figure out where we should place xticks."""
+    xticks = []
+    xticklabels = []
+    for dt in pd.date_range(f"{year}-04-11", f"{year}-06-20"):
+        if dt.dayofweek == 6:
+            xticks.append(dt)
+            xticklabels.append(f"{dt:%b %-d}")
+    return xticks, xticklabels
+
+
 @click.command()
 @click.option("--year", type=int, help="Year to plot")
 @click.option("--district", type=str, help="NASS District (lowercase)")
@@ -170,6 +192,8 @@ def main(year, district, state, crop):
     nass = get_nass(year, state, district, crop)
     fields = get_fields(year, district, state, crop)
 
+    xticks, xticklabels = get_labels(year)
+
     # Spatially filter fields that are inside the climdiv region
     daily_limits = compute_limits(fields["huc12"].unique(), year)
 
@@ -210,6 +234,7 @@ def main(year, district, state, crop):
         ]
         i += 1
     cmap = get_cmap("viridis")
+    norm = mpcolors.BoundaryNorm(np.arange(0, 101, 10), cmap.N)
     res = ax.imshow(
         data,
         extent=[
@@ -220,14 +245,15 @@ def main(year, district, state, crop):
         ],
         aspect="auto",
         interpolation="nearest",
-        vmin=0,
-        vmax=100,
+        norm=norm,
         cmap=cmap,
     )
     ax.set_yticks(range(4))
     ax.set_yticklabels(
         ["Combined", "Soil Moisture", "Soil Temp", "Precipitation"]
     )
+    ax.set_xticks(xticks)
+    ax.set_xticklabels(xticklabels)
     ax.text(0.02, 1.05, "Limiting Factors", transform=ax.transAxes)
     cax = fig.add_axes([0.92, 0.05, 0.02, 0.2])
     fig.colorbar(res, cax=cax, label="Percent of HUC12s Limited")
@@ -284,7 +310,8 @@ def main(year, district, state, crop):
 
     # ////////////////////////////////////////////////////////////
     ax2 = fig.add_axes([0.1, 0.4, 0.8, 0.5])
-    ax2.set_ylim(0, 100)
+    ax2.set_ylim(-2, 102)
+    ax2.set_yticks(range(0, 101, 10))
     ax2.set_ylabel("Percent of Acres Planted")
     ax2.grid(True)
     ax2.plot(
@@ -303,36 +330,40 @@ def main(year, district, state, crop):
     ax2.legend(loc=2)
 
     # create a table in the lower right corner with the values from
-    txt = ["Date       | NASS | DEP"]
+    txt = ["Weekly Progress", "Date   NASS DEP"]
     for valid, row in nass.iterrows():
         if valid not in accum.index:
             continue
         dep = accum.at[valid, "percent"]
         nass.at[valid, f"dep_{crop}_planted"] = dep
         val = row[f"{crop} planted"]
-        txt.append(f"{valid:%Y-%m-%d} | {val:.0f} | {dep:.1f}")
+        txt.append(f"{valid:%b %-2d}:  {val:3.0f} {dep:3.0f}")
     ax2.text(
-        0.99,
-        0.01,
+        1.05,
+        0.05,
         "\n".join(txt),
         transform=ax2.transAxes,
         fontsize=10,
+        fontfamily="monospace",
         va="bottom",
         ha="right",
         bbox=dict(facecolor="white", edgecolor="black", pad=5),
     )
+    ax2.set_xticks(xticks)
+    ax2.set_xticklabels(xticklabels)
 
     # Sync up the two xaxes
     ax.set_xlim(*ax2.get_xlim())
     ax3.set_xlim(*ax2.get_xlim())
 
     nass.to_csv(
-        f"{crop}_{year}_{district if district is not None else state}.csv"
+        f"plotsv2/{crop}_{year}_"
+        f"{district if district is not None else state}.csv"
     )
 
     ss = f"state_{state}"
     fig.savefig(
-        f"{crop}_progress_{year}_"
+        f"plotsv2/{crop}_progress_{year}_"
         f"{district if district is not None else ss}.png"
     )
 

src/pydep/tillage.py

@@ -54,7 +54,7 @@
 }
 
 
-def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year):
+def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year: int):
     """Read a block file and do replacements
 
     Args:
@@ -76,9 +76,11 @@ def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year):
     with open(blockfn, "r", encoding="utf8") as fh:
         data = fh.read()
     # Special consideration for planting alfalfa
-    if code == "P" and prevcode != "P":
-        # Best we can do now is plant it on Apr 15, sigh
+    if code == "P" and (prevcode != "P" or year == 1):
+        # Best we can do now is plant it on Apr 15 with a simple disk the
+        # day before
         data = (
+            f"4  14 {year} 1 Tillage OpCropDef.FCSTACDP      {{0.101600, 2}}\n"
             f"4  15 {year} 1 Plant-Perennial CropDef.ALFALFA  {{0.000000}}\n"
             f"{data}"
         )

src/regression_tests/CB_5percent_5till/answer.json

@@ -1,4 +1,5 @@
 {
     "av_det_2012": 4.499,
-    "av_det": 4.264
+    "av_det": 4.264,
+    "max_sw1": 59.61
 }