Merge pull request dailyerosion#308 from akrherz/241216

Omnibus

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

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

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

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

tests/io/test_rfactor.py → 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