diff --git a/scripts/dynamic_tillage/assemble_csvs.py b/scripts/dynamic_tillage/assemble_csvs.py index ce84c0d9..e1b1a5b4 100644 --- a/scripts/dynamic_tillage/assemble_csvs.py +++ b/scripts/dynamic_tillage/assemble_csvs.py @@ -8,19 +8,12 @@ def main(): """Go main Go.""" dfs = [] - for csvfn in glob.glob("*_??.csv"): + for csvfn in glob.glob("plotsv2/*_??.csv"): crop, _year, state = csvfn[:-4].split("_") if state not in ["MN", "IA"]: continue progress = pd.read_csv(csvfn, parse_dates=["valid"]) progress["state"] = state - # Need to fill out NASS with 100s - firstidx = progress[f"{crop} planted"].first_valid_index() - if firstidx > 0: - progress.loc[:firstidx, f"{crop} planted"] = 0 - lastidx = progress[f"{crop} planted"].last_valid_index() - if lastidx is not None: - progress.loc[(lastidx + 1) :, f"{crop} planted"] = 100 progress = progress.rename( columns={"dep_planted": f"{crop}_dep_planted"} ).set_index(["state", "valid"]) @@ -29,24 +22,17 @@ def main(): rectified = jumbo[ jumbo["corn planted"].notna() | jumbo["dep_corn_planted"].notna() ].copy() - soybeans = jumbo[ - jumbo["soybeans planted"].notna() - | jumbo["dep_soybeans_planted"].notna() - ].copy() - rectified["soybeans planted"] = soybeans["soybeans planted"] - rectified["dep_soybeans_planted"] = soybeans["dep_soybeans_planted"] # Drop things with all nulls rectified = rectified[ rectified["dep_corn_planted"].notna() - | rectified["dep_soybeans_planted"].notna() | rectified["dep_days_suitable"].notna() ] ( rectified.reset_index() .sort_values(["state", "valid"]) - .to_csv("IA_MN_dep_vs_nass_240715.csv", index=False) + .to_csv("IA_MN_dep_vs_nass_240828.csv", index=False) ) diff --git a/scripts/dynamic_tillage/fourpanel_limiter.py b/scripts/dynamic_tillage/fourpanel_limiter.py index 905d5ceb..4aae46c9 100644 --- a/scripts/dynamic_tillage/fourpanel_limiter.py +++ b/scripts/dynamic_tillage/fourpanel_limiter.py @@ -10,7 +10,9 @@ @click.command() -@click.option("--date", "dt", type=click.DateTime(), help="Date to plot.") +@click.option( + "--date", "dt", required=True, type=click.DateTime(), help="Date to plot." +) def main(dt): """Go Main Go.""" with get_sqlalchemy_conn("idep") as conn: @@ -61,7 +63,7 @@ def main(dt): states.plot(ax=ax, color="None", edgecolor="k", lw=0.5, aspect=None) ax.set_xlim(xlim) ax.set_ylim(ylim) - fig.savefig(f"huc12limiter_{dt:%Y%m%d}.png") + fig.savefig(f"plotsv2/huc12limiter_{dt:%Y%m%d}.png") if __name__ == "__main__": diff --git a/scripts/dynamic_tillage/huc12_timeseries.py b/scripts/dynamic_tillage/huc12_timeseries.py index ac257fd0..282d39d1 100644 --- a/scripts/dynamic_tillage/huc12_timeseries.py +++ b/scripts/dynamic_tillage/huc12_timeseries.py @@ -10,15 +10,21 @@ from sqlalchemy import text -def get_plastic_limit(huc12: str) -> pd.DataFrame: +def get_plastic_limit(huc12: str, year: int) -> pd.DataFrame: """Figure out what the plastic limit is.""" + charat = year - 2007 with get_sqlalchemy_conn("idep") as conn: # build up the cross reference of everyhing we need to know return pd.read_sql( text( """ - select o.ofe, p.fpath, o.fbndid, g.plastic_limit, - p.fpath || '_' || o.ofe as combo + select o.ofe, p.fpath, o.fbndid, + coalesce( + (g.wepp_max_sw + g.wepp_min_sw) / 2.0, + g.plastic_limit + ) as plastic_limit, + p.fpath || '_' || o.ofe as combo, + substr(landuse, :charat, 1) as crop from flowpaths p, flowpath_ofes o, gssurgo g WHERE o.flowpath = p.fid and p.huc_12 = :huc12 and p.scenario = 0 and o.gssurgo_id = g.id @@ -27,20 +33,21 @@ def get_plastic_limit(huc12: str) -> pd.DataFrame: conn, params={ "huc12": huc12, + "charat": charat, }, index_col=None, ) @click.command() -@click.option("--huc12", help="HUC12 To Plot") -@click.option("--year", type=int, help="Year to Plot") -def main(huc12, year): +@click.option("--huc12", required=True, help="HUC12 To Plot") +@click.option("--year", required=True, type=int, help="Year to Plot") +def main(huc12: str, year: int): """Go Main Go.""" - pldf = get_plastic_limit(huc12) + pldf = get_plastic_limit(huc12, year) smdfs = [] hucstatusdfs = [] - for dt in pd.date_range(f"{year}/04/14", f"{year}/05/31"): + for dt in pd.date_range(f"{year}/04/11", f"{year}/06/13"): fn = f"/mnt/idep2/data/smstate/{dt:%Y}/smstate{dt:%Y%m%d}.feather" if not os.path.isfile(fn): continue @@ -54,24 +61,35 @@ def main(huc12, year): statusdf["date"] = dt hucstatusdfs.append(statusdf.loc[[huc12]]) hucstatusdf = pd.concat(hucstatusdfs) - print(hucstatusdf) huc12sm = pd.concat(smdfs) huc12sm["combo"] = ( huc12sm["fpath"].astype(str) + "_" + huc12sm["ofe"].astype(str) ) - huc12sm["plastic_limit"] = ( + huc12sm["pl0.8"] = ( huc12sm["combo"].map(pldf.set_index("combo")["plastic_limit"]) * 0.8 ) - huc12sm["below_pl"] = huc12sm["sw1"] < huc12sm["plastic_limit"] + huc12sm["crop"] = huc12sm["combo"].map(pldf.set_index("combo")["crop"]) + huc12sm["below_pl"] = huc12sm["sw1"] < huc12sm["pl0.8"] fig = figure( - title=f"DEP: {huc12} Walnut Creek Dynamic Tillage Diagnostic", - subtitle="15 April - 23 April 2024", + title=f"DEP: {huc12} Brush Creek Dynamic Tillage Diagnostic", + subtitle=f"11 April - 15 June {year}", logo="dep", - figsize=(8, 6), + figsize=(8, 8), ) - yaxsize = 0.2 - ax = fig.add_axes([0.1, 0.68, 0.8, yaxsize]) + yaxsize = 0.15 + # Scatter plot of Soil Moisture + ax0 = fig.add_axes([0.1, 0.75, 0.8, yaxsize]) + past = huc12sm[huc12sm["crop"] == "P"] + ax0.scatter(past["date"], past["sw1"] - past["pl0.8"], c="r", s=2) + nonp = huc12sm[huc12sm["crop"] != "P"] + ax0.scatter(nonp["date"], nonp["sw1"] - nonp["pl0.8"], c="b", s=2) + ax0.xaxis.set_major_formatter(DateFormatter("%-d %b")) + ax0.xaxis.set_major_formatter(DateFormatter("%-d %b")) + ax0.grid(True) + + # Percent of OFEs below 0.8*PL + ax = fig.add_axes([0.1, 0.51, 0.8, yaxsize], sharex=ax0) ax.text( 0, 1, "* SM Evaluation using previous date", transform=ax.transAxes ) @@ -89,9 +107,8 @@ def main(huc12, year): ax.set_yticks([0, 10, 25, 50, 75, 90, 100]) ax.grid(True) ax.set_ylabel("Percent of OFEs\nBelow 0.8*PL") - ax.xaxis.set_major_formatter(DateFormatter("%-d %b")) - ax2 = fig.add_axes([0.1, 0.4, 0.8, yaxsize], sharex=ax) + ax2 = fig.add_axes([0.1, 0.3, 0.8, yaxsize], sharex=ax0) ax2.bar( hucstatusdf["date"], hucstatusdf["precip_mm"], @@ -103,7 +120,7 @@ def main(huc12, year): ax2.axhline(10, lw=2, color="k") ax2.grid(True) - ax3 = fig.add_axes([0.1, 0.1, 0.8, yaxsize], sharex=ax) + ax3 = fig.add_axes([0.1, 0.05, 0.8, yaxsize], sharex=ax0) ax3.bar( hucstatusdf["date"], hucstatusdf["tsoil_avg"], @@ -130,7 +147,7 @@ def main(huc12, year): color="red" if row[f"limited_by_{label}"] else "green", ) - fig.savefig(f"sm_{huc12}_{year}.png") + fig.savefig(f"plotsv2/sm_{huc12}_{year}.png") if __name__ == "__main__": diff --git a/scripts/dynamic_tillage/nass_comparison.py b/scripts/dynamic_tillage/nass_comparison.py index 380e6909..d37363e2 100644 --- a/scripts/dynamic_tillage/nass_comparison.py +++ b/scripts/dynamic_tillage/nass_comparison.py @@ -5,6 +5,7 @@ import click import geopandas as gpd +import matplotlib.colors as mpcolors import numpy as np import pandas as pd from matplotlib.patches import Rectangle @@ -123,6 +124,16 @@ def get_nass(year, state, district, crop) -> pd.DataFrame: ) nass["valid"] = pd.to_datetime(nass["valid"]) nass = nass.pivot(index="valid", columns="metric", values="num_value") + # NASS does not always include 0% and 100% in their data, so we need to + # fill in the blanks + col = f"{crop} planted" + if nass[col].min() > 0: + firstrow = nass[nass[col] > 0].index[0] + nass.loc[firstrow - pd.Timedelta(days=7), col] = 0 + # TODO, this may not always be right for in-progress years + if nass[col].max() < 100: + lastrow = nass[nass[col] > 0].index[-1] + nass.loc[lastrow + pd.Timedelta(days=7), col] = 100 # Create space for DEP to store things nass[f"dep_{crop}_planted"] = np.nan nass["dep_days_suitable"] = np.nan @@ -155,6 +166,17 @@ def get_fields(year, district, state, crop) -> gpd.GeoDataFrame: return fields +def get_labels(year): + """Figure out where we should place xticks.""" + xticks = [] + xticklabels = [] + for dt in pd.date_range(f"{year}-04-11", f"{year}-06-20"): + if dt.dayofweek == 6: + xticks.append(dt) + xticklabels.append(f"{dt:%b %-d}") + return xticks, xticklabels + + @click.command() @click.option("--year", type=int, help="Year to plot") @click.option("--district", type=str, help="NASS District (lowercase)") @@ -170,6 +192,8 @@ def main(year, district, state, crop): nass = get_nass(year, state, district, crop) fields = get_fields(year, district, state, crop) + xticks, xticklabels = get_labels(year) + # Spatially filter fields that are inside the climdiv region daily_limits = compute_limits(fields["huc12"].unique(), year) @@ -210,6 +234,7 @@ def main(year, district, state, crop): ] i += 1 cmap = get_cmap("viridis") + norm = mpcolors.BoundaryNorm(np.arange(0, 101, 10), cmap.N) res = ax.imshow( data, extent=[ @@ -220,14 +245,15 @@ def main(year, district, state, crop): ], aspect="auto", interpolation="nearest", - vmin=0, - vmax=100, + norm=norm, cmap=cmap, ) ax.set_yticks(range(4)) ax.set_yticklabels( ["Combined", "Soil Moisture", "Soil Temp", "Precipitation"] ) + ax.set_xticks(xticks) + ax.set_xticklabels(xticklabels) ax.text(0.02, 1.05, "Limiting Factors", transform=ax.transAxes) cax = fig.add_axes([0.92, 0.05, 0.02, 0.2]) fig.colorbar(res, cax=cax, label="Percent of HUC12s Limited") @@ -284,7 +310,8 @@ def main(year, district, state, crop): # //////////////////////////////////////////////////////////// ax2 = fig.add_axes([0.1, 0.4, 0.8, 0.5]) - ax2.set_ylim(0, 100) + ax2.set_ylim(-2, 102) + ax2.set_yticks(range(0, 101, 10)) ax2.set_ylabel("Percent of Acres Planted") ax2.grid(True) ax2.plot( @@ -303,36 +330,40 @@ def main(year, district, state, crop): ax2.legend(loc=2) # create a table in the lower right corner with the values from - txt = ["Date | NASS | DEP"] + txt = ["Weekly Progress", "Date NASS DEP"] for valid, row in nass.iterrows(): if valid not in accum.index: continue dep = accum.at[valid, "percent"] nass.at[valid, f"dep_{crop}_planted"] = dep val = row[f"{crop} planted"] - txt.append(f"{valid:%Y-%m-%d} | {val:.0f} | {dep:.1f}") + txt.append(f"{valid:%b %-2d}: {val:3.0f} {dep:3.0f}") ax2.text( - 0.99, - 0.01, + 1.05, + 0.05, "\n".join(txt), transform=ax2.transAxes, fontsize=10, + fontfamily="monospace", va="bottom", ha="right", bbox=dict(facecolor="white", edgecolor="black", pad=5), ) + ax2.set_xticks(xticks) + ax2.set_xticklabels(xticklabels) # Sync up the two xaxes ax.set_xlim(*ax2.get_xlim()) ax3.set_xlim(*ax2.get_xlim()) nass.to_csv( - f"{crop}_{year}_{district if district is not None else state}.csv" + f"plotsv2/{crop}_{year}_" + f"{district if district is not None else state}.csv" ) ss = f"state_{state}" fig.savefig( - f"{crop}_progress_{year}_" + f"plotsv2/{crop}_progress_{year}_" f"{district if district is not None else ss}.png" ) diff --git a/src/pydep/tillage.py b/src/pydep/tillage.py index 9f2dcc69..fc3f9ab8 100644 --- a/src/pydep/tillage.py +++ b/src/pydep/tillage.py @@ -54,7 +54,7 @@ } -def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year): +def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year: int): """Read a block file and do replacements Args: @@ -76,9 +76,11 @@ def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year): with open(blockfn, "r", encoding="utf8") as fh: data = fh.read() # Special consideration for planting alfalfa - if code == "P" and prevcode != "P": - # Best we can do now is plant it on Apr 15, sigh + if code == "P" and (prevcode != "P" or year == 1): + # Best we can do now is plant it on Apr 15 with a simple disk the + # day before data = ( + f"4 14 {year} 1 Tillage OpCropDef.FCSTACDP {{0.101600, 2}}\n" f"4 15 {year} 1 Plant-Perennial CropDef.ALFALFA {{0.000000}}\n" f"{data}" ) diff --git a/src/regression_tests/CB_5percent_5till/answer.json b/src/regression_tests/CB_5percent_5till/answer.json index 57e1e4c6..de4fc7a6 100644 --- a/src/regression_tests/CB_5percent_5till/answer.json +++ b/src/regression_tests/CB_5percent_5till/answer.json @@ -1,4 +1,5 @@ { "av_det_2012": 4.499, - "av_det": 4.264 + "av_det": 4.264, + "max_sw1": 59.61 } diff --git a/src/regression_tests/pasture/101800140707_651.man b/src/regression_tests/pasture/101800140707_651.man new file mode 100644 index 00000000..395ca100 --- /dev/null +++ b/src/regression_tests/pasture/101800140707_651.man @@ -0,0 +1,237 @@ +98.4 +# +# +# +# + +1 # number of OFE's +18 # (total) years in simulation + +####################### +# Plant Section # +####################### + +2 # Number of plant scenarios + + +bromegr1 +`Bromegrass-High Fertilization Level' +(from WEPP distribution database) + +1 #landuse +WeppWillSet +14.00000 23.00000 35.00000 10.00000 5.00000 30.00000 0.10000 0.15200 0.70000 0.00220 +0.85000 0.90000 0.65000 0.99000 12.00000 0.00000 0.90000 0.51000 +2 # mfo - +0.00900 0.00900 25.00000 0.00000 0.00600 0.30000 0.33000 0.34000 14 32.00000 +1.10000 9.00000 0.00000 + +ALFALFA +Alfalfa +J. M. Laflen +Set minimum temperature to 0.5 degrees C - dcf +1 #landuse +WeppWillSet +14.00000 23.00000 13.00159 4.00000 4.99968 30.00000 0.10000 0.14999 0.90000 0.00450 +0.85000 0.90000 0.65000 0.99000 12.00000 0.00000 0.90000 0.80011 +1 # mfo - +0.01500 0.01500 20.00000 0.25000 0.00600 2.40008 0.33000 0.60002 14 32.00000 +0.50000 6.00000 0.00000 + +####################### +# Operation Section # +####################### + +1 # Number of operation scenarios + + +FCSTACDP +`Field cultivator, secondary tillage, after duckfoot points +(from WEPP distribution database) +Maximum depth of 10 cm (4 inches) +1 #landuse +0.6000 0.3500 0 +4 # pcode - other +0.0250 0.3000 0.6000 0.3500 0.0150 1.0000 0.0500 + + + +############################### +# Initial Conditions Section # +############################### + +1 # Number of initial scenarios + + +gra_3425 +Initial conditions for a perennial grass strip already in existence +Can use this initial condition to run permanent grass strips +and have no disturbance, planting or other operations +1 #landuse +1.10000 0.50000 200 92 0.00000 0.50000 +1 # iresd +2 # mang perennial +500.00000 0.02000 0.50000 0.02000 0.00000 +1 # rtyp - temporary +0.00000 0.00000 0.10000 0.20000 0.00000 +0.20000 0.00000 + + + + +############################ +# Surface Effects Section # +############################ + +1 # Number of Surface Effects Scenarios + + +# +# Surface Effects Scenario 1 of 1 +# +Year 1 +From WEPP database +Your name, phone + +1 # landuse - cropland +1 # ntill - number of operations + 104 # mdate --- 4 / 14 + 1 # op --- FCSTACDP + 0.102 # depth + 2 # type + + +####################### +# Contouring Section # +####################### + +0 # Number of contour scenarios + +####################### +# Drainage Section # +####################### + +0 # Number of drainage scenarios + + +####################### +# Yearly Section # +####################### + +2 # looper; number of Yearly Scenarios +# +# Yearly scenario 1 of 2 +# +Year 1 + + + +1 # landuse +2 # plant growth scenario +1 # surface effect scenario +0 # contour scenario +0 # drainage scenario +2 # management + 0 # senescence date + 105 # perennial plant date --- 4 /15 + 0 # perennial stop growth date --- 0/0 + 0.0000 # row width + 1 # crop management - + 3 # number of cuttings + 152 # cutting date --- 6/1 + 196 # cutting date --- 7/15 + 244 # cutting date --- 9/1 +# +# Yearly scenario 2 of 2 +# +Year 2 + + + +1 # landuse +2 # plant growth scenario +0 # surface effect scenario +0 # contour scenario +0 # drainage scenario +2 # management + 0 # senescence date + 0 # perennial plant date --- 0 /0 + 0 # perennial stop growth date --- 0/0 + 0.0000 # row width + 1 # crop management - + 3 # number of cuttings + 152 # cutting date --- 6/1 + 196 # cutting date --- 7/15 + 244 # cutting date --- 9/1 + + +####################### +# Management Section # +####################### + +Manage +description 1 +description 2 +description 3 +1 # number of OFE's + 1 # initial condition index +1 # rotation repeats +18 # years in rotation + +# +# Rotation 1: year 1 to 18 +# + + 1 # - OFE: 1> + 1 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index +#----------------------------------- + 1 # - OFE: 1> + 2 # year index diff --git a/src/regression_tests/pasture/101800140707_651.slp b/src/regression_tests/pasture/101800140707_651.slp new file mode 100644 index 00000000..4ec93604 --- /dev/null +++ b/src/regression_tests/pasture/101800140707_651.slp @@ -0,0 +1,10 @@ +97.5 +# +# from flowpath2prj.py 2024-08-29 09:36:08.056228 +# +# +1 +111.9581 1.000 +9 27.7200 +0.000000,0.292145 0.164576,0.292145 0.329152,0.473558 0.493729,0.313349 0.583129,0.052918 0.681027,0.062842 0.784248,0.095914 0.890940,0.089303 1.000000,0.069455 + diff --git a/src/regression_tests/pasture/101800140707_651.sol b/src/regression_tests/pasture/101800140707_651.sol new file mode 100644 index 00000000..42afb49a --- /dev/null +++ b/src/regression_tests/pasture/101800140707_651.sol @@ -0,0 +1,8 @@ +2006.2 +comments: soil file +1 1 +'Valentine' 'FS' 3 0.230000 0.750000 5340560.000000 0.015800 2.080000 35.799999 + 130 94.0 3.0 0.750 2.7 0.0 + 300 95.0 3.0 0.500 2.6 0.0 + 2000 96.0 3.0 0.250 2.4 0.0 +0 0.000000 0 diff --git a/src/regression_tests/pasture/answer.json b/src/regression_tests/pasture/answer.json new file mode 100644 index 00000000..c901d680 --- /dev/null +++ b/src/regression_tests/pasture/answer.json @@ -0,0 +1,5 @@ +{ + "av_det_2012": 0.084, + "av_det": 0.084, + "max_sw1": 35.07 +} diff --git a/src/regression_tests/pasture/fp.run b/src/regression_tests/pasture/fp.run new file mode 100644 index 00000000..ce1cdc1f --- /dev/null +++ b/src/regression_tests/pasture/fp.run @@ -0,0 +1,31 @@ +E +Yes +1 +1 +No +1 +No +/dev/null +Yes +wepp_wb.txt +Yes +wepp_res.txt +Yes +wepp_soil.txt +No +Yes +wepp_grph.txt +Yes +wepp_env.txt +No +No +No +Yes +wepp_yld.txt +101800140707_651.man +101800140707_651.slp +../common/092.51x042.00.cli +101800140707_651.sol +0 +15 +0 diff --git a/src/regression_tests/run_regressions.py b/src/regression_tests/run_regressions.py index 27ffc0d1..63729f8b 100644 --- a/src/regression_tests/run_regressions.py +++ b/src/regression_tests/run_regressions.py @@ -5,6 +5,7 @@ import subprocess import sys +from pydep.io.dep import read_wb from pydep.io.wepp import read_env @@ -36,6 +37,17 @@ def main(): ) if failed: sys.exit(3) + + df = read_wb("wepp_wb.txt") + max_sw1 = df["sw1"].max() + failed = abs(max_sw1 - answer["max_sw1"]) > 1 # forgiving + sys.stdout.write( + f"{'FAILED' if failed else 'PASS'}: {fn}, " + f"expected {answer['max_sw1']:.3f}, got {max_sw1:.3f}\n" + ) + if failed: + sys.exit(3) + os.chdir("..") diff --git a/tests/test_tillage.py b/tests/test_tillage.py index 3d19897d..804a063d 100644 --- a/tests/test_tillage.py +++ b/tests/test_tillage.py @@ -3,6 +3,12 @@ from pydep import tillage +def test_gh266_pasture(): + """Test that an initial planting operation is created.""" + res = tillage.make_tillage(0, "IA_NORTH", "P", "P", "P", 1, 1) + assert "Plant-Perennial" in res + + def test_simple(): """Test import of API.""" plants = ["C", "B", "W", "P"]