diff --git a/scripts/tillage_timing/README.md b/scripts/tillage_timing/README.md index a30517fd..f64085f7 100644 --- a/scripts/tillage_timing/README.md +++ b/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. diff --git a/scripts/tillage_timing/plot_dailydeltas.py b/scripts/tillage_timing/plot_dailydeltas.py new file mode 100644 index 00000000..48ad1417 --- /dev/null +++ b/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) diff --git a/scripts/tillage_timing/plot_dailydeltas_huc12.py b/scripts/tillage_timing/plot_dailydeltas_huc12.py new file mode 100644 index 00000000..72cc19ad --- /dev/null +++ b/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) diff --git a/scripts/tillage_timing/plot_precip.py b/scripts/tillage_timing/plot_precip.py new file mode 100644 index 00000000..7efa38d6 --- /dev/null +++ b/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() diff --git a/scripts/tillage_timing/plot_precip_partitions.py b/scripts/tillage_timing/plot_precip_partitions.py new file mode 100644 index 00000000..c8e434f7 --- /dev/null +++ b/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()