From 84afdb5d7be56ef79be4a0db789128e38156a446 Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 19 Dec 2024 11:27:13 -0600 Subject: [PATCH 1/4] =?UTF-8?q?=E2=AC=86=EF=B8=8F=20Update=20ruff?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .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 09246f9eb..6ab457d70 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,7 +6,7 @@ repos: hooks: - id: eslint - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.8.3" + rev: "v0.8.4" hooks: - id: ruff args: [--fix, --exit-non-zero-on-fix] From 78a1501ac375c0441377d7b8e68c89490be5c621 Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 19 Dec 2024 14:26:27 -0600 Subject: [PATCH 2/4] =?UTF-8?q?=F0=9F=8E=A8=20More=20grid=20updates=20per?= =?UTF-8?q?=20pyIEM=20updates?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- pylib/iemweb/autoplot/scripts/p4.py | 2 +- pylib/iemweb/autoplot/scripts/p89.py | 2 +- pylib/iemweb/autoplot/scripts100/p175.py | 2 +- pylib/iemweb/autoplot/scripts100/p182.py | 2 +- pylib/iemweb/autoplot/scripts100/p185.py | 2 +- pylib/iemweb/iemre/daily.py | 5 +- pylib/iemweb/iemre/multiday.py | 6 +- pylib/iemweb/json/stage4.py | 42 +++++------ scripts/climodat/compute4regions.py | 2 +- scripts/climodat/era5land_extract.py | 26 +++---- scripts/climodat/merra_solarrad.py | 3 +- scripts/current/ifc_today_total.py | 2 +- scripts/current/mrms_today_total.py | 2 +- scripts/current/q3_xhour.py | 2 +- scripts/era5/fetch_era5.py | 94 ++++++++++++++++-------- scripts/era5/init_hourly.py | 40 +++++----- scripts/gfs/gfs2iemre.py | 17 ++--- scripts/hrrr/hrrr_ref2raster.py | 1 + scripts/iemre/grid_rsds.py | 9 +-- scripts/iemre/hourly_analysis.py | 17 +---- scripts/iemre/init_stage4_daily.py | 28 +++++++ scripts/iemre/init_stage4_dailyc.py | 28 +++++++ scripts/iemre/init_stage4_hourly.py | 28 ++++++- scripts/iemre/precip_ingest.py | 21 ++---- scripts/iemre/prism_adjust_stage4.py | 80 ++++++++++---------- scripts/iemre/stage4_12z_adjust.py | 13 ++-- scripts/mrms/init_daily_mrms.py | 40 ++++------ scripts/mrms/init_mrms_dailyc.py | 22 ++---- scripts/mrms/merge_mrms_q3.py | 10 +-- scripts/swat/swat_realtime.py | 6 +- 30 files changed, 303 insertions(+), 251 deletions(-) diff --git a/pylib/iemweb/autoplot/scripts/p4.py b/pylib/iemweb/autoplot/scripts/p4.py index ecf87f16f..f0aba9727 100644 --- a/pylib/iemweb/autoplot/scripts/p4.py +++ b/pylib/iemweb/autoplot/scripts/p4.py @@ -74,7 +74,7 @@ def plotter(fdict): raise NoDataFound("Data not available for year") with ncopen(ncfn) as nc: precip = nc.variables["p01d"] - czs = CachingZonalStats(iemre.AFFINE) + czs = CachingZonalStats(iemre.DOMAINS[""]["affine"]) hasdata = np.zeros( (nc.dimensions["lat"].size, nc.dimensions["lon"].size) ) diff --git a/pylib/iemweb/autoplot/scripts/p89.py b/pylib/iemweb/autoplot/scripts/p89.py index 7bdd81675..0d3047b15 100644 --- a/pylib/iemweb/autoplot/scripts/p89.py +++ b/pylib/iemweb/autoplot/scripts/p89.py @@ -99,7 +99,7 @@ def get_data(ctx): raise NoDataFound(f"Missing {ncfn}") with ncopen(ncfn) as nc: precip = nc.variables["p01d"] - czs = CachingZonalStats(iemre.AFFINE) + czs = CachingZonalStats(iemre.DOMAINS[""]["affine"]) hasdata = np.zeros( (nc.dimensions["lat"].size, nc.dimensions["lon"].size) ) diff --git a/pylib/iemweb/autoplot/scripts100/p175.py b/pylib/iemweb/autoplot/scripts100/p175.py index 716344b47..338f09228 100644 --- a/pylib/iemweb/autoplot/scripts100/p175.py +++ b/pylib/iemweb/autoplot/scripts100/p175.py @@ -81,7 +81,7 @@ def plotter(fdict): if not os.path.isfile(ncfn): raise NoDataFound(f"Data for year {sts.year} not found") with ncopen(ncfn) as nc: - czs = CachingZonalStats(iemre.AFFINE) + czs = CachingZonalStats(iemre.DOMAINS[""]["affine"]) hasdata = np.zeros( (nc.dimensions["lat"].size, nc.dimensions["lon"].size) ) diff --git a/pylib/iemweb/autoplot/scripts100/p182.py b/pylib/iemweb/autoplot/scripts100/p182.py index 02e86da32..e8b765135 100644 --- a/pylib/iemweb/autoplot/scripts100/p182.py +++ b/pylib/iemweb/autoplot/scripts100/p182.py @@ -17,7 +17,7 @@ from pyiem import iemre, util from pyiem.exceptions import NoDataFound from pyiem.grid.zs import CachingZonalStats -from pyiem.iemre import MRMS4IEMRE_AFFINE +from pyiem.mrms import MRMS4IEMRE_AFFINE from pyiem.plot import figure_axes from pyiem.reference import state_names diff --git a/pylib/iemweb/autoplot/scripts100/p185.py b/pylib/iemweb/autoplot/scripts100/p185.py index 72224b5ac..720654085 100644 --- a/pylib/iemweb/autoplot/scripts100/p185.py +++ b/pylib/iemweb/autoplot/scripts100/p185.py @@ -14,7 +14,7 @@ from pyiem import iemre, util from pyiem.exceptions import NoDataFound from pyiem.grid.zs import CachingZonalStats -from pyiem.iemre import MRMS4IEMRE_AFFINE +from pyiem.mrms import MRMS4IEMRE_AFFINE from pyiem.plot import get_cmap from pyiem.plot.geoplot import MapPlot diff --git a/pylib/iemweb/iemre/daily.py b/pylib/iemweb/iemre/daily.py index 4cc17eb64..82c98d492 100644 --- a/pylib/iemweb/iemre/daily.py +++ b/pylib/iemweb/iemre/daily.py @@ -26,7 +26,7 @@ import numpy as np import pyiem.prism as prismutil -from pyiem import iemre +from pyiem import iemre, mrms from pyiem.util import convert_value, ncopen from pyiem.webutil import iemapp @@ -85,8 +85,7 @@ def application(environ, start_response): if not os.path.isfile(ncfn): mrms_precip = None else: - j2 = int((lat - iemre.SOUTH) * 100.0) - i2 = int((lon - iemre.WEST) * 100.0) + i2, j2 = mrms.find_ij(lon, lat) with ncopen(ncfn) as nc: mrms_precip = nc.variables["p01d"][offset, j2, i2] / 25.4 else: diff --git a/pylib/iemweb/iemre/multiday.py b/pylib/iemweb/iemre/multiday.py index aac40768e..b8ef7f5a2 100644 --- a/pylib/iemweb/iemre/multiday.py +++ b/pylib/iemweb/iemre/multiday.py @@ -25,7 +25,7 @@ import numpy as np import pyiem.prism as prismutil from pydantic import Field -from pyiem import iemre +from pyiem import iemre, mrms from pyiem.util import convert_value, ncopen from pyiem.webutil import CGIModel, iemapp @@ -98,7 +98,6 @@ def application(environ, start_response): "Point not within any domain", ) ] - dom = iemre.DOMAINS[domain] i, j = iemre.find_ij(lon, lat, domain=domain) if i is None or j is None: return [ @@ -166,8 +165,7 @@ def application(environ, start_response): prism_precip = [None] * (offset2 - offset1) if ts1.year > 2000 and domain == "": - j2 = int((lat - dom["south"]) * 100.0) - i2 = int((lon - dom["west"]) * 100.0) + i2, j2 = mrms.find_ij(lon, lat) with ncopen(iemre.get_daily_mrms_ncname(ts1.year)) as nc: mrms_precip = nc.variables["p01d"][tslice, j2, i2] / 25.4 else: diff --git a/pylib/iemweb/json/stage4.py b/pylib/iemweb/json/stage4.py index f4b7b6aad..06d059e25 100644 --- a/pylib/iemweb/json/stage4.py +++ b/pylib/iemweb/json/stage4.py @@ -35,6 +35,7 @@ from pydantic import Field, field_validator from pyiem import iemre from pyiem.reference import ISO8601 +from pyiem.stage4 import find_ij from pyiem.util import mm2inch, ncopen, utc from pyiem.webutil import CGIModel, iemapp @@ -43,8 +44,8 @@ class Schema(CGIModel): """See how we are called.""" callback: str = Field(None, description="JSONP callback function name") - lat: float = Field(..., description="Latitude of point") - lon: float = Field(..., description="Longitude of point") + lat: float = Field(..., description="Latitude of point", ge=-90, le=90) + lon: float = Field(..., description="Longitude of point", ge=-180, le=180) valid: date = Field(..., description="Valid date of data") tz: str = Field("UTC", description="Timezone of valid date") @@ -88,26 +89,23 @@ def dowork(environ): } if not os.path.isfile(ncfn): return json.dumps(res) - with ncopen(ncfn) as nc: - dist = ( - (nc.variables["lon"][:] - environ["lon"]) ** 2 - + (nc.variables["lat"][:] - environ["lat"]) ** 2 - ) ** 0.5 - (j, i) = np.unravel_index(dist.argmin(), dist.shape) # skipcq - res["gridi"] = int(i) - res["gridj"] = int(j) - - ppt = nc.variables["p01m"][sidx:eidx, j, i] - - for tx, pt in enumerate(ppt): - valid = sts + timedelta(hours=tx) - utcnow = valid.astimezone(ZoneInfo("UTC")) - res["data"].append( - { - "end_valid": utcnow.strftime("%Y-%m-%dT%H:00:00Z"), - "precip_in": myrounder(mm2inch(pt), 2), - } - ) + i, j = find_ij(environ["lon"], environ["lat"]) + if i is not None: + with ncopen(ncfn) as nc: + res["gridi"] = i + res["gridj"] = j + + ppt = nc.variables["p01m"][sidx:eidx, j, i] + + for tx, pt in enumerate(ppt): + valid = sts + timedelta(hours=tx) + utcnow = valid.astimezone(ZoneInfo("UTC")) + res["data"].append( + { + "end_valid": utcnow.strftime("%Y-%m-%dT%H:00:00Z"), + "precip_in": myrounder(mm2inch(pt), 2), + } + ) return json.dumps(res) diff --git a/scripts/climodat/compute4regions.py b/scripts/climodat/compute4regions.py index 6ef4094f5..c68e4e7a8 100644 --- a/scripts/climodat/compute4regions.py +++ b/scripts/climodat/compute4regions.py @@ -83,7 +83,7 @@ def do_day(cursor, valid): index_col="id", geom_col="geom", ) - czs = CachingZonalStats(iemre.AFFINE) + czs = CachingZonalStats(iemre.DOMAINS[""]["affine"]) sthigh = czs.gen_stats(np.flipud(high), gdf["geom"]) stlow = czs.gen_stats(np.flipud(low), gdf["geom"]) stprecip = czs.gen_stats(np.flipud(precip), gdf["geom"]) diff --git a/scripts/climodat/era5land_extract.py b/scripts/climodat/era5land_extract.py index 5cc6cbf13..facdd4532 100644 --- a/scripts/climodat/era5land_extract.py +++ b/scripts/climodat/era5land_extract.py @@ -10,10 +10,10 @@ import geopandas as gpd import numpy as np import pandas as pd -from affine import Affine from pyiem.database import get_dbconn, get_sqlalchemy_conn +from pyiem.era5land import DOMAINS, find_ij from pyiem.grid.zs import CachingZonalStats -from pyiem.iemre import NORTH, WEST, hourly_offset +from pyiem.iemre import hourly_offset from pyiem.util import convert_value, logger, ncopen, utc LOG = logger() @@ -31,9 +31,8 @@ def compute_regions(data, varname, df): index_col="id", geom_col="geom", ) - affine = Affine(0.1, 0, WEST, 0, -0.1, NORTH) - czs = CachingZonalStats(affine) - data = czs.gen_stats(np.flipud(data), gdf["geom"]) + czs = CachingZonalStats(DOMAINS[""]["AFFINE_NC"]) + data = czs.gen_stats(data, gdf["geom"]) for i, sid in enumerate(gdf.index.values): df.at[sid, varname] = data[i] @@ -61,8 +60,6 @@ def build_stations(dt) -> pd.DataFrame: "era5land_soilm1m_avg", ]: df[col] = np.nan - df["i"] = np.nan - df["j"] = np.nan LOG.info("Found %s database entries", len(df.index)) return df @@ -79,8 +76,6 @@ def compute(df, sids, dt, do_regions=False): # Wm-2 to MJ factor = 3600.0 / 1_000_000.0 with ncopen(ncfn) as nc: - lons = nc.variables["lon"][:] - lats = nc.variables["lat"][:] if f"{dt:%m%d}" == "1231": rsds = np.sum(nc.variables["rsds"][idx0:], 0) * factor # Close enough @@ -104,20 +99,17 @@ def compute(df, sids, dt, do_regions=False): ) / 100.0 soilt = np.mean(nc.variables["soilt"][idx0:idx1, 0], 0) - df["i"] = np.digitize(df["lon"].values, lons) - df["j"] = np.digitize(df["lat"].values, lats) rsds = rsds.filled(np.nan) soilm = soilm.filled(np.nan) soilm1m = soilm1m.filled(np.nan) soilt = soilt.filled(np.nan) for sid, row in df.loc[sids].iterrows(): - df.at[sid, "era5land_srad"] = rsds[int(row["j"]), int(row["i"])] - df.at[sid, "era5land_soilt4_avg"] = soilt[int(row["j"]), int(row["i"])] - df.at[sid, "era5land_soilm4_avg"] = soilm[int(row["j"]), int(row["i"])] - df.at[sid, "era5land_soilm1m_avg"] = soilm1m[ - int(row["j"]), int(row["i"]) - ] + i, j = find_ij(row["lon"], row["lat"]) + df.at[sid, "era5land_srad"] = rsds[j, i] + df.at[sid, "era5land_soilt4_avg"] = soilt[j, i] + df.at[sid, "era5land_soilm4_avg"] = soilm[j, i] + df.at[sid, "era5land_soilm1m_avg"] = soilm1m[j, i] if do_regions: compute_regions(rsds, "era5land_srad", df) diff --git a/scripts/climodat/merra_solarrad.py b/scripts/climodat/merra_solarrad.py index be718f0b9..fbf83a426 100644 --- a/scripts/climodat/merra_solarrad.py +++ b/scripts/climodat/merra_solarrad.py @@ -38,7 +38,8 @@ def compute_regions(rsds, df): index_col="id", geom_col="geom", ) - affine = Affine(0.625, 0, -180.0, 0, -0.5, 90) + # This may not be exactly right, but alas + affine = Affine(0.625, 0, -180.3125, 0, -0.5, 90.5) czs = CachingZonalStats(affine) data = czs.gen_stats(np.flipud(rsds), gdf["geom"]) for i, sid in enumerate(gdf.index.values): diff --git a/scripts/current/ifc_today_total.py b/scripts/current/ifc_today_total.py index fc8009fa6..74a8ed046 100644 --- a/scripts/current/ifc_today_total.py +++ b/scripts/current/ifc_today_total.py @@ -13,7 +13,7 @@ def doday(ts, realtime): """ - Create a plot of precipitation stage4 estimates for some day + Create a plot of precipitation IFC estimates for some day """ idx = daily_offset(ts) with ncopen( diff --git a/scripts/current/mrms_today_total.py b/scripts/current/mrms_today_total.py index bc5bc7426..147a6c1eb 100644 --- a/scripts/current/mrms_today_total.py +++ b/scripts/current/mrms_today_total.py @@ -19,7 +19,7 @@ def doday(ts, realtime): """ - Create a plot of precipitation stage4 estimates for some day + Create a plot of precipitation MRMS estimates for some day """ lts = utc(ts.year, ts.month, ts.day, 12) lts = lts.astimezone(ZoneInfo("America/Chicago")) diff --git a/scripts/current/q3_xhour.py b/scripts/current/q3_xhour.py index 8c7c46d08..71f4a2e3c 100644 --- a/scripts/current/q3_xhour.py +++ b/scripts/current/q3_xhour.py @@ -22,7 +22,7 @@ def doit(ts, hours): """ - Create a plot of precipitation stage4 estimates for some day + Create a plot of precipitation MRMS estimates for some day """ # Start at 1 AM ts = ts.replace(minute=0, second=0, microsecond=0) diff --git a/scripts/era5/fetch_era5.py b/scripts/era5/fetch_era5.py index 0dc190503..85b78b5a2 100644 --- a/scripts/era5/fetch_era5.py +++ b/scripts/era5/fetch_era5.py @@ -11,7 +11,8 @@ import cdsapi import click import numpy as np -from pyiem import iemre +from netCDF4 import Dataset +from pyiem import era5land, iemre from pyiem.util import logger, ncopen, utc LOG = logger() @@ -33,23 +34,45 @@ warnings.simplefilter("ignore", RuntimeWarning) -def ingest(ncin, nc, valid, domain): +def ingest(ncin: Dataset, nc: Dataset, valid, domain): """Consume this grib file.""" + # Bad daryl had an off by one error here + ijs = slice(None, None) + if ncin.variables["longitude"].size != nc.variables["lon"].size: + LOG.info("Applying off by one logic") + ijs = slice(0, -1) + tidx = iemre.hourly_offset(valid) dd = "" if domain == "" else f"_{domain}" for ekey, key in VERBATIM.items(): - nc.variables[key][tidx] = np.flipud(ncin.variables[ekey][0]) + nc.variables[key][tidx, ijs, ijs] = np.flipud(ncin.variables[ekey][0]) - nc.variables["soilt"][tidx, 0] = np.flipud(ncin.variables["stl1"][0]) - nc.variables["soilt"][tidx, 1] = np.flipud(ncin.variables["stl2"][0]) - nc.variables["soilt"][tidx, 2] = np.flipud(ncin.variables["stl3"][0]) - nc.variables["soilt"][tidx, 3] = np.flipud(ncin.variables["stl4"][0]) + nc.variables["soilt"][tidx, 0, ijs, ijs] = np.flipud( + ncin.variables["stl1"][0] + ) + nc.variables["soilt"][tidx, 1, ijs, ijs] = np.flipud( + ncin.variables["stl2"][0] + ) + nc.variables["soilt"][tidx, 2, ijs, ijs] = np.flipud( + ncin.variables["stl3"][0] + ) + nc.variables["soilt"][tidx, 3, ijs, ijs] = np.flipud( + ncin.variables["stl4"][0] + ) - nc.variables["soilm"][tidx, 0] = np.flipud(ncin.variables["swvl1"][0]) - nc.variables["soilm"][tidx, 1] = np.flipud(ncin.variables["swvl2"][0]) - nc.variables["soilm"][tidx, 2] = np.flipud(ncin.variables["swvl3"][0]) - nc.variables["soilm"][tidx, 3] = np.flipud(ncin.variables["swvl4"][0]) + nc.variables["soilm"][tidx, 0, ijs, ijs] = np.flipud( + ncin.variables["swvl1"][0] + ) + nc.variables["soilm"][tidx, 1, ijs, ijs] = np.flipud( + ncin.variables["swvl2"][0] + ) + nc.variables["soilm"][tidx, 2, ijs, ijs] = np.flipud( + ncin.variables["swvl3"][0] + ) + nc.variables["soilm"][tidx, 3, ijs, ijs] = np.flipud( + ncin.variables["swvl4"][0] + ) # -- these vars are accumulated since 0z, so 0z is the 24hr sum rsds = nc.variables["rsds"] @@ -63,39 +86,44 @@ def ingest(ncin, nc, valid, domain): f"/mesonet/data/era5{dd}/{valid.year - 1}_era5land_hourly.nc" ) as nc2: tsolar = ( - np.sum(nc2.variables["rsds"][(tidx0 + 1) :], 0) * 3600.0 + np.sum(nc2.variables["rsds"][(tidx0 + 1) :, ijs, ijs], 0) + * 3600.0 + ) + tp01m = np.sum( + nc2.variables["p01m"][(tidx0 + 1) :, ijs, ijs], 0 + ) + tevap = np.sum( + nc2.variables["evap"][(tidx0 + 1) :, ijs, ijs], 0 ) - tp01m = np.sum(nc2.variables["p01m"][(tidx0 + 1) :], 0) - tevap = np.sum(nc2.variables["evap"][(tidx0 + 1) :], 0) else: - tsolar = np.sum(rsds[(tidx0 + 1) : tidx], 0) * 3600.0 - tp01m = np.sum(p01m[(tidx0 + 1) : tidx], 0) - tevap = np.sum(evap[(tidx0 + 1) : tidx], 0) + tsolar = np.sum(rsds[(tidx0 + 1) : tidx, ijs], 0) * 3600.0 + tp01m = np.sum(p01m[(tidx0 + 1) : tidx, ijs], 0) + tevap = np.sum(evap[(tidx0 + 1) : tidx, ijs], 0) elif valid.hour > 1: tidx0 = iemre.hourly_offset(valid.replace(hour=1)) - tsolar = np.sum(rsds[tidx0:tidx], 0) * 3600.0 - tp01m = np.sum(p01m[tidx0:tidx], 0) - tevap = np.sum(evap[tidx0:tidx], 0) + tsolar = np.sum(rsds[tidx0:tidx, ijs], 0) * 3600.0 + tp01m = np.sum(p01m[tidx0:tidx, ijs], 0) + tevap = np.sum(evap[tidx0:tidx, ijs], 0) else: - tsolar = np.zeros(rsds.shape[1:]) - tp01m = np.zeros(rsds.shape[1:]) - tevap = np.zeros(rsds.shape[1:]) + tsolar = np.zeros(ncin.variables["t2m"][0].shape) + tp01m = np.zeros(ncin.variables["t2m"][0].shape) + tevap = np.zeros(ncin.variables["t2m"][0].shape) # J m-2 to W/m2 val = np.flipud(ncin.variables["ssrd"][0]) newval = (val - tsolar) / 3600.0 # Le Sigh, sometimes things are negative, somehow? - nc.variables["rsds"][tidx] = np.ma.where(newval < 0, 0, newval) + rsds[tidx, ijs, ijs] = np.ma.where(newval < 0, 0, newval) # m to mm val = np.flipud(ncin.variables["e"][0]) newval = (val * 1000.0) - tevap - nc.variables["evap"][tidx] = np.ma.where(newval < 0, 0, newval) + evap[tidx, ijs, ijs] = np.ma.where(newval < 0, 0, newval) # m to mm val = np.flipud(ncin.variables["tp"][0]) newval = (val * 1000.0) - tp01m - nc.variables["p01m"][tidx] = np.ma.where(newval < 0, 0, newval) + p01m[tidx, ijs, ijs] = np.ma.where(newval < 0, 0, newval) -def run(valid, domain, force): +def run(valid, domain: str, force): """Run for the given valid time.""" dd = "" if domain == "" else f"_{domain}" ncoutfn = f"/mesonet/data/era5{dd}/{valid.year}_era5land_hourly.nc" @@ -111,7 +139,7 @@ def run(valid, domain, force): domain, ) return - dom = iemre.DOMAINS[domain] + dom = era5land.DOMAINS[domain] LOG.info("Running for %s[domain=%s]", valid, domain) ncfn = f"{domain}_{valid:%Y%m%d%H}.nc" @@ -125,11 +153,13 @@ def run(valid, domain, force): "month": f"{valid.month}", "day": f"{valid.day}", "time": f"{valid:%H}:00", + # Unclear what is happening here, but the download seems to + # conform to the requested grid "area": [ - dom["north_edge"], - dom["west_edge"], - dom["south_edge"], - dom["east_edge"], + dom["NORTH"], + dom["WEST"], + dom["SOUTH"], + dom["EAST"], ], "format": "netcdf", }, diff --git a/scripts/era5/init_hourly.py b/scripts/era5/init_hourly.py index 4f8df6095..b9eabe65a 100644 --- a/scripts/era5/init_hourly.py +++ b/scripts/era5/init_hourly.py @@ -1,22 +1,22 @@ """Generate the ERA5 hourly analysis file for a year""" -import datetime import os +from datetime import datetime import click import numpy as np -from pyiem import iemre +from pyiem import era5land from pyiem.util import logger, ncopen LOG = logger() -def init_year(ts, domain): +def init_year(ts: datetime, domain: str): """ Create a new NetCDF file for a year of our specification! """ + dom = era5land.DOMAINS[domain] dd = "" if domain == "" else f"_{domain}" - dom = iemre.DOMAINS[domain] ncfn = f"/mesonet/data/era5{dd}/{ts.year}_era5land_hourly.nc" if os.path.isfile(ncfn): LOG.info("Cowardly refusing to overwrite: %s", ncfn) @@ -31,19 +31,13 @@ def init_year(ts, domain): nc.realization = 1 nc.Conventions = "CF-1.0" nc.contact = "Daryl Herzmann, akrherz@iastate.edu, 515-294-5978" - nc.history = f"{datetime.datetime.now():%d %B %Y} Generated" + nc.history = f"{datetime.now():%d %B %Y} Generated" nc.comment = "No Comment at this time" - # Setup Dimensions - nc.createDimension( - "lat", - (dom["north"] - dom["south"]) * 10.0 + 1, - ) - nc.createDimension( - "lon", - (dom["east"] - dom["west"]) * 10.0 + 1, - ) - ts2 = datetime.datetime(ts.year + 1, 1, 1) + nc.createDimension("lat", dom["NY"]) + nc.createDimension("lon", dom["NX"]) + nc.createDimension("nv", 2) + ts2 = datetime(ts.year + 1, 1, 1) days = (ts2 - ts).days LOG.info("Year %s has %s days", ts.year, days) nc.createDimension("time", int(days) * 24) @@ -60,14 +54,24 @@ def init_year(ts, domain): lat.long_name = "Latitude" lat.standard_name = "latitude" lat.axis = "Y" - lat[:] = np.arange(dom["south"], dom["north"] + 0.001, 0.1) + lat.bounds = "lat_bnds" + lat[:] = dom["YAXIS"] + + lat_bnds = nc.createVariable("lat_bnds", float, ("lat", "nv")) + lat_bnds[:, 0] = dom["YAXIS"] - 0.05 + lat_bnds[:, 1] = dom["YAXIS"] + 0.05 lon = nc.createVariable("lon", float, ("lon",)) lon.units = "degrees_east" lon.long_name = "Longitude" lon.standard_name = "longitude" lon.axis = "X" - lon[:] = np.arange(dom["west"], dom["east"] + 0.001, 0.1) + lon.bounds = "lon_bnds" + lon[:] = dom["XAXIS"] + + lon_bnds = nc.createVariable("lon_bnds", float, ("lon", "nv")) + lon_bnds[:, 0] = dom["XAXIS"] - 0.05 + lon_bnds[:, 1] = dom["XAXIS"] + 0.05 tm = nc.createVariable("time", float, ("time",)) tm.units = f"Hours since {ts.year}-01-01 00:00:0.0" @@ -188,7 +192,7 @@ def init_year(ts, domain): @click.option("--domain", default="", help="IEMRE Domain") def main(year, domain): """Go Main Go""" - init_year(datetime.datetime(year, 1, 1), domain) + init_year(datetime(year, 1, 1), domain) if __name__ == "__main__": diff --git a/scripts/gfs/gfs2iemre.py b/scripts/gfs/gfs2iemre.py index 12c7a6419..1e8aacd30 100644 --- a/scripts/gfs/gfs2iemre.py +++ b/scripts/gfs/gfs2iemre.py @@ -125,11 +125,6 @@ def create(ts: datetime, domain: str, dom: dict) -> str: return fn -def shuffle(grid): - """Convert the 0-360 grid to -180 to 180.""" - return np.concatenate((grid[:, 1536:], grid[:, :1536]), axis=1) - - def merge_grib(nc, now, domain: str, dom: dict): """Merge what grib data we can find into the netcdf file.""" # Tricky here for how to compute a valid date @@ -164,7 +159,7 @@ def merge_grib(nc, now, domain: str, dom: dict): affine = Affine( dx, 0.0, - -180, # not right, but close enough? + 0 - dx / 2.0, 0.0, -dx, grb["latitudeOfFirstGridPointInDegrees"] + dx / 2.0, @@ -217,20 +212,20 @@ def merge_grib(nc, now, domain: str, dom: dict): days = (approxlocal.date() - now.date()).days if hits == 4: nc.variables["high_tmpk"][days] = iemre.reproject2iemre( - shuffle(tmaxgrid), affine, EPSG[4326], domain=domain + tmaxgrid, affine, EPSG[4326], domain=domain ) nc.variables["low_tmpk"][days] = iemre.reproject2iemre( - shuffle(tmingrid), affine, EPSG[4326], domain=domain + tmingrid, affine, EPSG[4326], domain=domain ) nc.variables["p01d"][days] = iemre.reproject2iemre( - shuffle(pgrid), affine, EPSG[4326], domain=domain + pgrid, affine, EPSG[4326], domain=domain ) pout = nc.variables["p01d"][days] nc.variables["tsoil"][days] = iemre.reproject2iemre( - shuffle(tsoilgrid / 4.0), affine, EPSG[4326], domain=domain + tsoilgrid / 4.0, affine, EPSG[4326], domain=domain ) nc.variables["srad"][days] = iemre.reproject2iemre( - shuffle(srad), affine, EPSG[4326], domain=domain + srad, affine, EPSG[4326], domain=domain ) sout = nc.variables["srad"][days] LOG.info( diff --git a/scripts/hrrr/hrrr_ref2raster.py b/scripts/hrrr/hrrr_ref2raster.py index f16d89379..f61b6118a 100644 --- a/scripts/hrrr/hrrr_ref2raster.py +++ b/scripts/hrrr/hrrr_ref2raster.py @@ -67,6 +67,7 @@ def do_step( lat1 = ds.refd.GRIB_latitudeOfFirstGridPointInDegrees lon1 = ds.refd.GRIB_longitudeOfFirstGridPointInDegrees llx, lly = pyproj.Proj(projparams)(lon1, lat1) + # This is the edge and not the corner hrrr_aff = Affine( ds.refd.GRIB_DxInMetres, 0.0, diff --git a/scripts/iemre/grid_rsds.py b/scripts/iemre/grid_rsds.py index aebd4af3a..b12056a90 100644 --- a/scripts/iemre/grid_rsds.py +++ b/scripts/iemre/grid_rsds.py @@ -21,7 +21,7 @@ import pyproj import xarray as xr from affine import Affine -from pyiem import iemre +from pyiem import era5land, iemre from pyiem.grid.util import grid_smear from pyiem.util import archive_fetch, logger, ncopen, utc @@ -63,10 +63,8 @@ def try_era5land(ts: datetime, domain: str, dom: dict) -> Optional[np.ndarray]: total = total / 24.0 # The affine defines the edges of the grid - aff = Affine(0.1, 0, dom["west"], 0, -0.1, dom["north"]) - vals = iemre.reproject2iemre( - np.flipud(total), aff, P4326.crs, domain=domain - ) + aff = era5land.DOMAINS[domain]["AFFINE_NC"] + vals = iemre.reproject2iemre(total, aff, P4326.crs, domain=domain) # ERA5Land is too tight to the coast, so we need to jitter the grid # 7 and 4 was found to be enough to appease DEP's needs @@ -196,6 +194,7 @@ def do_hrrr(ts: datetime) -> Optional[np.ndarray]: # We wanna store as W m-2, so we just average out the data by hour total = total / 24.0 + # This is the edge and not the center of the UL pixel affine_in = Affine( dx, 0.0, llcrnrx - dx / 2.0, 0.0, dy, llcrnry + dy / 2.0 ) diff --git a/scripts/iemre/hourly_analysis.py b/scripts/iemre/hourly_analysis.py index 917f72b04..81e47c934 100644 --- a/scripts/iemre/hourly_analysis.py +++ b/scripts/iemre/hourly_analysis.py @@ -12,11 +12,10 @@ import pandas as pd import pint import pygrib -from affine import Affine from metpy.calc import wind_components from metpy.interpolate import inverse_distance_to_grid from metpy.units import masked_array, units -from pyiem import iemre +from pyiem import era5land, iemre from pyiem.database import get_sqlalchemy_conn from pyiem.grid.util import grid_smear from pyiem.iemre import grb2iemre, hourly_offset, reproject2iemre @@ -40,21 +39,9 @@ def use_era5land(ts, kind, domain): LOG.warning("Failed to find %s", ncfn) return None tidx = hourly_offset(ts) + aff = era5land.DOMAINS[domain]["AFFINE_NC"] try: with ncopen(ncfn) as nc: - lats = nc.variables["lat"][:] - lons = nc.variables["lon"][:] - dx = lons[1] - lons[0] - dy = lats[1] - lats[0] - # This defines the SW edge of the grid - aff = Affine( - dx, - 0, - lons[0] - dx / 2.0, - 0, - dy, - lats[0] - dy / 2.0, - ) res = [] for task in tasks.get(kind, [kind]): if task == "soilt": diff --git a/scripts/iemre/init_stage4_daily.py b/scripts/iemre/init_stage4_daily.py index 92cd13d6e..ca8a01260 100644 --- a/scripts/iemre/init_stage4_daily.py +++ b/scripts/iemre/init_stage4_daily.py @@ -6,6 +6,7 @@ import click import numpy as np +from pyiem import stage4 from pyiem.util import logger, ncopen LOG = logger() @@ -33,14 +34,41 @@ def init_year(ts: datetime, ci: bool) -> Optional[str]: nc.contact = "Daryl Herzmann, akrherz@iastate.edu, 515-294-5978" nc.history = "%s Generated" % (datetime.now().strftime("%d %B %Y"),) nc.comment = "No Comment at this time" + # Store projection information + nc.proj4 = stage4.PROJ.srs # Setup Dimensions nc.createDimension("x", tmplnc.dimensions["x"].size) nc.createDimension("y", tmplnc.dimensions["y"].size) + nc.createDimension("nv", 2) ts2 = datetime(ts.year + 1, 1, 1) days = 2 if ci else (ts2 - ts).days nc.createDimension("time", int(days)) + x = nc.createVariable("x", float, ("x",)) + x.units = "m" + x.long_name = "X coordinate of projection" + x.standard_name = "projection_x_coordinate" + x.axis = "X" + x.bounds = "x_bounds" + x[:] = stage4.XAXIS + + x_bounds = nc.createVariable("x_bounds", float, ("x", "bnds")) + x_bounds[:, 0] = stage4.XAXIS - stage4.DX / 2.0 + x_bounds[:, 1] = stage4.XAXIS + stage4.DX / 2.0 + + y = nc.createVariable("y", float, ("y",)) + y.units = "m" + y.long_name = "Y coordinate of projection" + y.standard_name = "projection_y_coordinate" + y.axis = "Y" + y.bounds = "y_bounds" + y[:] = stage4.YAXIS + + y_bounds = nc.createVariable("y_bounds", float, ("y", "bnds")) + y_bounds[:, 0] = stage4.YAXIS - stage4.DY / 2.0 + y_bounds[:, 1] = stage4.YAXIS + stage4.DY / 2.0 + # Setup Coordinate Variables lat = nc.createVariable("lat", float, ("y", "x")) lat.units = "degrees_north" diff --git a/scripts/iemre/init_stage4_dailyc.py b/scripts/iemre/init_stage4_dailyc.py index f86316c2d..55a2d870e 100644 --- a/scripts/iemre/init_stage4_dailyc.py +++ b/scripts/iemre/init_stage4_dailyc.py @@ -4,6 +4,7 @@ import os import numpy as np +from pyiem import stage4 from pyiem.util import logger, ncopen LOG = logger() @@ -33,14 +34,41 @@ def init_year(ts): datetime.datetime.now().strftime("%d %B %Y"), ) nc.comment = "No Comment at this time" + # Store projection information + nc.proj4 = stage4.PROJ.srs # Setup Dimensions nc.createDimension("x", tmplnc.dimensions["x"].size) nc.createDimension("y", tmplnc.dimensions["y"].size) + nc.createDimension("nv", 2) ts2 = datetime.datetime(ts.year + 1, 1, 1) days = (ts2 - ts).days nc.createDimension("time", int(days)) + x = nc.createVariable("x", float, ("x",)) + x.units = "m" + x.long_name = "X coordinate of projection" + x.standard_name = "projection_x_coordinate" + x.axis = "X" + x.bounds = "x_bounds" + x[:] = stage4.XAXIS + + x_bounds = nc.createVariable("x_bounds", float, ("x", "bnds")) + x_bounds[:, 0] = stage4.XAXIS - stage4.DX / 2.0 + x_bounds[:, 1] = stage4.XAXIS + stage4.DX / 2.0 + + y = nc.createVariable("y", float, ("y",)) + y.units = "m" + y.long_name = "Y coordinate of projection" + y.standard_name = "projection_y_coordinate" + y.axis = "Y" + y.bounds = "y_bounds" + y[:] = stage4.YAXIS + + y_bounds = nc.createVariable("y_bounds", float, ("y", "bnds")) + y_bounds[:, 0] = stage4.YAXIS - stage4.DY / 2.0 + y_bounds[:, 1] = stage4.YAXIS + stage4.DY / 2.0 + # Setup Coordinate Variables lat = nc.createVariable("lat", float, ("y", "x")) lat.units = "degrees_north" diff --git a/scripts/iemre/init_stage4_hourly.py b/scripts/iemre/init_stage4_hourly.py index 15fab1f25..4da4a5b14 100644 --- a/scripts/iemre/init_stage4_hourly.py +++ b/scripts/iemre/init_stage4_hourly.py @@ -7,6 +7,7 @@ import click import numpy as np import pygrib +from pyiem import stage4 from pyiem.util import archive_fetch, logger, ncopen BASEDIR = "/mesonet/data/stage4" @@ -41,7 +42,8 @@ def init_year(ts: datetime, ci: bool) -> Optional[str]: nc.contact = "Daryl Herzmann, akrherz@iastate.edu, 515-294-5978" nc.history = ("%s Generated") % (datetime.now().strftime("%d %B %Y"),) nc.comment = "No Comment at this time" - + # Store projection information + nc.proj4 = stage4.PROJ.srs # Setup Dimensions nc.createDimension("x", lats.shape[1]) nc.createDimension("y", lats.shape[0]) @@ -51,6 +53,30 @@ def init_year(ts: datetime, ci: bool) -> Optional[str]: LOG.info("Year %s has %s days", ts.year, days) nc.createDimension("time", int(days) * 24) + x = nc.createVariable("x", float, ("x",)) + x.units = "m" + x.long_name = "X coordinate of projection" + x.standard_name = "projection_x_coordinate" + x.axis = "X" + x.bounds = "x_bounds" + x[:] = stage4.XAXIS + + x_bounds = nc.createVariable("x_bounds", float, ("x", "bnds")) + x_bounds[:, 0] = stage4.XAXIS - stage4.DX / 2.0 + x_bounds[:, 1] = stage4.XAXIS + stage4.DX / 2.0 + + y = nc.createVariable("y", float, ("y",)) + y.units = "m" + y.long_name = "Y coordinate of projection" + y.standard_name = "projection_y_coordinate" + y.axis = "Y" + y.bounds = "y_bounds" + y[:] = stage4.YAXIS + + y_bounds = nc.createVariable("y_bounds", float, ("y", "bnds")) + y_bounds[:, 0] = stage4.YAXIS - stage4.DY / 2.0 + y_bounds[:, 1] = stage4.YAXIS + stage4.DY / 2.0 + # Setup Coordinate Variables lat = nc.createVariable("lat", float, ("y", "x")) lat.units = "degrees_north" diff --git a/scripts/iemre/precip_ingest.py b/scripts/iemre/precip_ingest.py index 1b5f732c0..4f45f5f6d 100644 --- a/scripts/iemre/precip_ingest.py +++ b/scripts/iemre/precip_ingest.py @@ -11,8 +11,7 @@ import click import numpy as np import pygrib -from affine import Affine -from pyiem import iemre +from pyiem import era5land, iemre, stage4 from pyiem.stage4 import AFFINE_NATIVE as STAGE4_AFFINE from pyiem.stage4 import AFFINE_NATIVE_OLD as STAGE4_AFFINE_OLD from pyiem.stage4 import ARCHIVE_FLIP as STAGE4_ARCHIVE_FLIP @@ -66,7 +65,7 @@ def ingest_hourly_grib(valid): else: p01m[tidx, :, :] = val nc.variables["p01m_status"][tidx] = 1 - LOG.debug( + LOG.info( "write p01m to stage4 netcdf min: %.2f avg: %.2f max: %.2f", np.min(val), np.mean(val), @@ -83,18 +82,9 @@ def copy_to_iemre(valid): val = nc.variables["p01m"][tidx] LOG.info("stage4 mean: %.2f max: %.2f", np.mean(val), np.max(val)) - # Define stage4 projection and affine - projparams = { - "a": 6371200.0, - "b": 6371200.0, - "proj": "stere", - "lat_ts": 60.0, - "lat_0": 90.0, - "lon_0": 255.0 - 360.0, - } # Reproject to IEMRE aff = STAGE4_AFFINE if valid >= STAGE4_ARCHIVE_FLIP else STAGE4_AFFINE_OLD - res = iemre.reproject2iemre(val, aff, projparams, domain="") + res = iemre.reproject2iemre(val, aff, stage4.PROJPARMS, domain="") LOG.info("iemre mean: %.2f max: %.2f", np.mean(res), np.max(res)) # Lets clip bad data @@ -123,9 +113,8 @@ def era5workflow(valid): p01m = nc.variables["p01m"][idx] # Convert trace/drizzle to 0, values < 0.01in or .254mm p01m[p01m < 0.254] = 0 - dom = iemre.DOMAINS[""] - affine_in = Affine(0.1, 0, dom["west"], 0, -0.1, dom["north"]) - val = iemre.reproject2iemre(np.flipud(p01m), affine_in, "EPSG:4326") + affine_in = era5land.DOMAINS[""]["AFFINE_NC"] + val = iemre.reproject2iemre(p01m, affine_in, "EPSG:4326") with ncopen( iemre.get_hourly_ncname(valid.year, domain=""), "a", timeout=300 ) as nc: diff --git a/scripts/iemre/prism_adjust_stage4.py b/scripts/iemre/prism_adjust_stage4.py index b2a5a8b12..7683f5718 100644 --- a/scripts/iemre/prism_adjust_stage4.py +++ b/scripts/iemre/prism_adjust_stage4.py @@ -5,30 +5,31 @@ RUN_NOON.sh for 1 day ago """ -import datetime import os +from datetime import datetime, timedelta import click import numpy as np import pandas as pd from pyiem import prism as prismutil +from pyiem import stage4 as stage4util from pyiem.iemre import daily_offset, hourly_offset -from pyiem.util import find_ij, logger, ncopen, utc -from scipy.interpolate import NearestNDInterpolator +from pyiem.util import logger, ncopen, utc +from rasterio.warp import reproject DEBUGLON = -93.89 DEBUGLAT = 42.04 LOG = logger() -def compute_s4total(valid): +def compute_s4total(valid: datetime): """Figure out the 24 hour total for this timestamp.""" # Split into two, so to support 1 January without yucky code # 13z yesterday, arrears - sts = valid - datetime.timedelta(hours=23) + sts = valid - timedelta(hours=23) tidx0 = hourly_offset(sts) # 23z yesterday - tidx1 = hourly_offset(sts + datetime.timedelta(hours=10)) + tidx1 = hourly_offset(sts + timedelta(hours=10)) ncfn = f"/mesonet/data/stage4/{sts.year}_stage4_hourly.nc" stage4total = None if os.path.isfile(ncfn): @@ -38,7 +39,7 @@ def compute_s4total(valid): stage4total = np.sum(p01m[tidx0 : (tidx1 + 1), :, :], axis=0) # 0z today - sts = valid - datetime.timedelta(hours=12) + sts = valid - timedelta(hours=12) tidx0 = hourly_offset(sts) # 11z today tidx1 = hourly_offset(valid) @@ -52,56 +53,50 @@ def compute_s4total(valid): return stage4total -def workflow(valid): +def workflow(valid: datetime): """Our workflow""" LOG.info("Processing %s", valid) # read prism tidx = daily_offset(valid) with ncopen(f"/mesonet/data/prism/{valid.year}_daily.nc", "r") as nc: - ppt = nc.variables["ppt"][tidx, :, :] - # missing as zero - lons = nc.variables["lon"][:] - lats = nc.variables["lat"][:] - ppt = np.where(ppt.mask, 0, ppt) - (lons, lats) = np.meshgrid(lons, lats) + # rasterio freaks out if we have masked arrays + ppt = nc.variables["ppt"][tidx].filled(0) (pi, pj) = prismutil.find_ij(DEBUGLON, DEBUGLAT) - LOG.info("prism debug point ppt: %.3f", ppt[pj, pi]) s4total = compute_s4total(valid) - with ncopen(f"/mesonet/data/stage4/{valid.year}_stage4_hourly.nc") as nc: - s4lons = nc.variables["lon"][:] - s4lats = nc.variables["lat"][:] - sj, si = find_ij(s4lons, s4lats, DEBUGLON, DEBUGLAT) - LOG.info( - "stage4 s4total: %.3f lon: %.2f (%.2f) lat: %.2f (%.2f)", - s4total[sj, si], - s4lons[sj, si], - DEBUGLON, - s4lats[sj, si], - DEBUGLAT, - ) + (si, sj) = stage4util.find_ij(DEBUGLON, DEBUGLAT) # make sure the s4total does not have zeros s4total = np.where(s4total < 0.001, 0.001, s4total) - nn = NearestNDInterpolator((lons.flatten(), lats.flatten()), ppt.flat) - prism_on_s4grid = nn(s4lons, s4lats) - multiplier = prism_on_s4grid / s4total - LOG.info( - "gridavgs: prism: %.3f stageIV: %.3f prismons4grid: %.3f mul: %.3f", - np.mean(ppt), - np.mean(s4total), - np.mean(prism_on_s4grid), - np.mean(multiplier), + # reproject prism onto the stage4 grid + prism_on_s4grid = np.zeros_like(s4total) + reproject( + ppt, + prism_on_s4grid, + src_transform=prismutil.AFFINE_NC, + src_crs="EPSG:4326", + dst_transform=stage4util.AFFINE_NC, + dst_crs=stage4util.PROJPARMS, + dst_nodata=0, ) + + multiplier = prism_on_s4grid / s4total LOG.info( - "prism: %.3f stageIV: %.4f prismons4grid: %.3f mul: %.3f", + "DEBUGPT prism: %.3f s4total: %.3f prism_on_s4grid: %.3f mul: %.3f", ppt[pj, pi], s4total[sj, si], prism_on_s4grid[sj, si], multiplier[sj, si], ) + LOG.info( + "gridavgs: prism: %.3f stageIV: %.3f prism_on_s4grid: %.3f mul: %.3f", + np.mean(ppt), + np.mean(s4total), + np.mean(prism_on_s4grid), + np.mean(multiplier), + ) # 13 z yesterday - sts = valid - datetime.timedelta(hours=23) + sts = valid - timedelta(hours=23) # Through 12z file ets = valid for now in pd.date_range(sts, ets, freq="1h"): @@ -119,7 +114,7 @@ def workflow(valid): newval = np.where(newval < 0.01, 0, newval) nc.variables["p01m"][idx, :, :] = newval LOG.info( - "adjust %s[%s] oldval: %.3f newval: %.3f", + "adjust %s[%s] DEBUGPT: oldval: %.3f newval: %.3f", now.strftime("%Y%m%d%H"), idx, oldval[sj, si], @@ -130,10 +125,9 @@ def workflow(valid): nc.variables["p01m_status"][idx] = 3 else: LOG.warning( - "NOOP for time %s[idx:%s]", + "Not setting p01m_status due to zeroes for %s[idx:%s]", ( - datetime.datetime(valid.year, 1, 1, 0) - + datetime.timedelta(hours=idx) + datetime(valid.year, 1, 1, 0) + timedelta(hours=idx) ).strftime("%Y-%m-%dT%H"), idx, ) @@ -141,7 +135,7 @@ def workflow(valid): @click.command() @click.option("--date", "dt", type=click.DateTime(), required=True) -def main(dt: datetime.datetime): +def main(dt: datetime): """Go Main""" workflow(utc(dt.year, dt.month, dt.day, 12)) diff --git a/scripts/iemre/stage4_12z_adjust.py b/scripts/iemre/stage4_12z_adjust.py index 2422dd60a..7fbd7ff81 100644 --- a/scripts/iemre/stage4_12z_adjust.py +++ b/scripts/iemre/stage4_12z_adjust.py @@ -45,10 +45,8 @@ def merge(ts): if fn is None: LOG.info("stage4_12z_adjust %s is missing", ppath) return - grbs = pygrib.open(fn) - grb = grbs[1] - val = grb.values - grbs.close() + with pygrib.open(fn) as grbs: + val = grbs[1].values save12z(ts, val) # storage is in the arrears @@ -64,17 +62,18 @@ def merge(ts): ncfn = f"/mesonet/data/stage4/{pair[0].year}_stage4_hourly.nc" idx0 = iemre.hourly_offset(pair[0]) idx1 = iemre.hourly_offset(pair[1]) + tslice = slice(idx0, idx1 + 1) LOG.info("%s [%s thru %s]", ncfn, idx0, idx1) with ncopen(ncfn, timeout=60) as nc: # Check that the status value is (1,2), otherwise prism ran sentinel = np.nanmax( - nc.variables["p01m_status"][idx0 : (idx1 + 1)] + nc.variables["p01m_status"][tslice], ) if sentinel == 3: - LOG.warning("Aborting as p01m_status[%s] == 3", sentinel) + LOG.warning("Aborting as p01m_status[%s] == 3", tslice) return p01m = nc.variables["p01m"] - total = np.nansum(p01m[idx0 : (idx1 + 1)], axis=0) + total = np.nansum(p01m[tslice], axis=0) if hourly_total is None: hourly_total = total else: diff --git a/scripts/mrms/init_daily_mrms.py b/scripts/mrms/init_daily_mrms.py index 11f11a561..2ed752670 100644 --- a/scripts/mrms/init_daily_mrms.py +++ b/scripts/mrms/init_daily_mrms.py @@ -5,13 +5,13 @@ import click import numpy as np -from pyiem import iemre +from pyiem import iemre, mrms from pyiem.util import logger, ncopen LOG = logger() -def init_year(ts, for_dep=False, ci=False): +def init_year(ts: datetime, for_dep=False, ci=False): """Create the netcdf file for a given year. Args: @@ -22,16 +22,6 @@ def init_year(ts, for_dep=False, ci=False): fn = iemre.get_daily_mrms_ncname(ts.year) os.makedirs(os.path.dirname(fn), exist_ok=True) - # Centroid of the upper left pixel - ulx = -125.995 - uly = 49.995 - ny = 2700 - nx = 6100 - dx = 0.01 - dy = 0.01 - # Centroid of the lower right pixel - lry = 23.005 - lrx = -65.005 if for_dep: fn = fn.replace("daily", "dep") if os.path.isfile(fn): @@ -51,8 +41,8 @@ def init_year(ts, for_dep=False, ci=False): nc.comment = "No Comment at this time" # Setup Dimensions - nc.createDimension("lat", ny) - nc.createDimension("lon", nx) + nc.createDimension("lat", mrms.MRMS4IEMRE_NX) + nc.createDimension("lon", mrms.MRMS4IEMRE_NY) days = 2 if ci else ((ts.replace(year=ts.year + 1)) - ts).days day_axis = np.arange(0, int(days)) if for_dep: @@ -71,11 +61,11 @@ def init_year(ts, for_dep=False, ci=False): lat.bounds = "lat_bnds" lat.axis = "Y" # Grid centers - lat[:] = np.arange(lry, uly + 0.0001, dy) + lat[:] = mrms.MRMS4IEMRE_YAXIS lat_bnds = nc.createVariable("lat_bnds", float, ("lat", "nv")) - lat_bnds[:, 0] = lat[:] - dy / 2.0 - lat_bnds[:, 1] = lat[:] + dy / 2.0 + lat_bnds[:, 0] = lat[:] - mrms.MRMS4IEMRE_DY / 2.0 + lat_bnds[:, 1] = lat[:] + mrms.MRMS4IEMRE_DY / 2.0 lon = nc.createVariable("lon", float, ("lon",)) lon.units = "degrees_east" @@ -83,11 +73,11 @@ def init_year(ts, for_dep=False, ci=False): lon.standard_name = "longitude" lon.bounds = "lon_bnds" lon.axis = "X" - lon[:] = np.arange(ulx, lrx + 0.0001, dx) + lon[:] = mrms.MRMS4IEMRE_XAXIS lon_bnds = nc.createVariable("lon_bnds", float, ("lon", "nv")) - lon_bnds[:, 0] = lon[:] - dx / 2.0 - lon_bnds[:, 1] = lon[:] + dx / 2.0 + lon_bnds[:, 0] = lon[:] - mrms.MRMS4IEMRE_DX / 2.0 + lon_bnds[:, 1] = lon[:] + mrms.MRMS4IEMRE_DX / 2.0 tm = nc.createVariable("time", float, ("time",)) tm.units = f"Days since {ts.year}-01-01 00:00:0.0" @@ -110,11 +100,13 @@ def init_year(ts, for_dep=False, ci=False): # Add some faked data if ci: - p01d[0] = np.meshgrid(np.linspace(0, 25, nx), np.linspace(0, 25, ny))[ - 0 - ] + p01d[0] = np.meshgrid( + np.linspace(0, 25, mrms.MRMS4IEMRE_NX), + np.linspace(0, 25, mrms.MRMS4IEMRE_NY), + )[0] p01d[1] = np.meshgrid( - np.linspace(0, 100, nx), np.linspace(0, 100, ny) + np.linspace(0, 100, mrms.MRMS4IEMRE_NX), + np.linspace(0, 100, mrms.MRMS4IEMRE_NY), )[0] nc.close() diff --git a/scripts/mrms/init_mrms_dailyc.py b/scripts/mrms/init_mrms_dailyc.py index 949b4df4d..219de7d0c 100644 --- a/scripts/mrms/init_mrms_dailyc.py +++ b/scripts/mrms/init_mrms_dailyc.py @@ -4,7 +4,7 @@ import os import numpy as np -from pyiem import iemre +from pyiem import iemre, mrms from pyiem.util import logger, ncopen LOG = logger() @@ -15,21 +15,15 @@ def init_year(ts): Create a new NetCDF file for a year of our specification! """ fn = iemre.get_dailyc_mrms_ncname() - # Centroid of the upper left pixel - ulx = -125.995 - uly = 49.995 - ny = 2700 - nx = 6100 - dx = 0.01 - dy = 0.01 - # Centroid of the lower right pixel - lry = 23.005 - lrx = -65.005 + ny = mrms.MRMS4IEMRE_NY + nx = mrms.MRMS4IEMRE_NX + dx = mrms.MRMS4IEMRE_DX + dy = mrms.MRMS4IEMRE_DY if os.path.isfile(fn): LOG.warning("Cowardly refusing to create fn: %s", fn) return nc = ncopen(fn, "w") - nc.title = "IEM Daily Reanalysis Climatology %s" % (ts.year,) + nc.title = f"IEM Daily Reanalysis Climatology {ts.year}" nc.platform = "Grided Climatology" nc.description = "IEM daily analysis on a 0.01 degree grid" nc.institution = "Iowa State University, Ames, IA, USA" @@ -59,7 +53,7 @@ def init_year(ts): lat.bounds = "lat_bnds" lat.axis = "Y" # Grid centers - lat[:] = np.arange(lry, uly + 0.0001, dy) + lat[:] = mrms.MRMS4IEMRE_YAXIS lat_bnds = nc.createVariable("lat_bnds", float, ("lat", "nv")) lat_bnds[:, 0] = lat[:] - dy / 2.0 @@ -71,7 +65,7 @@ def init_year(ts): lon.standard_name = "longitude" lon.bounds = "lon_bnds" lon.axis = "X" - lon[:] = np.arange(ulx, lrx + 0.0001, dx) + lon[:] = mrms.MRMS4IEMRE_XAXIS lon_bnds = nc.createVariable("lon_bnds", float, ("lon", "nv")) lon_bnds[:, 0] = lon[:] - dx / 2.0 diff --git a/scripts/mrms/merge_mrms_q3.py b/scripts/mrms/merge_mrms_q3.py index 19c8583ab..25659f08f 100644 --- a/scripts/mrms/merge_mrms_q3.py +++ b/scripts/mrms/merge_mrms_q3.py @@ -68,13 +68,13 @@ def run(dt: datetime.date, for_dep: bool): if utcdt < utcnow: LOG.info("MISSING %s", utcdt) continue - grbs = pygrib.open(gribfn) - grb = grbs[1] - if lats is None: - lats, _ = grb.latlons() + with pygrib.open(gribfn) as grbs: + grb = grbs[1] + if lats is None: + lats, _ = grb.latlons() + val = grb["values"] os.unlink(gribfn) - val = grb["values"] # Anything less than zero, we set to zero val = np.where(val < 0, 0, val) if total is None: diff --git a/scripts/swat/swat_realtime.py b/scripts/swat/swat_realtime.py index 6902c3996..18382b9fa 100644 --- a/scripts/swat/swat_realtime.py +++ b/scripts/swat/swat_realtime.py @@ -65,10 +65,8 @@ def workflow(page, huc12s): # default time domain, only run for one date if the files exist. ets = date.today() - timedelta(days=1) LOG.info("Running for time domain %s to %s", now, ets) - mrmsaffine = Affine(0.01, 0.0, mrms.WEST, 0.0, -0.01, mrms.NORTH) - mrms_czs = CachingZonalStats( - mrmsaffine, - ) + # This is N to S + mrms_czs = CachingZonalStats(mrms.AFFINE) dom = iemre.DOMAINS[""] iemreaffine = Affine( iemre.DX, 0.0, dom["west"], 0.0, 0 - iemre.DY, dom["north"] From f23d6b112f7fb0c6f09aef2ca7b3d4a4a2234d4b Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 19 Dec 2024 15:47:23 -0600 Subject: [PATCH 3/4] =?UTF-8?q?=F0=9F=93=9D=20Address=20lint=20and=20inter?= =?UTF-8?q?mediate=20bugs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- scripts/climodat/merra_solarrad.py | 4 ++-- scripts/iemre/grid_rsds.py | 4 ++-- scripts/mrms/init_daily_mrms.py | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/scripts/climodat/merra_solarrad.py b/scripts/climodat/merra_solarrad.py index fbf83a426..a33e958eb 100644 --- a/scripts/climodat/merra_solarrad.py +++ b/scripts/climodat/merra_solarrad.py @@ -155,8 +155,8 @@ def do(dt: datetime): cursor = pgconn.cursor() for sid, row in df[df[COL].notna()].iterrows(): cursor.execute( - f"UPDATE alldata set {COL} = %s where station = %s and " - "day = %s", + f"UPDATE alldata set {COL} = %s where station = %s " # skipcq + "and day = %s", (row[COL], sid, dt), ) diff --git a/scripts/iemre/grid_rsds.py b/scripts/iemre/grid_rsds.py index b12056a90..0f0ae0359 100644 --- a/scripts/iemre/grid_rsds.py +++ b/scripts/iemre/grid_rsds.py @@ -30,7 +30,7 @@ SWITCH_DATE = utc(2014, 10, 10, 20) -def try_era5land(ts: datetime, domain: str, dom: dict) -> Optional[np.ndarray]: +def try_era5land(ts: datetime, domain: str) -> Optional[np.ndarray]: """Attempt to use ERA5Land data.""" dd = "" if domain == "" else f"_{domain}" # inbound `ts` represents noon local time, we want values from 1 AM @@ -248,7 +248,7 @@ def main(dt: Optional[datetime], year, month): for sts_in in queue: for domain, dom in iemre.DOMAINS.items(): sts = sts_in.replace(tzinfo=dom["tzinfo"]) - srad = try_era5land(sts, domain, dom) + srad = try_era5land(sts, domain) if srad is not None and postprocess(srad, sts, domain): continue LOG.info("try_era5land failed to find data") diff --git a/scripts/mrms/init_daily_mrms.py b/scripts/mrms/init_daily_mrms.py index 2ed752670..7d1097e22 100644 --- a/scripts/mrms/init_daily_mrms.py +++ b/scripts/mrms/init_daily_mrms.py @@ -41,8 +41,8 @@ def init_year(ts: datetime, for_dep=False, ci=False): nc.comment = "No Comment at this time" # Setup Dimensions - nc.createDimension("lat", mrms.MRMS4IEMRE_NX) - nc.createDimension("lon", mrms.MRMS4IEMRE_NY) + nc.createDimension("lat", mrms.MRMS4IEMRE_NY) + nc.createDimension("lon", mrms.MRMS4IEMRE_NX) days = 2 if ci else ((ts.replace(year=ts.year + 1)) - ts).days day_axis = np.arange(0, int(days)) if for_dep: From 0f2b0699719d0368fa96c68eed74aeb58fc37075 Mon Sep 17 00:00:00 2001 From: akrherz Date: Thu, 19 Dec 2024 16:00:42 -0600 Subject: [PATCH 4/4] =?UTF-8?q?=F0=9F=90=9B=20Fix=20intermediate=20bug?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- scripts/iemre/init_stage4_daily.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/iemre/init_stage4_daily.py b/scripts/iemre/init_stage4_daily.py index ca8a01260..78528d6c6 100644 --- a/scripts/iemre/init_stage4_daily.py +++ b/scripts/iemre/init_stage4_daily.py @@ -53,7 +53,7 @@ def init_year(ts: datetime, ci: bool) -> Optional[str]: x.bounds = "x_bounds" x[:] = stage4.XAXIS - x_bounds = nc.createVariable("x_bounds", float, ("x", "bnds")) + x_bounds = nc.createVariable("x_bounds", float, ("x", "nv")) x_bounds[:, 0] = stage4.XAXIS - stage4.DX / 2.0 x_bounds[:, 1] = stage4.XAXIS + stage4.DX / 2.0 @@ -65,7 +65,7 @@ def init_year(ts: datetime, ci: bool) -> Optional[str]: y.bounds = "y_bounds" y[:] = stage4.YAXIS - y_bounds = nc.createVariable("y_bounds", float, ("y", "bnds")) + y_bounds = nc.createVariable("y_bounds", float, ("y", "nv")) y_bounds[:, 0] = stage4.YAXIS - stage4.DY / 2.0 y_bounds[:, 1] = stage4.YAXIS + stage4.DY / 2.0