more tillage timing work #67

scripts/tillage_timing/README.md

@@ -23,17 +23,17 @@ ID | Name | Status
 67 | May 20 Tillage, Constant Climate | done
 68 | May 25 Tillage, Constant Climate | done
 69 | May 30 Tillage, Constant Climate | done
-70 | April 10 Tillage, Local Climate | done
-71 | April 15 Tillage, Local Climate | done
-72 | April 20 Tillage, Local Climate | done
-73 | April 25 Tillage, Local Climate | done
-74 | April 30 Tillage, Local Climate | done
-75 | May 5 Tillage, Local Climate | done
-76 | May 10 Tillage, Local Climate | done
-77 | May 15 Tillage, Local Climate | done
-78 | May 20 Tillage, Local Climate | done
-79 | May 25 Tillage, Local Climate | done
-80 | May 30 Tillage, Local Climate | done
+70 | April 10 Tillage, Local Climate | no wb
+71 | April 15 Tillage, Local Climate | no wb
+72 | April 20 Tillage, Local Climate | no wb
+73 | April 25 Tillage, Local Climate | no wb
+74 | April 30 Tillage, Local Climate | no wb
+75 | May 5 Tillage, Local Climate | no wb
+76 | May 10 Tillage, Local Climate | no wb
+77 | May 15 Tillage, Local Climate | no wb
+78 | May 20 Tillage, Local Climate | no wb
+79 | May 25 Tillage, Local Climate | no wb
+80 | May 30 Tillage, Local Climate | no wb
 81 | Dynamic Tillage April 15 after 45% | done
 82 | Dynamic Tillage April 15 after 40% | done
 83 | Dynamic Tillage April 15 after 35% | done
@@ -44,3 +44,17 @@ ID | Name | Status
 88 | Dynamic Tillage April 15 after Plastic Limit 1.0 | done
 89 | Dynamic Tillage April 15 after Plastic Limit 0.9 | done
 90 | Dynamic Tillage April 15 after Plastic Limit Fx | done
+
+21 Jan 2020 Meeting Notes
+------
+
+- [x] plot averaged 5 day precipitation for the 30 HUC12s
+- [x] plot showing how anticedent and after precip changes by date
+- [ ] Water content on days of tillage for 1a and 1b
+- [ ] Am I plotting gravimetric or VWC for soil moisture?
+- [ ] seperate each MLRA and tillage run
+- [ ] yikes, do dynamic tillage runs for every year
+- [ ] how to compute vulernability windows
+- [ ] can I plot LAI with time and see what WEPP does here
+- [ ] can the window of vulnerability be isolated by date of year?
+- [ ] clarify that dates are generally planting date, not tillage.

scripts/tillage_timing/plot_dailydeltas.py

@@ -0,0 +1,70 @@
+"""A plot of the deltas for erosion between scenarios."""
+import datetime
+import sys
+
+from pyiem.dep import read_env
+from pyiem.plot.use_agg import plt
+import matplotlib.dates as mdates
+
+
+def main(argv):
+    """Go Main Go."""
+    huc12 = argv[1]
+    fpath = argv[2]
+    year = int(argv[3])
+    prop_cycle = plt.rcParams["axes.prop_cycle"]
+    colors = prop_cycle.by_key()["color"]
+
+    data = {}
+    for scenario in range(59, 70):
+        data[scenario] = read_env(
+            "/i/%s/env/%s/%s/%s_%s.env"
+            % (scenario, huc12[:8], huc12[8:], huc12, fpath)
+        ).set_index("date")
+        print(data[scenario]["av_det"].sum())
+    ax = plt.axes([0.2, 0.1, 0.75, 0.75])
+    baseline = data[59][data[59].index.year == year]
+    yticklabels = []
+    for scenario in range(60, 70):
+        color = colors[scenario - 60]
+        date = datetime.date(2000, 4, 15) + datetime.timedelta(
+            days=(scenario - 60) * 5
+        )
+        scendata = data[scenario][data[scenario].index.year == year]
+        delta = scendata["sed_del"] - baseline["sed_del"]
+        delta = delta[delta != 0]
+        total = (
+            (scendata["sed_del"].sum() - baseline["sed_del"].sum())
+            / baseline["sed_del"].sum()
+        ) * 100.0
+        yticklabels.append("%s %4.2f%%" % (date.strftime("%b %d"), total))
+        x = delta.index.to_pydatetime()
+        # res = ax.scatter(x, delta.values + (scenario - 60))
+        for idx, val in enumerate(delta):
+            ax.arrow(
+                x[idx],
+                scenario - 60,
+                0,
+                val,
+                head_width=0.5,
+                head_length=0.1,
+                fc=color,
+                ec=color,
+            )
+        ax.axhline(scenario - 60, color=color)
+    ax.set_xlim(datetime.date(year, 1, 1), datetime.date(year + 1, 1, 1))
+    ax.set_ylim(-0.5, 10)
+    ax.xaxis.set_major_locator(mdates.DayLocator([1]))
+    ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
+    ax.set_title(
+        "huc12: %s fpath: %s\n%s Daily Change in Delivery vs Apr 10 Planting"
+        % (huc12, fpath, year)
+    )
+    ax.grid(axis="x")
+    ax.set_yticks(range(10))
+    ax.set_yticklabels(yticklabels)
+    plt.gcf().savefig("test.png")
+
+
+if __name__ == "__main__":
+    main(sys.argv)

scripts/tillage_timing/plot_dailydeltas_huc12.py

@@ -0,0 +1,74 @@
+"""A plot of the deltas for erosion between scenarios."""
+import datetime
+import sys
+
+from pyiem.plot.use_agg import plt
+from pyiem.util import get_dbconn
+import pandas as pd
+from pandas.io.sql import read_sql
+import matplotlib.dates as mdates
+
+
+def main(argv):
+    """Go Main Go."""
+    huc12 = argv[1]
+    year = int(argv[2])
+    prop_cycle = plt.rcParams["axes.prop_cycle"]
+    colors = prop_cycle.by_key()["color"]
+    pgconn = get_dbconn("idep")
+    df = read_sql(
+        """
+        SELECT scenario, huc_12, avg_delivery, valid from results_by_huc12
+        WHERE scenario >= 59 and scenario < 70 and
+        extract(year from valid) = %s and huc_12 = %s ORDER by valid ASC
+    """,
+        pgconn,
+        params=(year, huc12),
+    )
+    df["valid"] = pd.to_datetime(df["valid"])
+    ax = plt.axes([0.2, 0.1, 0.75, 0.75])
+    baseline = df[df["scenario"] == 59].copy().set_index("valid")
+    yticklabels = []
+    col = "avg_delivery"
+    for scenario in range(60, 70):
+        color = colors[scenario - 60]
+        date = datetime.date(2000, 4, 15) + datetime.timedelta(
+            days=(scenario - 60) * 5
+        )
+        scendata = df[df["scenario"] == scenario].copy().set_index("valid")
+        delta = scendata[col] - baseline[col]
+        delta = delta[delta != 0]
+        total = (
+            (scendata[col].sum() - baseline[col].sum()) / baseline[col].sum()
+        ) * 100.0
+        yticklabels.append("%s %4.2f%%" % (date.strftime("%b %d"), total))
+        x = delta.index.to_pydatetime()
+        # res = ax.scatter(x, delta.values + (scenario - 60))
+        for idx, val in enumerate(delta):
+            ax.arrow(
+                x[idx],
+                scenario - 60,
+                0,
+                val * 10.0,
+                head_width=4,
+                head_length=0.1,
+                fc=color,
+                ec=color,
+            )
+        ax.axhline(scenario - 60, color=color)
+    ax.set_xlim(datetime.date(year, 1, 1), datetime.date(year + 1, 1, 1))
+    ax.set_ylim(-0.5, 10)
+    ax.xaxis.set_major_locator(mdates.DayLocator([1]))
+    ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
+    ax.set_title(
+        "huc12: %s \n%s Daily Change in Delivery vs Apr 10 Planting"
+        % (huc12, year)
+    )
+    ax.grid(axis="x")
+    ax.set_yticks(range(10))
+    ax.set_yticklabels(yticklabels)
+    plt.gcf().savefig("test.png")
+
+
+if __name__ == "__main__":
+    main(sys.argv)

scripts/tillage_timing/plot_precip.py

@@ -0,0 +1,42 @@
+"""Plot of 5 day averaged precipitation."""
+import calendar
+
+from pyiem.plot.use_agg import plt
+from pyiem.util import get_dbconn
+from pandas.io.sql import read_sql
+
+
+def main():
+    """Go Main Go."""
+    myhucs = [x.strip() for x in open("myhucs.txt").readlines()]
+    pgconn = get_dbconn("idep")
+    df = read_sql(
+        """
+        SELECT extract(year from valid),
+        (extract(doy from valid)::int / 5)::int as datum,
+        max_precip
+        from results_by_huc12 where scenario = 0 and
+        valid < '2019-01-01' and valid >= '2008-01-01'
+        and huc_12 in %s and max_precip >= 50
+    """,
+        pgconn,
+        params=(tuple(myhucs),),
+    )
+    # each datum should have 30 hucs * 11 years * 5 days of data
+    gdf = df.groupby("datum").count() / (30.0 * 11.0 * 5.0)
+
+    (fig, ax) = plt.subplots(1, 1)
+    ax.bar(gdf.index.values * 5, gdf["max_precip"].values * 100.0, width=5)
+
+    ax.set_xticks([1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335])
+    ax.set_xticklabels(calendar.month_abbr[1:])
+    ax.grid(True)
+    ax.set_title("2008-2018 HUC12 Frequency of Daily Precipitation >= 50mm")
+    ax.set_ylabel("Frequency [%]")
+    ax.set_xlim(-1, 367)
+
+    fig.savefig("test.png")
+
+
+if __name__ == "__main__":
+    main()

scripts/tillage_timing/plot_precip_partitions.py

@@ -0,0 +1,73 @@
+"""Plot of 5 day averaged precipitation."""
+import datetime
+
+from pyiem.plot.use_agg import plt
+from pyiem.util import get_dbconn
+import pandas as pd
+from pandas.io.sql import read_sql
+
+
+def main():
+    """Go Main Go."""
+    myhucs = [x.strip() for x in open("myhucs.txt").readlines()]
+    pgconn = get_dbconn("idep")
+    df = read_sql(
+        """
+        SELECT extract(year from valid),
+        extract(doy from valid)::int as doy,
+        qc_precip
+        from results_by_huc12 where scenario = 0 and
+        valid < '2019-01-01' and valid >= '2008-01-01'
+        and huc_12 in %s
+    """,
+        pgconn,
+        params=(tuple(myhucs),),
+    )
+
+    (fig, ax) = plt.subplots(1, 1)
+    # How many to average by
+    total = 30.0 * 11.0
+    # Create two week windows before and inclusive after a given date
+    top = []
+    bottom = []
+    idxs = []
+    for date in pd.date_range("2000/4/1", "2000/5/31"):
+        idx = int(date.strftime("%j"))
+        idxs.append(idx)
+        sdoy = int((date - datetime.timedelta(days=20)).strftime("%j"))
+        edoy = int((date + datetime.timedelta(days=20)).strftime("%j"))
+        before = df[(df["doy"] >= sdoy) & (df["doy"] < idx)].sum() / total
+        bottom.append(0 - before["qc_precip"])
+        after = df[(df["doy"] >= idx) & (df["doy"] < edoy)].sum() / total
+        top.append(after["qc_precip"])
+
+    ax.bar(idxs, top, label="After Date")
+    ax.bar(idxs, bottom, label="Before Date")
+    ax.legend(loc=2, ncol=2)
+    xticks = []
+    xticklabels = []
+    for date in pd.date_range("2000/4/5", "2000/5/30", freq="5d"):
+        xticks.append(int(date.strftime("%j")))
+        xticklabels.append(date.strftime("%b\n%-d"))
+    ax.set_xticks(xticks)
+    ax.set_xticklabels(xticklabels)
+    ax.set_xlim(xticks[0] - 6, xticks[-1] + 3)
+    # Force labels to update
+    ax.set_ylim(-105, 105)
+    ax.set_yticklabels(ax.get_yticks())
+    labels = []
+    for y in ax.get_yticklabels():
+        val = float(y.get_text())
+        if val < 0:
+            val = 0 - val
+        labels.append("%.0f" % (val,))
+    ax.set_yticklabels(labels)
+    ax.grid(True)
+    ax.set_title("2008-2018 HUC12 Precipitation for 14 Day Periods")
+    ax.set_ylabel("Precipitation per Year [mm]")
+
+    fig.savefig("test.png")
+
+
+if __name__ == "__main__":
+    main()