diff --git a/scripts/cligen/daily_clifile_editor.py b/scripts/cligen/daily_clifile_editor.py index 71ec4b46..1f0c400c 100644 --- a/scripts/cligen/daily_clifile_editor.py +++ b/scripts/cligen/daily_clifile_editor.py @@ -393,15 +393,15 @@ def load_precip(data, dt: date, tile_bounds: BOUNDS): now = midnight # Require at least 75% data coverage, if not, we will abort back to legacy - data["quorum"] = ts * 0.75 + quorum = {"cnt": ts * 0.75} def _cb_lp(args): """callback""" tidx, grid = args # If tidx is outside of bounds, just make this a no-op if grid is not None and tidx < data["precip"].shape[2]: - data["quorum"] -= 1 data["precip"][:, :, tidx] = grid + quorum["cnt"] -= 1 tidx = 0 shp = data["solar"].shape @@ -424,10 +424,10 @@ def _cb_lp(args): tidx += 1 pool.close() pool.join() - if data["quorum"] > 0: + if quorum["cnt"] > 0: LOG.warning( "Failed 75%% quorum with MRMS a2m %.1f, loading legacy", - data["quorum"], + quorum["cnt"], ) load_precip_legacy(data, dt, tile_bounds) LOG.info("finished precip calculation") diff --git a/scripts/cligen/r_factor.py b/scripts/cligen/r_factor.py index 86a12f97..b2875f3c 100644 --- a/scripts/cligen/r_factor.py +++ b/scripts/cligen/r_factor.py @@ -1,17 +1,21 @@ """R factor work.""" +import os +from multiprocessing import Pool + import cartopy.crs as ccrs import geopandas as gpd import matplotlib.colors as mpcolors import numpy as np import pandas as pd from matplotlib.patches import Polygon +from pyiem.database import get_sqlalchemy_conn from pyiem.plot import MapPlot from pyiem.plot.use_agg import plt -from pyiem.util import get_sqlalchemy_conn +from sqlalchemy import text from tqdm import tqdm -from pydep.io.wepp import read_cli +from pydep.rfactor import compute_rfactor_from_cli def plot(): @@ -66,45 +70,59 @@ def plot(): mp.fig.savefig("test.png") +def job(args) -> pd.DataFrame: + """Process.""" + (cliid, row) = args + if not os.path.isfile(row["filepath"]): + return cliid, pd.DataFrame() + return cliid, compute_rfactor_from_cli(row["filepath"]) + + def dump_data(): """Go main Go.""" - with get_sqlalchemy_conn("idep") as conn: - df = pd.read_sql( + inserts = 0 + with get_sqlalchemy_conn("idep") as conn, Pool() as pool: + clidf = pd.read_sql( """ - SELECT huc_12, max(climate_file) as cli from flowpaths where - scenario = 0 GROUP by huc_12 + SELECT id, filepath from climate_files where + scenario = 0 """, conn, - index_col="huc_12", + index_col="id", ) - data = { - "huc12": [], - "rfactor_yr_avg": [], - "rfactor_yr_max": [], - "rfactor_yr_min": [], - "rfactor_min_year": [], - "rfactor_max_year": [], - } - for year in range(2007, 2022): - data[f"rfactor_{year}"] = [] - for huc_12, row in tqdm(df.iterrows(), total=len(df.index)): - fn = row["cli"].replace("/i/", "/mnt/idep2/2/") - cdf = read_cli(fn, compute_rfactor=True) - data["huc12"].append(huc_12) - # drop 2020 - yearly = cdf["rfactor"].groupby(cdf.index.year).sum().iloc[:-1] - data["rfactor_yr_avg"].append(yearly.mean()) - data["rfactor_yr_min"].append(yearly.min()) - data["rfactor_yr_max"].append(yearly.max()) - data["rfactor_min_year"].append(yearly.idxmin()) - data["rfactor_max_year"].append(yearly.idxmax()) - for year in range(2007, 2022): - data[f"rfactor_{year}"].append(yearly.loc[year]) - - df = pd.DataFrame(data) - df.to_csv("/tmp/data.csv", index=False) + for clid, resdf in tqdm( + pool.imap(job, clidf.iterrows()), total=len(clidf) + ): + if resdf.empty: + continue + conn.execute( + text(""" + delete from climate_file_yearly_summary where climate_file_id = :clid + """), + {"clid": clid}, + ) + for year in range(2007, 2025): + conn.execute( + text( + """ + insert into climate_file_yearly_summary + (climate_file_id, year, rfactor, rfactor_storms) values + (:clid, :year, :rfactor, :rfactor_storms) + """ + ), + { + "clid": clid, + "year": year, + "rfactor": resdf.at[year, "rfactor"], + "rfactor_storms": resdf.at[year, "storm_count"], + }, + ) + inserts += 1 + if inserts > 1_000: + conn.commit() + inserts = 0 + conn.commit() if __name__ == "__main__": - # dump_data() - plot() + dump_data() diff --git a/scripts/dynamic_tillage/map_huc12_progress.py b/scripts/dynamic_tillage/map_huc12_progress.py new file mode 100644 index 00000000..13cfe279 --- /dev/null +++ b/scripts/dynamic_tillage/map_huc12_progress.py @@ -0,0 +1,133 @@ +"""Create a lapse of progress for a given year and HUC12.""" + +from datetime import date + +import click +import geopandas as gpd +import pandas as pd +from matplotlib.patches import Rectangle +from pyiem.database import get_sqlalchemy_conn +from pyiem.plot import MapPlot +from pyiem.util import logger +from sqlalchemy import text +from tqdm import tqdm + +LOG = logger() + + +def plot_map_progress( + huc12: str, + huc12df: gpd.GeoDataFrame, + fieldsdf: gpd.GeoDataFrame, + dt: date, + i: int, +): + """Make a map diagnostic""" + # Get the fields and their status + + minx, miny, maxx, maxy = huc12df["geom"].total_bounds + buffer = 0.01 + + mp = MapPlot( + title=f"DEP Planting Progress {dt:%Y %b %d}", + subtitle=f"For HUC12: {huc12df.iloc[0]['name']} [{huc12}]", + logo="dep", + sector="custom", + west=minx - buffer, + north=maxy + buffer, + south=miny - buffer, + east=maxx + buffer, + caption="Daily Erosion Project", + continentalcolor="white", + ) + huc12df.to_crs(mp.panels[0].crs).plot( + ax=mp.panels[0].ax, + aspect=None, + facecolor="None", + edgecolor="b", + linewidth=2, + zorder=100, + ) + fdf = fieldsdf + fdf["color"] = "yellow" + te = "tillage_events" + fdf.loc[fdf["plant"].notna(), "color"] = "white" + fdf.loc[fdf["till1"] <= dt, "color"] = "r" + fdf.loc[(fdf["till1"] <= dt) & (fdf[te] == 1), "color"] = "tan" + fdf.loc[(fdf["till2"] <= dt) & (fdf[te] == 2), "color"] = "tan" + fdf.loc[(fdf["till3"] <= dt) & (fdf[te] == 3), "color"] = "tan" + fieldsdf.loc[fieldsdf["plant"] <= dt, "color"] = "g" + fieldsdf.to_crs(mp.panels[0].crs).plot( + ax=mp.panels[0].ax, + aspect=None, + facecolor=fieldsdf["color"], + edgecolor="k", + linewidth=1, + zorder=100, + ) + pn = Rectangle((0, 0), 1, 1, fc="yellow", ec="k") + p0 = Rectangle((0, 0), 1, 1, fc="white", ec="k") + p1 = Rectangle((0, 0), 1, 1, fc="r", ec="k") + p2 = Rectangle((0, 0), 1, 1, fc="tan", ec="k") + p3 = Rectangle((0, 0), 1, 1, fc="g", ec="k") + mp.panels[0].ax.legend( + (pn, p0, p1, p2, p3), + ( + "Non-Ag", + "Awaiting", + "Tillage Started", + "Tillage Complete", + "Planted", + ), + ncol=2, + fontsize=11, + loc=2, + ) + + mp.fig.savefig(f"{i:04.0f}.png") + mp.close() + + +@click.command() +@click.option("--huc12", required=True, type=str, help="HUC12 to plot") +@click.option("--year", required=True, type=int, help="Year to plot") +def main(huc12: str, year: int): + """Go Main Go.""" + with get_sqlalchemy_conn("idep") as conn: + huc12df = gpd.read_postgis( + text(""" + select huc_12, name, ST_Transform(geom, 4326) as geom + from huc12 where huc_12 = :huc12 and scenario = 0 + """), + conn, + params={"huc12": huc12}, + index_col="huc_12", + geom_col="geom", + crs="EPSG:4326", + ) + fieldsdf = gpd.read_postgis( + text(""" + select o.*, ST_Transform(f.geom, 4326) as geom, 0 as tillage_events + from fields f LEFT JOIN field_operations o + on (f.field_id = o.field_id and o.year = :year) + where huc12 = :huc12 + """), + conn, + params={"year": year, "huc12": huc12}, + index_col="field_id", + parse_dates=["till1", "till2", "till3", "plant"], + geom_col="geom", + crs="EPSG:4326", + ) + fieldsdf.loc[fieldsdf["till1"].notna(), "tillage_events"] = 1 + fieldsdf.loc[fieldsdf["till2"].notna(), "tillage_events"] = 2 + fieldsdf.loc[fieldsdf["till3"].notna(), "tillage_events"] = 3 + LOG.info("Found %s fields for plotting", len(fieldsdf)) + progress = tqdm(pd.date_range(f"{year}-04-11", f"{year}-06-14")) + for i, dt in enumerate(progress): + progress.set_description(f"Plotting {dt}") + plot_map_progress(huc12, huc12df, fieldsdf, dt, i) + + +if __name__ == "__main__": + main() diff --git a/src/pydep/rfactor.py b/src/pydep/rfactor.py index fc71829b..33e4eb9c 100644 --- a/src/pydep/rfactor.py +++ b/src/pydep/rfactor.py @@ -1,5 +1,7 @@ """pydep's R-factor calculations.""" +from datetime import date + import numpy as np import pandas as pd from scipy.interpolate import interp1d @@ -63,8 +65,8 @@ def _compute_accum(lines): )(np.arange(0.5, 24.01, 0.5)) -def _rfactor_window(accum: np.ndarray) -> list: - """Compute R-factors for the last day in the 10 day window.""" +def _rfactor_year(accum: np.ndarray) -> list: + """Compute R-factors for a year.""" storm_start_stop = [] storm_start = -1 x = 12 @@ -80,9 +82,6 @@ def _rfactor_window(accum: np.ndarray) -> list: x += 1 rfactors = [] for start, stop in storm_start_stop: - # For a storm to be considered, it must stop within the tenth day - if stop < (9 * 24 * 2): - continue event = accum[start:stop] - accum[start] rate_mmhr = (event[1:] - event[0:-1]) * 2.0 # sum of E x I @@ -119,26 +118,29 @@ def compute_rfactor_from_cli( range(int(year1), int(year1) + int(years_simuluated)), name="year" ), ) - # Life Choice, store 10 days of 30 minute accum - window_accum = np.zeros(10 * 24 * 2) + # We compute one year at a time. + year_accum = np.array([]) while lnum < len(lines): - (_da, _mo, year, breakpoints, *_bah) = lines[lnum].strip().split() + (da, mo, year, breakpoints, *_bah) = lines[lnum].strip().split() + dt = date(int(year), int(mo), int(da)) + doy = dt.timetuple().tm_yday + idx = (doy - 1) * 24 * 2 + if doy == 1: + year_accum = np.zeros(366 * 24 * 2) breakpoints = int(breakpoints) if breakpoints == 0: accum = np.zeros(24 * 2) else: accum = _compute_accum(lines[lnum + 1 : lnum + 1 + breakpoints]) - # Chop off the oldest day and add the new day - window_accum = np.concatenate( - ( - window_accum[48:] - window_accum[48], - accum + window_accum[-1] - window_accum[48], - ) + # Place this accum into the yearly accum and apply offset to keep mono + year_accum[idx : idx + 24 * 2] = accum + ( + 0 if doy == 1 else year_accum[idx - 1] ) - rfactors = _rfactor_window(window_accum) - if rfactors: - resultdf.at[int(year), "rfactor"] += np.sum(rfactors) - resultdf.at[int(year), "storm_count"] += len(rfactors) + if f"{dt:%m%d}" == "1231": + rfactors = _rfactor_year(year_accum) + if rfactors: + resultdf.at[dt.year, "rfactor"] = np.sum(rfactors) + resultdf.at[dt.year, "storm_count"] = len(rfactors) lnum += breakpoints + 1 return resultdf diff --git a/tests/io/test_rfactor.py b/tests/test_rfactor.py similarity index 67% rename from tests/io/test_rfactor.py rename to tests/test_rfactor.py index e75ce2c0..a6e222f6 100644 --- a/tests/io/test_rfactor.py +++ b/tests/test_rfactor.py @@ -8,10 +8,12 @@ def get_path(name): """helper""" basedir = os.path.dirname(__file__) - return os.path.join(basedir, "..", "data", name) + return os.path.join(basedir, "data", name) def test_rfactor(): """Walk before we run.""" resultdf = compute_rfactor_from_cli(get_path("cli.txt")) + # There is no right answer, but we can check for changes + assert abs(resultdf.at[2007, "rfactor"] - 1252.09) < 0.1 assert resultdf.at[2007, "storm_count"] == 76