From 3b265e5a229fe620b1538cc192ef346cbf4c91dd Mon Sep 17 00:00:00 2001 From: akrherz Date: Mon, 7 Oct 2019 15:36:36 -0500 Subject: [PATCH] more tillage timing work around #67 --- scripts/plots/huc12_flowpath_balance.py | 64 ++++++++++ scripts/tillage_timing/README.md | 19 +++ .../tillage_timing/dynamic_tillage_mod_rot.py | 109 ++++++++++++++++++ scripts/tillage_timing/plot_tillage_dates.py | 40 +++++++ 4 files changed, 232 insertions(+) create mode 100644 scripts/plots/huc12_flowpath_balance.py create mode 100644 scripts/tillage_timing/dynamic_tillage_mod_rot.py create mode 100644 scripts/tillage_timing/plot_tillage_dates.py diff --git a/scripts/plots/huc12_flowpath_balance.py b/scripts/plots/huc12_flowpath_balance.py new file mode 100644 index 00000000..214722cd --- /dev/null +++ b/scripts/plots/huc12_flowpath_balance.py @@ -0,0 +1,64 @@ +"""Diagnostic on HUC12 flowpath balance.""" + +import numpy as np +import cartopy.crs as ccrs +from matplotlib.patches import Polygon +import matplotlib.colors as mpcolors +from geopandas import read_postgis +from pyiem.util import get_dbconn +from pyiem.plot.use_agg import plt +from pyiem.plot.geoplot import MapPlot + + +def main(): + """Go Main Go.""" + pgconn = get_dbconn('idep') + df = read_postgis(""" + with centroids as ( + select huc_12, st_centroid(geom) as center, simple_geom from huc12 + where scenario = 0), + agg as ( + select c.huc_12, + sum(case when st_y(center) < st_ymax(geom) then 1 else 0 end) as west, + count(*) from flowpaths f JOIN centroids c on + (f.huc_12 = c.huc_12) WHERE f.scenario = 0 + GROUP by c.huc_12) + select a.huc_12, st_transform(c.simple_geom, 4326) as geo, + a.west, a.count from agg a JOIN centroids c + ON (a.huc_12 = c.huc_12) + """, pgconn, index_col=None, geom_col='geo') + df['percent'] = df['west'] / df['count'] * 100. + bins = np.arange(0, 101, 10) + cmap = plt.get_cmap('RdBu') + norm = mpcolors.BoundaryNorm(bins, cmap.N) + mp = MapPlot( + continentalcolor='thistle', nologo=True, + sector='custom', + south=36.8, north=48.0, west=-99.2, east=-88.9, + subtitle='', + title=('DEP Flowpaths North of HUC12 Centroid (%.0f/%.0f %.2f%%)' % ( + df['west'].sum(), df['count'].sum(), + df['west'].sum() / df['count'].sum() * 100. + ))) + for _i, row in df.iterrows(): + c = cmap(norm([row['percent'], ]))[0] + arr = np.asarray(row['geo'].exterior) + points = mp.ax.projection.transform_points( + ccrs.Geodetic(), arr[:, 0], arr[:, 1]) + p = Polygon(points[:, :2], fc=c, ec='None', zorder=2, lw=0.1) + mp.ax.add_patch(p) + mp.drawcounties() + mp.draw_colorbar( + bins, cmap, norm, + title='Percent', extend='neither') + mp.postprocess(filename='/tmp/huc12_north.png') + + # gdf = df.groupby('fps').count() + # gdf.columns = ['count', ] + # gdf['cumsum'] = gdf['count'].cumsum() + # gdf['percent'] = gdf['cumsum'] / gdf['count'].sum() * 100. + # gdf.to_csv('/tmp/huc12_flowpath_cnts.csv') + + +if __name__ == '__main__': + main() diff --git a/scripts/tillage_timing/README.md b/scripts/tillage_timing/README.md index a2d95557..623fca5a 100644 --- a/scripts/tillage_timing/README.md +++ b/scripts/tillage_timing/README.md @@ -23,3 +23,22 @@ 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 | +75 | May 5 Tillage, Local Climate | +76 | May 10 Tillage, Local Climate | +77 | May 15 Tillage, Local Climate | +78 | May 20 Tillage, Local Climate | +79 | May 25 Tillage, Local Climate | +80 | May 30 Tillage, Local Climate | +81 | Dynamic Tillage April 15 after 45% | +82 | Dynamic Tillage April 15 after 40% | +83 | Dynamic Tillage April 15 after 35% | +84 | Dynamic Tillage April 15 after 30% | +85 | Dynamic Tillage April 15 after 25% | +86 | Dynamic Tillage April 15 after 20% | +87 | Dynamic Tillage April 15 after 15% | +88 | Dynamic Tillage April 15 after Plastic Limit | diff --git a/scripts/tillage_timing/dynamic_tillage_mod_rot.py b/scripts/tillage_timing/dynamic_tillage_mod_rot.py new file mode 100644 index 00000000..9b2f3bc3 --- /dev/null +++ b/scripts/tillage_timing/dynamic_tillage_mod_rot.py @@ -0,0 +1,109 @@ +"""Yikes, inspect WB file, do dynamic tillage dates for 2018.""" +import sys +import datetime + +import pandas as pd +from pandas.io.sql import read_sql +from pyiem.util import get_dbconn +from pyiem.dep import read_wb +from tqdm import tqdm + +APR15 = pd.Timestamp(year=2018, month=4, day=15) +MAY30 = pd.Timestamp(year=2018, month=5, day=30) +THRESHOLDS = { + 81: 45, + 82: 40, + 83: 35, + 84: 30, + 85: 25, + 86: 20, + 87: 15, +} + + +def date_diagnostic(scenario, dates): + """Look at some stats.""" + running = 0 + with open("scenario%s_dates.txt" % (scenario, ), 'w') as fp: + fp.write("date,total\n") + for date in pd.date_range(APR15, MAY30): + hits = [d for d in dates if d == date] + running += len(hits) + fp.write("%s,%s\n" % (date, running)) + + +def main(argv): + """Go Main Go.""" + scenario = int(argv[1]) + dbconn = get_dbconn('idep') + # 0. Loop over flowpaths + df = read_sql(""" + SELECT huc_12, fpath from flowpaths where scenario = %s + """, dbconn, params=(scenario, )) + print("Found %s flowpaths" % (len(df.index), )) + # dates = {} + # for s in range(81, 88): + # dates[s] = [] + for i, row in tqdm(df.iterrows(), total=len(df.index)): + # 1. Load up scenario 0 WB file + wbfn = "/i/0/wb/%s/%s/%s_%s.wb" % ( + row["huc_12"][:8], row["huc_12"][8:], row['huc_12'], row['fpath']) + wbdf = read_wb(wbfn) + wbdf2 = wbdf[( + (wbdf['ofe'] == 1) & (wbdf['date'] >= APR15) & + (wbdf['date'] <= MAY30))] + wbdf3 = wbdf2[wbdf2['sw1'] < THRESHOLDS[scenario]] + if len(wbdf3.index) > 0: + tillage_date = wbdf3.iloc[0]['date'] + else: + tillage_date = MAY30 + # 2. load up prj file, figuring out which .rot file to edit + prjfn = "/i/59/prj/%s/%s/%s_%s.prj" % ( + row["huc_12"][:8], row["huc_12"][8:], row['huc_12'], row['fpath']) + newprjfn = "/i/%s/%s" % (scenario, prjfn[5:]) + with open(newprjfn, 'w') as fp: + i = 0 + for line in open(prjfn): + if line.find(".rot") == -1: + fp.write(line) + continue + if not line.strip().startswith("File = "): + fp.write(line) + continue + rotfn = ( + "/home/akrherz/projects/dep/prj2wepp/wepp/data/" + "managements/SCEN59/%s" + ) % (line.split("SCEN59/")[1][:-2], ) + newrotfn = "/i/%s/custom_rot/%s_%s_%s.rot" % ( + scenario, row['huc_12'], row['fpath'], i + ) + # 5. mod .prj file to see this new .rot + fp.write(" File = \"%s\"\n" % (newrotfn, )) + i += 1 + # 3. Modify 2018 (year 12) .rot dates + # 4. Save out .rot file somewhere per-flowpath specific + with open(newrotfn, 'w') as fp2: + for line2 in open(rotfn): + tokens = line2.split() + if (len(tokens) > 4 and tokens[2] == "12" + and tokens[4] in ["Tillage", "Plant-Annual"] + and int(tokens[0]) < 7): + date = datetime.date( + 2018, int(tokens[0]), int(tokens[1]) + ) + offset = datetime.date(2018, 4, 10) - date + nd = tillage_date - offset + newline = "%02i %02i %s" % ( + nd.month, nd.day, " ".join(tokens[2:])) + fp2.write(newline+"\n") + else: + fp2.write(line2) + + # 6. cry + # for s in THRESHOLDS: + # date_diagnostic(s, dates[s]) + + + +if __name__ == '__main__': + main(sys.argv) diff --git a/scripts/tillage_timing/plot_tillage_dates.py b/scripts/tillage_timing/plot_tillage_dates.py new file mode 100644 index 00000000..6759c56d --- /dev/null +++ b/scripts/tillage_timing/plot_tillage_dates.py @@ -0,0 +1,40 @@ +"""Generate a plot.""" + +import pandas as pd +import matplotlib.dates as mdates +from pyiem.plot.use_agg import plt +THRESHOLDS = { + 81: 45, + 82: 40, + 83: 35, + 84: 30, + 85: 25, + 86: 20, + 87: 15, +} + + +def main(): + """Go Main Go.""" + (fig, ax) = plt.subplots(1, 1) + for scenario in THRESHOLDS: + df = pd.read_csv( + "scenario%s_dates.txt" % (scenario, ), header=None, + names=('date', 'total')) + df['date'] = pd.to_datetime(df['date']) + ax.plot( + df['date'].values, df['total'].values / df['total'].max() * 100., + label='%.0f%%' % (THRESHOLDS[scenario], )) + ax.grid(True) + ax.xaxis.set_major_formatter(mdates.DateFormatter("%-d %b")) + ax.set_xlabel("Date of 2018, tillage done on 30 May if threshold unmeet.") + ax.set_ylabel("Percentage of Flowpaths Tilled") + ax.set_title("DEP Tillage Timing based on 0-10cm VWC") + ax.legend(loc=4, ncol=3) + ax.set_yticks([0, 5, 10, 25, 50, 75, 90, 95, 100]) + ax.set_ylim(0, 100.1) + fig.savefig('test.png') + + +if __name__ == '__main__': + main()