From db1e7c4363ac1617b3e78962537eed69f33d4429 Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 11 Apr 2024 12:06:32 -0500 Subject: [PATCH 1/7] set dyntillage mud it in date to 10 June per #167 --- scripts/dynamic_tillage/workflow.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py index 4d976a93..59d724f0 100644 --- a/scripts/dynamic_tillage/workflow.py +++ b/scripts/dynamic_tillage/workflow.py @@ -176,7 +176,8 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: if pd.isnull(acres_tilled): acres_tilled = 0 - if dt.month < 6: + mud_it_in = f"{dt:%m%d}" >= "0610" + if not mud_it_in: # Soil Moisture check. We are not modeling every field, so we are a # bit generous here, are 50% of **model** fields above plastic_limit? fields["sw1"] = pd.Series(fieldsw) @@ -198,7 +199,7 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: fields["tillage"] = fields["management"].str.slice(char_at - 1, char_at) total_acres = fields["acres"].sum() - limit = total_acres / 10.0 if dt.month < 6 else total_acres + 1 + limit = (total_acres / 10.0) if not mud_it_in else total_acres + 1 # Work on tillage first, so to avoid planting on tilled fields for fbndid, row in fields[fields["till_needed"]].iterrows(): @@ -244,11 +245,12 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: ) cursor.close() conn.commit() - return acres_planted, acres_tilled # Update the .rot files, when necessary - df2 = df[(df["field_id"].isin(df2["field_id"])) & (df["ofe"].notna())] - for _, row in df2.iterrows(): - edit_rotfile(dt.year, huc12, row) + if dt.month == 6 and dt.day == 14: + df2 = df[(df["field_id"].isin(df2["field_id"])) & (df["ofe"].notna())] + for _, row in df2.iterrows(): + edit_rotfile(dt.year, huc12, row) + return acres_planted, acres_tilled for fpath in df2["fpath"].unique(): prj2wepp(huc12, fpath) @@ -309,8 +311,8 @@ def main(scenario, dt, huc12): geom_col="geom", index_col="huc_12", ) - # In June, we mud it in. - if dt.month < 6: + # After June 10th, we mud it in. + if f"{dt:%m%d}" < "0610": # Estimate today's rainfall so to delay triggers estimate_rainfall(huc12df, dt) # Any HUC12 over 10mm gets skipped From 57f083f08564d428340ad413723e8c1424676d1e Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 11 Apr 2024 12:08:08 -0500 Subject: [PATCH 2/7] set dyntill plastic_limit check to 0.8*PL --- scripts/dynamic_tillage/workflow.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py index 59d724f0..83ca0dd9 100644 --- a/scripts/dynamic_tillage/workflow.py +++ b/scripts/dynamic_tillage/workflow.py @@ -179,11 +179,12 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: mud_it_in = f"{dt:%m%d}" >= "0610" if not mud_it_in: # Soil Moisture check. We are not modeling every field, so we are a - # bit generous here, are 50% of **model** fields above plastic_limit? + # bit generous here, are 50% of **model** fields above + # (plastic_limit * 0.8) ? fields["sw1"] = pd.Series(fieldsw) modelled = (fields["sw1"] > 0).sum() if modelled > 1: - hits = (fields["sw1"] > fields["plastic_limit"]).sum() + hits = (fields["sw1"] > (fields["plastic_limit"] * 0.8)).sum() if hits > modelled / 2: LOG.debug( "HUC12: %s has %s/%s fields above plastic limit", From e6a67584467c9b1d8d771e679b7300f7c9dc4440 Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 11 Apr 2024 22:49:19 -0500 Subject: [PATCH 3/7] fix: [dyntill] gfs soilt robustness --- scripts/dynamic_tillage/workflow.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py index 83ca0dd9..6b91d441 100644 --- a/scripts/dynamic_tillage/workflow.py +++ b/scripts/dynamic_tillage/workflow.py @@ -13,7 +13,7 @@ import shutil import subprocess import threading -from datetime import timedelta +from datetime import datetime, timedelta from multiprocessing import cpu_count from multiprocessing.pool import Pool from typing import Tuple @@ -258,13 +258,22 @@ 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}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"][:]) - x = np.digitize(huc12df["lon"].values, nc.variables["lon"][:]) - for i, idx in enumerate(huc12df.index.values): - huc12df.at[idx, "tsoil"] = tsoil[y[i], x[i]] + # Look back from 12z at most 24 hours for a file to use + z12 = datetime(dt.year, dt.month, dt.day, 12) + for offset in range(0, 25, 6): + now = z12 - timedelta(hours=offset) + ppath = f"{now:%Y/%m/%d}/model/gfs/gfs_{now:%Y%m%d%H}_iemre.nc" + with archive_fetch(ppath) as ncfn: + if ncfn is None: + continue + LOG.info("Using GFS@%s for soil temperature", now) + with 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"][:]) + x = np.digitize(huc12df["lon"].values, nc.variables["lon"][:]) + for i, idx in enumerate(huc12df.index.values): + huc12df.at[idx, "tsoil"] = tsoil[y[i], x[i]] + break def estimate_rainfall(huc12df, dt): @@ -303,6 +312,7 @@ def main(scenario, dt, huc12): false as limited_by_soiltemp, false as limited_by_soilmoisture, false as limited, + 99 as tsoil, ST_x(st_transform(st_centroid(geom), 4326)) as lon, ST_y(st_transform(st_centroid(geom), 4326)) as lat from huc12 where scenario = :scenario {huc12_filter} From 0e96c3881c017a84470e517f3df3f470a78d07fb Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 11 Apr 2024 22:50:21 -0500 Subject: [PATCH 4/7] feat: [dyntill] set max rate to 6% --- scripts/dynamic_tillage/workflow.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py index 6b91d441..e5114348 100644 --- a/scripts/dynamic_tillage/workflow.py +++ b/scripts/dynamic_tillage/workflow.py @@ -200,7 +200,8 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: fields["tillage"] = fields["management"].str.slice(char_at - 1, char_at) total_acres = fields["acres"].sum() - limit = (total_acres / 10.0) if not mud_it_in else total_acres + 1 + # NB: Crude assessment of NASS peak daily planting rate, was 10% + limit = (total_acres * 0.06) if not mud_it_in else total_acres + 1 # Work on tillage first, so to avoid planting on tilled fields for fbndid, row in fields[fields["till_needed"]].iterrows(): From 091933023339d913a94a8a11667a7d8eca36c762 Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 12 Apr 2024 16:03:01 -0500 Subject: [PATCH 5/7] [WEPP] increase mxtlsq=300 for more tillage dates --- src/README.md | 1 + src/wepp20230117/pmxtls.inc | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/README.md b/src/README.md index 95bd8d44..872bd91a 100644 --- a/src/README.md +++ b/src/README.md @@ -13,3 +13,4 @@ The following code changes were made to the baseline WEPP model. with high resolution time precision?!? - Change WEPP `.env` files to have explicit years (#183). - [sedout.for] Make `.env` output single space delimited (#208). +- [pmxtls.inc] Increase `mxtlsq=300` to support dynamic tillage dates. diff --git a/src/wepp20230117/pmxtls.inc b/src/wepp20230117/pmxtls.inc index ac2f3595..ba5b084e 100644 --- a/src/wepp20230117/pmxtls.inc +++ b/src/wepp20230117/pmxtls.inc @@ -4,7 +4,7 @@ c begin include file pmxtls.inc c + + + PARAMETER DECLARATIONS + + + integer mxtlsq - parameter (mxtlsq=30) + parameter (mxtlsq=300) c + + + PARAMETER DEFINITIONS + + + From eafba53e58130c5900a9bc3534d00073dd85e27b Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 12 Apr 2024 16:03:39 -0500 Subject: [PATCH 6/7] feat: better support june dates for dyn tillage #167 --- scripts/dynamic_tillage/workflow.py | 39 ++++++++++++++++------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py index e5114348..322ab41e 100644 --- a/scripts/dynamic_tillage/workflow.py +++ b/scripts/dynamic_tillage/workflow.py @@ -156,7 +156,11 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: return 0, 0 # Screen out when there are no fields that need planting or tilling - if not df["plant_needed"].any() and not df["till_needed"].any(): + if ( + not df["plant_needed"].any() + and not df["till_needed"].any() + and f"{dt:%m%d}" < "0614" + ): LOG.debug("No fields need planting or tilling for %s", huc12) return 0, 0 @@ -235,23 +239,24 @@ def do_huc12(dt, huc12, fieldsw) -> Tuple[int, int]: # Update the database df2 = fields[fields["operation_done"]].reset_index() - if df2.empty: - return acres_planted, acres_tilled - conn, cursor = get_dbconnc("idep") - cursor.executemany( - "update field_operations set plant = %s, till1 = %s, " - "till2 = %s, till3 = %s where field_id = %s and year = %s", - df2[ - ["plant", "till1", "till2", "till3", "field_id", "year"] - ].itertuples(index=False), - ) - cursor.close() - conn.commit() - # Update the .rot files, when necessary + if not df2.empty: + conn, cursor = get_dbconnc("idep") + cursor.executemany( + "update field_operations set plant = %s, till1 = %s, " + "till2 = %s, till3 = %s where field_id = %s and year = %s", + df2[ + ["plant", "till1", "till2", "till3", "field_id", "year"] + ].itertuples(index=False), + ) + cursor.close() + conn.commit() + # Update all the .rot files if dt.month == 6 and dt.day == 14: + df2 = df[df["ofe"].notna()] + elif not df2.empty: df2 = df[(df["field_id"].isin(df2["field_id"])) & (df["ofe"].notna())] - for _, row in df2.iterrows(): - edit_rotfile(dt.year, huc12, row) + for _, row in df2.iterrows(): + edit_rotfile(dt.year, huc12, row) return acres_planted, acres_tilled for fpath in df2["fpath"].unique(): prj2wepp(huc12, fpath) @@ -313,7 +318,7 @@ def main(scenario, dt, huc12): false as limited_by_soiltemp, false as limited_by_soilmoisture, false as limited, - 99 as tsoil, + 99.9 as tsoil, ST_x(st_transform(st_centroid(geom), 4326)) as lon, ST_y(st_transform(st_centroid(geom), 4326)) as lat from huc12 where scenario = :scenario {huc12_filter} From ebf5775b59d9b69d478804ef187c8e309e96a749 Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 12 Apr 2024 16:04:06 -0500 Subject: [PATCH 7/7] chore: bump ruff --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7fc4355c..793b8658 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,6 +1,6 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.3.5" + rev: "v0.3.7" hooks: - id: ruff args: [--fix, --exit-non-zero-on-fix]