more tillage timing work around #67

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()

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 |

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)

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()