diff --git a/pyproject.toml b/pyproject.toml index 882e1aea..39f8e7ea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,3 +9,8 @@ lint.select = [ ] target-version = "py39" +# per memory issues microsoft/pylance-release/issues/5739 +[tool.pyright] +autoSearchPaths = false +exclude = ["data/", "prj2wepp/"] +ignore = ["data/", "prj2wepp/"] diff --git a/scripts/dynamic_tillage/accum_acres_plot.py b/scripts/dynamic_tillage/accum_acres_plot.py index 2e179b1b..21c56a8c 100644 --- a/scripts/dynamic_tillage/accum_acres_plot.py +++ b/scripts/dynamic_tillage/accum_acres_plot.py @@ -222,5 +222,5 @@ def main(): if __name__ == "__main__": - for _i, _dt in enumerate(pd.date_range("2023-04-15", "2023-06-04")): + for _i, _dt in enumerate(pd.date_range("2008-04-15", "2008-06-05")): plot_map_progress(_i, _dt) diff --git a/scripts/dynamic_tillage/fourpanel_limiter.py b/scripts/dynamic_tillage/fourpanel_limiter.py new file mode 100644 index 00000000..905d5ceb --- /dev/null +++ b/scripts/dynamic_tillage/fourpanel_limiter.py @@ -0,0 +1,68 @@ +"""A four panel plot diagnostic showing the limiting factor.""" + +import click +import geopandas as gpd +import numpy as np +import pandas as pd +from pyiem.database import get_sqlalchemy_conn +from pyiem.plot import figure +from pyiem.util import load_geodf + + +@click.command() +@click.option("--date", "dt", type=click.DateTime(), help="Date to plot.") +def main(dt): + """Go Main Go.""" + with get_sqlalchemy_conn("idep") as conn: + huc12df = gpd.read_postgis( + "SELECT huc_12, st_transform(simple_geom, 4326) as geom " + "from huc12 where scenario = 0", + conn, + geom_col="geom", + index_col="huc_12", + ) + + huc12data = pd.read_feather( + f"/mnt/idep2/data/huc12status/{dt:%Y}/{dt:%Y%m%d}.feather" + ) + huc12df = huc12df.join(huc12data, how="inner") + + states = load_geodf("us_states", 4326) + + fig = figure( + title=f"{dt:%d %b %Y} Dynamic Tillage Limiting Factor", + subtitle="Red indicates limiting factor", + logo="dep", + figsize=(10.24, 7.68), + ) + opts = [ + [0.03, 0.03, "Plastic Limit", "limited_by_soilmoisture"], + [0.5, 0.03, "Combined Limit", "limited"], + [0.03, 0.48, "Soil Temp Limit", "limited_by_soiltemp"], + [0.5, 0.48, "Rainfall Limit", "limited_by_precip"], + ] + for opt in opts: + ax = fig.add_axes([opt[0], opt[1], 0.42, 0.38]) + ax.set_xticklabels([]) + ax.set_yticklabels([]) + title = ( + f"{opt[2]} {huc12df[opt[3]].sum()}/{len(huc12df)} " + f"{huc12df[opt[3]].sum() / len(huc12df):.1%}" + ) + ax.set_title(title, fontsize=16) + # plot in red when limited is true, else blue + huc12df.plot( + ax=ax, + aspect=None, + color=np.where(huc12df[opt[3]], "r", "b"), + ) + xlim = ax.get_xlim() + ylim = ax.get_ylim() + 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") + + +if __name__ == "__main__": + main() diff --git a/scripts/dynamic_tillage/nass_comparison.py b/scripts/dynamic_tillage/nass_comparison.py new file mode 100644 index 00000000..94431bec --- /dev/null +++ b/scripts/dynamic_tillage/nass_comparison.py @@ -0,0 +1,131 @@ +"""Plot our planting progress vs that of NASS.""" + +import click +import geopandas as gpd +import pandas as pd +from pyiem.database import get_sqlalchemy_conn +from pyiem.plot import figure_axes +from sqlalchemy import text + +XREF = { + "nw": "IAC001", + "nc": "IAC002", + "ne": "IAC003", + "wc": "IAC004", + "c": "IAC005", + "ec": "IAC006", + "sw": "IAC007", + "sc": "IAC008", + "se": "IAC009", +} + + +@click.command() +@click.option("--year", type=int, help="Year to plot") +@click.option("--district", type=str, help="NASS District") +def main(year, district): + """Go Main Go.""" + charidx = year - 2007 + 1 + with get_sqlalchemy_conn("coop") as conn: + climdiv = gpd.read_postgis( + text( + """ + select t.iemid, r.geom from climodat_regions r JOIN stations t + on (r.iemid = t.iemid) where t.id = :sid + """ + ), + conn, + params={"sid": XREF[district]}, + geom_col="geom", + index_col="iemid", + ) + + # Get the NASS data + with get_sqlalchemy_conn("coop") as conn: + nass = pd.read_sql( + text( + """ + select * from nass_iowa where metric = 'corn planted' + and extract(year from valid) = :year ORDER by valid ASc + """ + ), + conn, + params={"year": year}, + ) + + # Figure out the planting dates + with get_sqlalchemy_conn("idep") as conn: + fields = gpd.read_postgis( + text(""" + select st_transform(geom, 4326) as geom, plant, huc12, fbndid, + acres from fields f LEFT JOIN field_operations o + on (f.field_id = o.field_id and o.year = :year) + where scenario = 0 and + substr(landuse, :charidx, 1) = 'C' + """), + conn, + params={"year": year, "charidx": charidx}, + geom_col="geom", + parse_dates="plant", + ) + + # Spatially filter fields that are inside the climdiv region + fields = gpd.sjoin(fields, climdiv, predicate="intersects") + + # accumulate acres planted by date + accum = fields[["plant", "acres"]].groupby("plant").sum().cumsum() + accum["percent"] = accum["acres"] / fields["acres"].sum() * 100.0 + accum = accum.reindex( + pd.date_range(f"{year}-04-15", f"{year}-06-23") + ).ffill() + + fig, ax = figure_axes( + logo="dep", + title=( + f"{year} Corn Planting Progress vs " + f"NASS Iowa Crop District {district.upper()}" + ), + figsize=(10.24, 7.68), + ) + ax.set_ylim(0, 100) + ax.set_ylabel("Percent of Acres Planted") + ax.set_xlabel("Date") + ax.grid(True) + ax.plot( + accum.index, + accum["acres"] / fields["acres"].sum() * 100.0, + label="DEP", + ) + + ax.scatter( + nass["valid"], + nass[district], + marker="*", + s=30, + color="r", + label="NASS", + ) + ax.legend(loc=2) + + # create a table in the lower right corner with the values from + txt = ["Date | NASS | DEP"] + for _, row in nass.iterrows(): + dep = accum.at[pd.Timestamp(row["valid"]), "percent"] + txt.append(f"{row['valid']} | {row['se']:.0f} | {dep:.1f}") + print("\n".join(txt)) + ax.text( + 0.99, + 0.01, + "\n".join(txt), + transform=ax.transAxes, + fontsize=10, + va="bottom", + ha="right", + bbox=dict(facecolor="white", edgecolor="black", pad=5), + ) + + fig.savefig(f"corn_progress_{year}_{district}.png") + + +if __name__ == "__main__": + main() diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py index 35bc3861..4d976a93 100644 --- a/scripts/dynamic_tillage/workflow.py +++ b/scripts/dynamic_tillage/workflow.py @@ -31,6 +31,7 @@ pd.set_option("future.no_silent_downcasting", True) LOG = logger() CPU_COUNT = min([4, cpu_count() / 2]) +HUC12STATUSDIR = "/mnt/idep2/data/huc12status" PROJDIR = "/opt/dep/prj2wepp" EXE = f"{PROJDIR}/prj2wepp" @@ -112,6 +113,7 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: Returns the number of acres planted and tilled. A negative value is a sentinel that the HUC12 failed due to soil moisture constraint. """ + dbcolidx = dt.year - 2007 + 1 with get_sqlalchemy_conn("idep") as conn: # build up the cross reference of everyhing we need to know df = pd.read_sql( @@ -126,7 +128,9 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: SELECT ofe, fpath, f.fbndid, plastic_limit, f.landuse, f.management, f.acres, f.field_id from fields f LEFT JOIN myofes m on (f.fbndid = m.fbndid) - WHERE f.huc12 = :huc12 and f.isag ORDER by fpath, ofe) + WHERE f.huc12 = :huc12 and + substr(landuse, :dbcolidx, 1) = ANY(:crops) + ORDER by fpath, ofe) select a.*, o.till1, o.till2, o.till3, o.plant, o.year, plant >= :dt as plant_needed, coalesce( @@ -138,7 +142,13 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: """ ), conn, - params={"huc12": huc12, "year": dt.year, "dt": dt}, + params={ + "huc12": huc12, + "year": dt.year, + "dt": dt, + "crops": ["B", "C", "L", "W"], + "dbcolidx": dbcolidx, + }, index_col=None, ) if df.empty: @@ -171,15 +181,16 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: # bit generous here, are 50% of **model** fields above plastic_limit? fields["sw1"] = pd.Series(fieldsw) modelled = (fields["sw1"] > 0).sum() - hits = (fields["sw1"] > fields["plastic_limit"]).sum() - if hits > modelled / 2: - LOG.debug( - "HUC12: %s has %s/%s fields above plastic limit", - huc12, - hits, - modelled, - ) - return -1, -1 + if modelled > 1: + hits = (fields["sw1"] > fields["plastic_limit"]).sum() + if hits > modelled / 2: + LOG.debug( + "HUC12: %s has %s/%s fields above plastic limit", + huc12, + hits, + modelled, + ) + return -1, -1 # Compute things for year char_at = (dt.year - 2007) + 1 @@ -244,7 +255,7 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: def estimate_soiltemp(huc12df, dt): """Write mean GFS soil temperature C into huc12df.""" - ppath = f"{dt:%Y/%m/%d}/model/gfs/gfs_{dt:%Y%m%d}00_iemre.nc" + ppath = f"{dt:%Y/%m/%d}/model/gfs/gfs_{dt:%Y%m%d}12_iemre.nc" with archive_fetch(ppath) as ncfn, ncopen(ncfn) as nc: tsoil = np.mean(nc.variables["tsoil"][1:7], axis=0) - 273.15 y = np.digitize(huc12df["lat"].values, nc.variables["lat"][:]) @@ -349,6 +360,7 @@ def main(scenario, dt, huc12): progress.update() if planted < 0: huc12df.at[huc12, "limited_by_soilmoisture"] = True + huc12df.at[huc12, "limited"] = True else: total_planted += planted total_tilled += tilled @@ -365,6 +377,13 @@ def main(scenario, dt, huc12): huc12df["limited_by_soiltemp"].sum(), ) + # save a debug output file + mydir = f"{HUC12STATUSDIR}/{dt:%Y}" + os.makedirs(mydir, exist_ok=True) + if len(huc12df.index) == 1: # Don't save when run for single huc12 + return + huc12df.drop(columns="geom").to_feather(f"{mydir}/{dt:%Y%m%d}.feather") + if __name__ == "__main__": main()