Merge pull request dailyerosion#271 from akrherz/240829-2

Omnibus

scripts/dynamic_tillage/assemble_csvs.py

@@ -8,15 +8,15 @@
 def main():
     """Go main Go."""
     dfs = []
-    for csvfn in glob.glob("plotsv2/*_??.csv"):
-        crop, _year, state = csvfn[:-4].split("_")
-        if state not in ["MN", "IA"]:
+    for csvfn in glob.glob("plotsv2/corn*.csv"):
+        crop, _year, datum = csvfn.split("/")[-1][:-4].split("_")
+        if datum in ["IA", "MN", "KS", "NE"]:
             continue
         progress = pd.read_csv(csvfn, parse_dates=["valid"])
-        progress["state"] = state
+        progress["district"] = datum
         progress = progress.rename(
             columns={"dep_planted": f"{crop}_dep_planted"}
-        ).set_index(["state", "valid"])
+        ).set_index(["district", "valid"])
         dfs.append(progress)
     jumbo = pd.concat(dfs)
     rectified = jumbo[
@@ -31,8 +31,8 @@ def main():
 
     (
         rectified.reset_index()
-        .sort_values(["state", "valid"])
-        .to_csv("IA_MN_dep_vs_nass_240828.csv", index=False)
+        .sort_values(["district", "valid"])
+        .to_csv("IA_district_dep_vs_nass_240830.csv", index=False)
     )
 
 

scripts/dynamic_tillage/nass_comparison.py

@@ -205,6 +205,9 @@ def main(year, district, state, crop):
     if ldate >= f"{date.today():%Y-%m-%d}":
         ldate = f"{(date.today() - timedelta(days=1)):%Y-%m-%d}"
     accum = accum.reindex(pd.date_range(f"{year}-04-11", ldate)).ffill()
+    nass = nass.reindex(accum.index)
+    nass.index.name = "valid"
+    nass["dep_corn_planted"] = accum["percent"]
 
     if state is not None:
         title = f"state of {state_names[state]}"
@@ -335,7 +338,6 @@ def main(year, district, state, crop):
         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:%b %-2d}:  {val:3.0f} {dep:3.0f}")
     ax2.text(
@@ -353,12 +355,14 @@ def main(year, district, state, crop):
     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"plotsv2/{crop}_{year}_"
-        f"{district if district is not None else state}.csv"
+        f"{district if district is not None else state}.csv",
+        float_format="%.2f",
     )
 
     ss = f"state_{state}"

scripts/dynamic_tillage/stamp_plot.py

@@ -0,0 +1,133 @@
+"""Create a Postage Stamp Plot of all the years of data."""
+
+import glob
+
+import click
+import matplotlib.cm as mpcm
+import matplotlib.colors as mpcolors
+import numpy
+import pandas as pd
+from matplotlib.patches import Rectangle
+from pyiem.plot import figure_axes, get_cmap
+from pyiem.plot.use_agg import plt
+
+
+@click.command()
+@click.option("--crop", default="corn", help="Which crop to process")
+def main(crop: str):
+    """Go main Go."""
+    dfs = []
+    for csvfn in glob.glob(f"plotsv2/{crop}_*.csv"):
+        _, _, datum = csvfn[:-4].split("_")
+        progress = pd.read_csv(csvfn, parse_dates=["valid"])
+        progress["datum"] = datum
+        dfs.append(progress)
+    jumbo = pd.concat(dfs)
+    jumbo["dayofyear"] = jumbo["valid"].dt.dayofyear
+    fig, ax = figure_axes(
+        title="Corn Planting Progress (DEP Lines, NASS Dots)",
+        logo="dep",
+        figsize=(10, 8),
+    )
+    ax.set_position([0.05, 0.07, 0.85, 0.83])
+    # Remove all splines
+    for side in ["top", "right", "left", "bottom"]:
+        ax.spines[side].set_visible(False)
+
+    xtilesize = 1.0 / 18.0
+    ytilesize = 1.0 / 11.0
+    cmap = get_cmap("RdYlGn_r")
+    norm = mpcolors.BoundaryNorm([0, 10, 40, 70, 100], cmap.N)
+    for x, year in enumerate(range(2007, 2025)):
+        # Psuedo axis definition
+        x0 = x / 18.0
+        xmin = pd.Timestamp(f"{year}-04-11").dayofyear
+        xmax = pd.Timestamp(f"{year}-06-15").dayofyear
+        for y, (district, _color) in enumerate(
+            zip(
+                "sw sc se wc c ec nw nc ne IA MN".split(),
+                (
+                    "red blue green orange tan purple "
+                    "brown pink skyblue black gray"
+                ).split(),
+            )
+        ):
+            # Psuedo axis definition
+            y0 = y / 11.0
+            df = jumbo[
+                (jumbo["valid"].dt.year == year)
+                & (jumbo["datum"] == district)
+                & (jumbo["dayofyear"] >= xmin)
+                & (jumbo["dayofyear"] <= xmax)
+            ]
+            nass = df[df["corn planted"].notna()]
+            # Calculate RMSE
+            rmse = (
+                (nass["corn planted"] - nass["dep_corn_planted"]).pow(2).mean()
+            ) ** 0.5
+            print(f"{year} {district} RMSE: {rmse:.2f}")
+            background_color = cmap(norm(rmse))
+            ax.add_patch(
+                Rectangle(
+                    (x0, y0),
+                    xtilesize,
+                    ytilesize,
+                    color=background_color,
+                    alpha=0.5,
+                )
+            )
+            ax.plot(
+                x0 + (df["dayofyear"] - xmin) / (xmax - xmin) * xtilesize,
+                y0 + df["dep_corn_planted"] / 100.0 * ytilesize,
+                color="k",
+            )
+            ax.scatter(
+                x0 + (nass["dayofyear"] - xmin) / (xmax - xmin) * xtilesize,
+                y0 + nass["corn planted"] / 100.0 * ytilesize,
+                color="k",
+                s=5,
+            )
+            ax.text(
+                x0 + 0.5 * xtilesize,
+                y0 + 0.5 * ytilesize,
+                f"{rmse:.0f}",
+                ha="center",
+                va="center",
+                color="b",
+                fontsize=14,
+            )
+
+    ax.set_xticks(numpy.arange(0, 0.99, 1.0 / 18.0) + xtilesize / 2.0)
+    ax.set_xticklabels(range(2007, 2025), rotation=45)
+    ax.set_yticks(numpy.arange(0, 0.99, 1.0 / 11.0) + ytilesize / 2.0)
+    ax.set_yticklabels(
+        (
+            "IA\nSW IA\nSC IA\nSE IA\nWC IA\nC IA\nEC "
+            "IA\nNW IA\nNC IA\nNE IA MN"
+        ).split(" "),
+        rotation=45,
+        ha="right",
+    )
+    # Add gridlines
+    for x in numpy.arange(0, 1.1, 1.0 / 18.0):
+        ax.axvline(x, color="k", lw=0.5)
+    for y in numpy.arange(0, 1.01, 1.0 / 11.0):
+        ax.axhline(y, color="k", lw=0.5)
+    ax.set_xlim(-0.01, 1.01)
+    ax.set_ylim(-0.01, 1.01)
+
+    # Add colorbar
+    ax2 = fig.add_axes([0.9, 0.1, 0.02, 0.8])
+    cb = plt.colorbar(
+        mpcm.ScalarMappable(norm=norm, cmap=cmap),
+        cax=ax2,
+        orientation="vertical",
+        alpha=0.5,
+    )
+    cb.set_label("RMSE of NASS vs DEP Planting Progress [%]")
+
+    fig.savefig("test.png")
+
+
+if __name__ == "__main__":
+    main()

scripts/dynamic_tillage/workflow.py

@@ -372,6 +372,8 @@ def main(scenario, dt, huc12, edr, run_prj2wepp):
             f"/mnt/idep2/data/smstate/{yesterday:%Y}/"
             f"smstate{yesterday:%Y%m%d}.feather"
         )
+        # Only consider rows that are in either Corn or Soybean
+        smdf = smdf[smdf["crop"].isin(["C", "S"])]
         for huc12, gdf in smdf.groupby("huc12"):
             if huc12 not in huc12df.index:
                 continue

scripts/import/flowpath2prj.py

@@ -7,6 +7,7 @@
   F - forest            F1          F25    IniCropDef.Tre_2239
   G - Sorghum
   P - Pasture           P1          P25    IniCropDef.gra_3425
+      (If all P, then pasture.  If P intermixed, then Alfalfa see GH270)
   C - Corn              C1          C25    IniCropDef.Default
   R - Other crops (Barley?)  R1          R25    IniCropDef.Aft_12889
   T - Water (could have flooded out for one year, wetlands)
@@ -75,7 +76,7 @@
 }
 
 
-def do_rotation(scenario, zone, rotfn, landuse, management):
+def do_rotation(scenario, zone, rotfn, landuse: str, management):
     """Create the rotation file.
 
     Args:
@@ -92,13 +93,22 @@ def do_rotation(scenario, zone, rotfn, landuse, management):
     data = {"yearly": ""}
     data["name"] = f"{landuse}-{management}"
     # Special hack for forest
-    if landuse[0] == "F" and landuse == len(landuse) * landuse[0]:
+    if landuse[0] == "F" and len(set(landuse)) == 1:
         data["initcond"] = FOREST[management[0]]
         for i in range(1, YEARS):
             # Reset roughness each year
             data["yearly"] += (
                 f"1 2 {i} 1 Tillage OpCropDef.Old_6807 {{0.001, 2}}\n"
             )
+    # Special hack for pasture
+    elif landuse[0] == "P" and len(set(landuse)) == 1:
+        data["initcond"] = INITIAL_COND[landuse[0]]
+        # For WEPP purposes and to our knowledge, we still need to initially
+        # do a light tillage and plant this pasture to alfalfa
+        data["yearly"] += (
+            "4  14 1 1 Tillage OpCropDef.FCSTACDP      {0.101600, 2}\n"
+            "4  15 1 1 Plant-Perennial CropDef.ALFALFA  {0.000000}\n"
+        )
     else:
         data["initcond"] = INITIAL_COND.get(landuse[0], INITIAL_COND_DEFAULT)
         prevcode = "C" if landuse[0] not in INITIAL_COND else landuse[0]