more tillage scenario work #67

scripts/tillage_timing/README.md

@@ -27,18 +27,20 @@ ID | Name | Status
 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 |
+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
+81 | Dynamic Tillage April 15 after 45% | done
+82 | Dynamic Tillage April 15 after 40% | done
+83 | Dynamic Tillage April 15 after 35% | done
+84 | Dynamic Tillage April 15 after 30% | done
+85 | Dynamic Tillage April 15 after 25% | done
+86 | Dynamic Tillage April 15 after 20% | done
+87 | Dynamic Tillage April 15 after 15% | done
+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

scripts/tillage_timing/dynamic_tillage_mod_rot.py

@@ -19,6 +19,11 @@
     86: 20,
     87: 15,
 }
+PL_THRESHOLDS = {
+    88: 1.,
+    89: 0.9,
+    90: None,
+}
 
 
 def date_diagnostic(scenario, dates):
@@ -32,6 +37,46 @@ def date_diagnostic(scenario, dates):
             fp.write("%s,%s\n" % (date, running))
 
 
+def get_threshold_bypl(scenario, row):
+    """Inspect the soil file for what we need."""
+    fn = "/i/0/sol/%s/%s/%s_%s.sol" % (
+        row["huc_12"][:8], row["huc_12"][8:], row['huc_12'], row['fpath'])
+    # Take the top layer of the first soil, good enough for now.
+    line = open(fn).readlines()[4]
+    tokens = line.strip().split()
+    # depth, sand, clay, OM, CEC, Rock
+    clay = float(tokens[2])
+    om = float(tokens[3])
+    pl = 14.22 + 0.005 * clay**2 + 3.63 * om - 0.048 * clay * om
+    if PL_THRESHOLDS[scenario] is None:
+        line = open(fn).readlines()[3]
+        tokens = line.strip().split()
+        satv = float(tokens[4])
+        satm = satv / 2.65 * (1 - satv)
+        return ((pl / 100. - 0.4 * satm) / 0.6) * 100.
+    return pl * PL_THRESHOLDS[scenario]
+
+
+def compute_tillage_date(scenario, row):
+    """Get the tillage date."""
+    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))]
+    if scenario in THRESHOLDS:
+        the_threshold = THRESHOLDS[scenario]
+    else:
+        the_threshold = get_threshold_bypl(scenario, row)
+    wbdf3 = wbdf2[wbdf2['sw1'] < the_threshold]
+    if len(wbdf3.index) > 0:
+        tillage_date = wbdf3.iloc[0]['date']
+    else:
+        tillage_date = MAY30
+    return tillage_date
+
+
 def main(argv):
     """Go Main Go."""
     scenario = int(argv[1])
@@ -41,22 +86,11 @@ def main(argv):
         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] = []
+    dates = []
     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
+        tillage_date = compute_tillage_date(scenario, row)
+        dates.append(tillage_date)
         # 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'])
@@ -100,9 +134,7 @@ def main(argv):
                             fp2.write(line2)
 
     # 6. cry
-    # for s in THRESHOLDS:
-    #    date_diagnostic(s, dates[s])
-
+    date_diagnostic(scenario, dates)
 
 
 if __name__ == '__main__':

scripts/tillage_timing/plot_tillage_dates.py

@@ -12,25 +12,39 @@
     86: 20,
     87: 15,
 }
+PL_THRESHOLDS = {
+    88: 1.0,
+    89: 0.9,
+    90: None,
+}
 
 
 def main():
     """Go Main Go."""
     (fig, ax) = plt.subplots(1, 1)
+    ax.set_position([0.1, 0.2, 0.85, 0.75])
     for scenario in THRESHOLDS:
-        df = pd.read_csv(
-            "scenario%s_dates.txt" % (scenario, ), header=None,
-            names=('date', 'total'))
+        df = pd.read_csv("scenario%s_dates.txt" % (scenario, ))
         df['date'] = pd.to_datetime(df['date'])
         ax.plot(
             df['date'].values, df['total'].values / df['total'].max() * 100.,
             label='%.0f%%' % (THRESHOLDS[scenario], ))
+    for scenario in PL_THRESHOLDS:
+        df = pd.read_csv("scenario%s_dates.txt" % (scenario, ))
+        df['date'] = pd.to_datetime(df['date'])
+        if scenario != 90:
+            label = 'PL%.1f' % (PL_THRESHOLDS[scenario], )
+        else:
+            label = 'PL(SAT)'
+        ax.plot(
+            df['date'].values, df['total'].values / df['total'].max() * 100.,
+            label=label, linestyle='-.')
     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.legend(loc=(-0.02, -0.26), ncol=5)
     ax.set_yticks([0, 5, 10, 25, 50, 75, 90, 95, 100])
     ax.set_ylim(0, 100.1)
     fig.savefig('test.png')