From de3b4fff46d7d0c89dba38b8d4ca81242fd963a9 Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 9 Feb 2024 11:55:13 -0600 Subject: [PATCH 1/4] add min_rh,max_rh to IEMRE daily netcdf --- htdocs/plotting/auto/scripts/p86.py | 10 ++++++++-- scripts/iemre/db_to_netcdf.py | 1 + scripts/iemre/init_daily.py | 18 ++++++++++++++++++ 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/htdocs/plotting/auto/scripts/p86.py b/htdocs/plotting/auto/scripts/p86.py index 1845c35232..0dd7fc40ee 100644 --- a/htdocs/plotting/auto/scripts/p86.py +++ b/htdocs/plotting/auto/scripts/p86.py @@ -28,6 +28,8 @@ "high_tmpk_12z": "Maximum Temperature at 12 UTC", "high_soil4t": "Maximum 4 Inch Soil Temperature", "low_soil4t": "Minimum 4 Inch Soil Temperature", + "min_rh": "Minimum Relative Humidity", + "max_rh": "Maximum Relative Humidity", "power_swdn": "NASA POWER :: Incident Shortwave Down", "rsds": "Solar Radiation", "avg_dwpk": "Average Dew Point", @@ -80,7 +82,9 @@ def unit_convert(nc, varname, idx0, jslice, islice): data = None if not varname.startswith("range"): data = nc.variables[varname][idx0, jslice, islice] - if varname in ["rsds", "power_swdn"]: + if varname in ["min_rh", "max_rh"]: + pass + elif varname in ["rsds", "power_swdn"]: # Value is in W m**-2, we want MJ multi = (86400.0 / 1000000.0) if varname == "rsds" else 1 data = data * multi @@ -105,7 +109,6 @@ def unit_convert(nc, varname, idx0, jslice, islice): "high_soil4t", "low_soil4t", ]: - # Value is in W m**-2, we want MJ data = masked_array(data, units("degK")).to(units("degF")).m else: # range_tmpk range_tmpk_12z vname2 = f"low_tmpk{'_12z' if varname == 'range_tmpk_12z' else ''}" @@ -161,6 +164,9 @@ def plotter(fdict): plot_units = "mph" clevs = pretty_bins(0, ptiles[1]) clevs[0] = 0.01 + elif varname in ["min_rh", "max_rh"]: + plot_units = "%" + clevs = pretty_bins(0, 100) elif varname in ["p01d", "p01d_12z", "snow_12z", "snowd_12z"]: # Value is in W m**-2, we want MJ plot_units = "inch" diff --git a/scripts/iemre/db_to_netcdf.py b/scripts/iemre/db_to_netcdf.py index 2be43e4a3a..72b108c38b 100644 --- a/scripts/iemre/db_to_netcdf.py +++ b/scripts/iemre/db_to_netcdf.py @@ -33,6 +33,7 @@ def main(argv): with ncopen(ncfn, "a", timeout=600) as nc: for vname in ds: if vname not in nc.variables: + LOG.warning("Variable %s not in netcdf file, skipping", vname) continue # Careful here, ds could contain NaN values nc.variables[vname][idx] = np.ma.array( diff --git a/scripts/iemre/init_daily.py b/scripts/iemre/init_daily.py index d3e2eea18f..96a1837ca2 100644 --- a/scripts/iemre/init_daily.py +++ b/scripts/iemre/init_daily.py @@ -207,6 +207,24 @@ def init_year(ts): v1.standard_name = "4inch Soil Temperature" v1.coordinates = "lon lat" + v1 = nc.createVariable( + "min_rh", np.ubyte, ("time", "lat", "lon"), fill_value=255 + ) + v1.units = "%" + v1.scale_factor = 0.5 + v1.long_name = "Minimum Relative Humidity" + v1.standard_name = "Minimum Relative Humidity" + v1.coordinates = "lon lat" + + v1 = nc.createVariable( + "max_rh", np.ubyte, ("time", "lat", "lon"), fill_value=255 + ) + v1.units = "%" + v1.scale_factor = 0.5 + v1.long_name = "Maximum Relative Humidity" + v1.standard_name = "Maximum Relative Humidity" + v1.coordinates = "lon lat" + nc.close() From f5529527efd50a97a94afc4efd3c8dd679ed3d5f Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 9 Feb 2024 12:01:24 -0600 Subject: [PATCH 2/4] t(h)rash IEMRE code some more --- htdocs/iemre/index.phtml | 101 +++++++--- scripts/RUN_0Z.sh | 2 +- scripts/RUN_10_AFTER.sh | 2 +- scripts/RUN_MIDNIGHT.sh | 2 +- scripts/RUN_NOON.sh | 6 +- scripts/iemre/daily_analysis.py | 311 ++++++++++++------------------- scripts/iemre/hourly_analysis.py | 92 +++++---- 7 files changed, 237 insertions(+), 279 deletions(-) diff --git a/htdocs/iemre/index.phtml b/htdocs/iemre/index.phtml index 330bbdad58..a83b51a44e 100644 --- a/htdocs/iemre/index.phtml +++ b/htdocs/iemre/index.phtml @@ -115,38 +115,77 @@ degree. - +
-

Data Availability

- - - - - - - - +

Hourly Data Sources

+ +

Yearly netcdf files since 1950, with updates +at about 10 after the hour. Since we want real-time data, some of these +model sources are not available immediately, so the denoted intermediate source +fills the gap until the final source is available.

+ +
Hourly IntervalDaily Interval
-
    -
  • 2 Meter Air and Dew Point Temperature
    • -
  • 10 Meter Wind Speed and Direction
    • -
  • Precipitation
    • -
  • Approx 4inch Soil Temperature
    • -
  • Sky Coverage Percentage
    • -
-		 -
    -
  • High and Low 2 Meter Air Temperature
    • -
  • High and Low ~4 inch soil temperature
    • -
  • Average wind speed magnitude
    • -
  • Precipitation
    • -
  • Snowfall at 12 UTC
    • -
  • Snow Depth at 12 UTC
    • -
  • Solar Radiation
    • -
  • Solar Radiation via NASA POWER
    • -
  • Average Dew Point Temperature
    • -
-
+ + + + + + + + + + + + + + + + + + + + + + + + + + + +
VariableIntermediate SourceFinal Source
PrecipitationStage IV1997+ Stage IV adjusted by PRISM +
1950-1997 ERA5-Land
2m Air Temperature +
2 m Dew Point Temperature +
10m Wind Speed +		Crude ASOS grid analysis2010+: RTMA, once available +
1950-2010 ERA5-Land
4 inch (cough 0-10 or 0-7cm) Soil TemperatureHRRR1950+: ERA5-Land
Sky CoverageCrude ASOS grid analysisCrude ASOS grid analysis
+ + +

Daily Data Sources

+ +

Yearly netcdf files since 1893. Some of these +variables are directly sourced from the hourly grids above. The concept of +a daily variable is a bit nebulous, but is intended to be a 6 to 6 UTC period. +

+ + + + + + + + + + + + + + + + + +
VariableIntermediate SourceFinal Source
High/Low 4 inch (cough 0-10 or 0-7cm) Soil TemperatureHRRR1950+: ERA5-Land
Average 2m Dew Point TemperatureSimple average of hourly RTMA2010+: RTMA, once available +
1950-2010 ERA5-Land

Data Access

@@ -177,7 +216,7 @@ Example: https://mesonet.agron.iastate.edu/iemre/hourly/2010-05-01/42.54/-96.40/

Raw Data Access

- +

The IEM makes extensive use of the NetCDF file format to store these gridded datasets. You can find individual web directories for various IEMRE related datasets: diff --git a/scripts/RUN_0Z.sh b/scripts/RUN_0Z.sh index 7ea9ebedbb..8aea7d9903 100644 --- a/scripts/RUN_0Z.sh +++ b/scripts/RUN_0Z.sh @@ -13,7 +13,7 @@ python ingest_climdiv.py & cd ../iemre # need to run daily analysis for climodat estimator to then work -python daily_analysis.py $(date +'%Y %m %d') +python daily_analysis.py --date=$(date +'%Y-%m-%d') cd ../climodat python sync_coop_updates.py diff --git a/scripts/RUN_10_AFTER.sh b/scripts/RUN_10_AFTER.sh index 5b83bac3ae..e8a4fb0e38 100644 --- a/scripts/RUN_10_AFTER.sh +++ b/scripts/RUN_10_AFTER.sh @@ -35,7 +35,7 @@ fi # we run at 11 PM for today as well if [ $LHH -eq "23" ] then - python daily_analysis.py $TODAY + python daily_analysis.py --date=$(date +'%Y-%m-%d') python grid_rsds.py $TODAY fi diff --git a/scripts/RUN_MIDNIGHT.sh b/scripts/RUN_MIDNIGHT.sh index 15c36765a8..ceb78105df 100644 --- a/scripts/RUN_MIDNIGHT.sh +++ b/scripts/RUN_MIDNIGHT.sh @@ -12,7 +12,7 @@ python update_daily_srad.py $(date --date '1 day ago' +'%Y %m %d') # Need this done so that IEMRE daily grids are there for DEP cd ../iemre -python daily_analysis.py $(date --date '1 day ago' +'%Y %m %d') +python daily_analysis.py --date=$(date --date '1 day ago' +'%Y-%m-%d') cd ../asos python adjust_report_type.py $(date -u --date '1 day ago' +'%Y %m %d') diff --git a/scripts/RUN_NOON.sh b/scripts/RUN_NOON.sh index a77c43a105..be2d2326ed 100644 --- a/scripts/RUN_NOON.sh +++ b/scripts/RUN_NOON.sh @@ -9,7 +9,7 @@ cd ../iemre python stage4_12z_adjust.py $(date +'%Y %m %d') # Copy the data to IEMRE hourly python precip_ingest.py $(date +'%Y %m %d') -python daily_analysis.py $(date +'%Y %m %d') +python daily_analysis.py --date=$(date +'%Y-%m-%d') cd ../prism python ingest_prism.py $(date --date '1 days ago' +'%Y %m %d') @@ -34,9 +34,9 @@ python sync_coop_updates.py cd ../iemre # Since we have now adjusted the 12z precip 1 day ago, we should rerun # iemre for two days ago -python daily_analysis.py $(date --date '2 days ago' +'%Y %m %d') +python daily_analysis.py --date=$(date --date '2 days ago' +'%Y-%m-%d') # Updated soil temperature data from ERA5 -python daily_analysis.py $(date --date '10 days ago' +'%Y %m %d') +python daily_analysis.py --date=$(date --date '10 days ago' +'%Y-%m-%d') # and now recompute climodat statewide/climate from two days ago cd ../climodat diff --git a/scripts/iemre/daily_analysis.py b/scripts/iemre/daily_analysis.py index 4c56340fe0..a50a6aaae1 100644 --- a/scripts/iemre/daily_analysis.py +++ b/scripts/iemre/daily_analysis.py @@ -8,17 +8,18 @@ RUN_10_AFTER.sh """ import datetime -import os import subprocess -import sys +import click import numpy as np import pandas as pd +from metpy.calc import relative_humidity_from_dewpoint from metpy.interpolate import inverse_distance_to_grid +from metpy.units import units from pyiem import iemre +from pyiem.database import get_sqlalchemy_conn from pyiem.util import ( convert_value, - get_sqlalchemy_conn, logger, ncopen, utc, @@ -83,7 +84,7 @@ def generic_gridder(df, idx): def copy_iemre_hourly(ts, ds): - """Compute the 6 UTC to 6 UTC totals via IEMRE hourly values.""" + """Lots of work to do here...""" sts = utc(ts.year, ts.month, ts.day, 6) ets = sts + datetime.timedelta(hours=23) pairs = [(sts, ets)] @@ -94,110 +95,105 @@ def copy_iemre_hourly(ts, ds): (sts, sts + datetime.timedelta(hours=17)), (sts + datetime.timedelta(hours=18), ets), ] - windhours = 0 - max_tmpk = None - min_tmpk = None - hi_soil4t = None - lo_soil4t = None - sped = None - totalprecip = None - for pair in pairs: - ncfn = iemre.get_hourly_ncname(pair[0].year) - if not os.path.isfile(ncfn): - LOG.warning("Missing %s", ncfn) - continue - tidx1 = iemre.hourly_offset(pair[0]) - tidx2 = iemre.hourly_offset(pair[1]) - LOG.info("Using %s for [%s thru %s]", ncfn, tidx1, tidx2) - with ncopen(ncfn, timeout=600) as nc: - nc_p01m = nc.variables["p01m"] - nc_uwnd = nc.variables["uwnd"] - nc_vwnd = nc.variables["vwnd"] - nc_t4 = nc.variables["soil4t"] - nc_tmpk = nc.variables["tmpk"] - # Caution: precip is stored arrears, so we have ugliness - _idx1 = 0 if tidx1 == 0 else tidx1 + 1 - LOG.info("precip read [%s<%s]", _idx1, tidx2 + 2) - precip = np.nansum(nc_p01m[_idx1 : (tidx2 + 2), :, :], axis=0) - if totalprecip is None: - totalprecip = precip - else: - totalprecip = np.nansum([precip, totalprecip], axis=0) - maxt4 = np.ma.max(nc_t4[tidx1 : (tidx2 + 1), :, :], axis=0) - if hi_soil4t is None: - hi_soil4t = maxt4 - else: - hi_soil4t = np.ma.maximum(hi_soil4t, maxt4) - mint4 = np.ma.min(nc_t4[tidx1 : (tidx2 + 1), :, :], axis=0) - if lo_soil4t is None: - lo_soil4t = mint4 - else: - lo_soil4t = np.ma.maximum(lo_soil4t, mint4) - maxt = np.ma.max(nc_tmpk[tidx1 : (tidx2 + 1), :, :], axis=0) - if max_tmpk is None: - max_tmpk = maxt - else: - max_tmpk = np.ma.maximum(max_tmpk, maxt) - mint = np.ma.min(nc_tmpk[tidx1 : (tidx2 + 1), :, :], axis=0) - if min_tmpk is None: - min_tmpk = mint - else: - min_tmpk = np.ma.minimum(min_tmpk, mint) - for offset in range(tidx1, tidx2 + 1): - uwnd = nc_uwnd[offset, :, :] - vwnd = nc_vwnd[offset, :, :] - if uwnd.mask.all(): - LOG.info("No wind for offset: %s", offset) - continue - mag = (uwnd**2 + vwnd**2) ** 0.5 - windhours += 1 - if sped is None: - sped = mag - else: - sped += mag - if hi_soil4t is not None: - ds["high_soil4t"].values = hi_soil4t - ds["low_soil4t"].values = lo_soil4t - # NB: these are crude bias offsets involved when using hourly data - ds["high_tmpk"].values = max_tmpk + 0.8 - ds["low_tmpk"].values = min_tmpk - 0.8 - ds["p01d"].values = totalprecip - if windhours > 0: - ds["wind_speed"].values = sped / windhours - - -def copy_iemre_12z(ts, ds): - """Compute the 24 Hour precip at 12 UTC.""" - # arrears, so inclusive end ets = utc(ts.year, ts.month, ts.day, 12) # 13z yesterday sts = ets - datetime.timedelta(hours=23) - pairs = [(sts, ets)] + pairs12z = [(sts, ets)] if sts.year != ets.year: # 13z to 23z (inclusve) # 0z to 12z (inclusive) - pairs = [ + pairs12z = [ (sts, sts + datetime.timedelta(hours=10)), (sts + datetime.timedelta(hours=11), ets), ] - totalprecip = None - for pair in pairs: - ncfn = iemre.get_hourly_ncname(pair[0].year) - if not os.path.isfile(ncfn): - LOG.warning("Missing %s", ncfn) - continue - tidx1 = iemre.hourly_offset(pair[0]) - tidx2 = iemre.hourly_offset(pair[1]) - LOG.info("Using %s for [%s thru %s]", ncfn, tidx1, tidx2) - with ncopen(ncfn, timeout=600) as nc: - nc_p01m = nc.variables["p01m"] - # Caution: precip is stored arrears, so we have ugliness - precip = np.nansum(nc_p01m[tidx1 : (tidx2 + 1), :, :], axis=0) - if totalprecip is None: - totalprecip = precip - else: - totalprecip = np.nansum([precip, totalprecip], axis=0) - ds["p01d_12z"].values = totalprecip + + # One Off + for vname in ["min_rh", "max_rh"]: + aggfunc = np.ma.max if vname == "max_rh" else np.ma.min + res = None + for pair in pairs: + ncfn = iemre.get_hourly_ncname(pair[0].year) + tidx1 = iemre.hourly_offset(pair[0]) + tidx2 = iemre.hourly_offset(pair[1]) + with ncopen(ncfn, timeout=600) as nc: + tmpk = nc.variables["tmpk"] + dwpk = nc.variables["dwpk"] + for offset in range(tidx1, tidx2 + 1): + rh = ( + relative_humidity_from_dewpoint( + units("degK") * tmpk[offset], + units("degK") * dwpk[offset], + ).m + * 100.0 + ) + if res is None: + res = rh + else: + res = aggfunc([res, rh], axis=0) + ds[vname].values = res + + # One off + for vname in ["wind_speed"]: + hours = 0 + runningsum = None + for pair in pairs: + ncfn = iemre.get_hourly_ncname(pair[0].year) + tidx1 = iemre.hourly_offset(pair[0]) + tidx2 = iemre.hourly_offset(pair[1]) + with ncopen(ncfn, timeout=600) as nc: + uwnd = nc.variables["uwnd"] + vwnd = nc.variables["vwnd"] + for offset in range(tidx1, tidx2 + 1): + val = (uwnd[offset] ** 2 + vwnd[offset] ** 2) ** 0.5 + if val.mask.all(): + LOG.info("No wind for offset: %s", offset) + continue + if runningsum is None: + runningsum = val + else: + runningsum += val + hours += 1 + if hours > 0: + ds["wind_speed"].values = runningsum / hours + for vname in ( + "high_tmpk low_tmpk p01d high_soil4t avg_dwpk " + "low_soil4t high_tmpk_12z low_tmpk_12z p01d_12z" + ).split(): + res = None + aggfunc = np.ma.max + if vname.startswith("p01d"): + aggfunc = np.ma.sum # was np.nansum, better check this + if vname == "avg_dwpk": + aggfunc = np.ma.mean + ncvarname = ( + vname.replace("high_", "") + .replace("low_", "") + .replace("_12z", "") + .replace("avg_", "") + ) + ncvarname = "p01m" if ncvarname == "p01d" else ncvarname + for pair in pairs12z if vname.endswith("12z") else pairs: + ncfn = iemre.get_hourly_ncname(pair[0].year) + tidx1 = iemre.hourly_offset(pair[0]) + tidx2 = iemre.hourly_offset(pair[1]) + tslice = slice(tidx1, tidx2 + 1) + if vname.startswith("p01d"): + # Caution: precip is stored arrears, so we have ugliness + tslice = slice(0 if tidx1 == 0 else tidx1 + 1, tidx2 + 2) + LOG.info("Using %s[%s] for %s", ncfn, tslice, vname) + with ncopen(ncfn, timeout=600) as nc: + ncvar = nc.variables[ncvarname] + ncval = aggfunc(ncvar[tslice], axis=0) + if res is None: + res = ncval + else: + res = aggfunc([res, ncval], axis=0) + # NB: these are crude bias offsets involved when using hourly data + if vname in ["high_tmpk", "high_tmpk_12z"]: + res += 0.8 + if vname in ["low_tmpk", "low_tmpk_12z"]: + res -= 0.8 + ds[vname].values = res def use_climodat_12z(ts, ds): @@ -246,81 +242,18 @@ def use_climodat_12z(ts, ds): if len(df.index) < 50: LOG.warning("Failed quorum") return - res = generic_gridder(df, "highdata") - ds["high_tmpk_12z"].values = convert_value(res, "degF", "degK") - - res = generic_gridder(df, "lowdata") - ds["low_tmpk_12z"].values = convert_value(res, "degF", "degK") - - res = generic_gridder(df, "snowdata") - ds["snow_12z"].values = convert_value(res, "inch", "millimeter") - - res = generic_gridder(df, "snowddata") - ds["snowd_12z"].values = convert_value(res, "inch", "millimeter") - - -def use_coop_12z(ts, ds): - """Use the COOP data for gridding""" - LOG.info("12z hi/lo for %s", ts) - mybuf = 2.0 - giswkt = "SRID=4326;POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))" % ( - iemre.WEST - mybuf, - iemre.SOUTH - mybuf, - iemre.WEST - mybuf, - iemre.NORTH + mybuf, - iemre.EAST + mybuf, - iemre.NORTH + mybuf, - iemre.EAST + mybuf, - iemre.SOUTH - mybuf, - iemre.WEST - mybuf, - iemre.SOUTH - mybuf, - ) - with get_sqlalchemy_conn("iem") as conn: - df = pd.read_sql( - f""" - SELECT ST_x(s.geom) as lon, ST_y(s.geom) as lat, s.state, - s.id as station, s.name as name, - (CASE WHEN pday >= 0 then pday else null end) as precipdata, - (CASE WHEN snow >= 0 then snow else null end) as snowdata, - (CASE WHEN snowd >= 0 then snowd else null end) as snowddata, - (CASE WHEN max_tmpf > -50 and max_tmpf < 130 - then max_tmpf else null end) as highdata, - (CASE WHEN min_tmpf > -50 and min_tmpf < 95 - then min_tmpf else null end) as lowdata - from summary_{ts.year} c, stations s WHERE day = %s and - ST_Contains(ST_GeomFromEWKT(%s), geom) and s.network ~* 'COOP' - and c.iemid = s.iemid and - extract(hour from c.coop_valid at time zone s.tzname) between 4 and 11 - """, - conn, - params=( - ts, - giswkt, - ), - ) - LOG.info("loaded %s rows from iemaccess database", len(df.index)) - - # Require that high > low before any gridding, accounts for some COOP - # sites that only report TOB and not 24 hour high/low - df.loc[df["highdata"] <= df["lowdata"], ["highdata", "lowdata"]] = None - - if len(df.index) < 4: - LOG.warning("Failed quorum") - res = generic_gridder(df, "highdata") - if res is not None: + if ts.year < 1951: + res = generic_gridder(df, "highdata") ds["high_tmpk_12z"].values = convert_value(res, "degF", "degK") - res = generic_gridder(df, "lowdata") - if res is not None: + res = generic_gridder(df, "lowdata") ds["low_tmpk_12z"].values = convert_value(res, "degF", "degK") res = generic_gridder(df, "snowdata") - if res is not None: - ds["snow_12z"].values = convert_value(res, "inch", "millimeter") + ds["snow_12z"].values = convert_value(res, "inch", "millimeter") res = generic_gridder(df, "snowddata") - if res is not None: - ds["snowd_12z"].values = convert_value(res, "inch", "millimeter") + ds["snowd_12z"].values = convert_value(res, "inch", "millimeter") def use_asos_daily(ts, ds): @@ -455,37 +388,24 @@ def use_climodat_daily(ts, ds): def workflow(ts): """Do Work""" - today = datetime.date.today() # load up our current data ds = iemre.get_grids(ts) + LOG.info("loaded %s variables from IEMRE database", len(ds)) + # rsds -> grid_rsds.py # power_swdn -> TODO - LOG.info("loaded %s variables from IEMRE database", len(ds)) - if ts.year > 1927: + # high_tmpk, low_tmpk, p01d, wind_speed, min_rh, max_rh, high_soil4t, + # low_soil4t, high_tmpk_12z low_tmpk_12z p01d_12z + if ts.year > 1949: + copy_iemre_hourly(ts, ds) + else: LOG.info("Using ASOS for daily summary variables") - # high_tmpk, low_tmpk, p01d - # avg_dwpk, wind_speed, min_rh, max_rh use_asos_daily(ts, ds) - if ts < today: - # high_tmpk, low_tmpk, p01d use_climodat_daily(ts, ds) - if ts.year > 1996: - # high_soil4t, low_soil4t, wind_speed, p01d - copy_iemre_hourly(ts, ds) - if ts.year > 2011: - # high_tmpk_12z low_tmpk_12z p01d_12z snow_12z snowd_12z - use_coop_12z(ts, ds) - else: - # high_tmpk_12z low_tmpk_12z p01d_12z snow_12z snowd_12z - use_climodat_12z(ts, ds) - if ts.year > 1996: - # p01d_12z - copy_iemre_12z(ts, ds) - if ts.year < 1951: - LOG.info("Verbatim copy of daily hi,lo,pcpn to 12z") - for col in ["high_tmpk", "low_tmpk", "p01d"]: - ds[f"{col}_12z"].values = ds[col] + + # snow_12z snowd_12z + use_climodat_12z(ts, ds) for vname in list(ds.keys()): # Some grids are too tight to the CONUS boundary, so we smear things @@ -505,13 +425,14 @@ def workflow(ts): ) -def main(argv): +@click.command() +@click.option("--date", "valid", type=click.DateTime(), help="Date") +def main(valid: datetime.datetime): """Go Main Go""" - ts = datetime.date(int(argv[1]), int(argv[2]), int(argv[3])) - LOG.info("Run %s", ts) - workflow(ts) + LOG.info("Run %s", valid) + workflow(valid.date()) LOG.info("Done.") if __name__ == "__main__": - main(sys.argv) + main() diff --git a/scripts/iemre/hourly_analysis.py b/scripts/iemre/hourly_analysis.py index d921b5becf..bbb9876ce4 100644 --- a/scripts/iemre/hourly_analysis.py +++ b/scripts/iemre/hourly_analysis.py @@ -9,6 +9,7 @@ import click import numpy as np import pandas as pd +import pint import pygrib from affine import Affine from metpy.calc import wind_components @@ -17,7 +18,7 @@ from pyiem import iemre from pyiem.database import get_sqlalchemy_conn from pyiem.iemre import grb2iemre, hourly_offset, reproject2iemre -from pyiem.util import logger, ncopen +from pyiem.util import archive_fetch, logger, ncopen # Prevent invalid value encountered in cast warnings.simplefilter("ignore", RuntimeWarning) @@ -79,61 +80,57 @@ def use_hrrr_soilt(ts): """Verbatim copy HRRR, if it exists.""" # We may be running close to real-time, so it makes some sense to take # files from the recent past. - grbfn = None for offset in range(5): fn = (ts - datetime.timedelta(hours=offset)).strftime( - "/mesonet/ARCHIVE/data/%Y/%m/%d/model/hrrr/%H/" - "hrrr.t%Hz.3kmf00.grib2" + "%Y/%m/%d/model/hrrr/%H/hrrr.t%Hz.3kmf00.grib2" ) - if os.path.isfile(fn): - grbfn = fn - break - if grbfn is None: - LOG.info("Failed to find any HRRR grib files") - return None - try: - grbs = pygrib.open(fn) - for grb in grbs: - if grb.shortName != "st" or str(grb).find("level 0.1 m") == -1: + with archive_fetch(fn) as fn: + if fn is None: continue - return grb2iemre(grb) - grbs.close() - except Exception as exp: - LOG.info("%s exp:%s", fn, exp) + try: + grbs = pygrib.open(fn) + for grb in grbs: + if ( + grb.shortName != "st" + or str(grb).find("level 0.1 m") == -1 + ): + continue + return grb2iemre(grb) + grbs.close() + except Exception as exp: + LOG.info("%s exp:%s", fn, exp) return None def use_rtma(ts, kind): """Verbatim copy RTMA, if it exists.""" - fn = ts.strftime( - "/mesonet/ARCHIVE/data/%Y/%m/%d/model/rtma/%H/" - "rtma.t%Hz.awp2p5f000.grib2" - ) - tasks = { - "wind": [ - "10u", - "10v", - ], - "tmp": [ - "2t", - ], - "dwp": [ - "2d", - ], - } - if not os.path.isfile(fn): - LOG.info("Failed to find %s", fn) + if ts.year < 2011: return None - res = [] - try: - grbs = pygrib.open(fn) - for task in tasks[kind]: - grb = grbs.select(shortName=task)[0] - res.append(grb2iemre(grb)) - return res - except Exception as exp: - LOG.info("%s exp:%s", fn, exp) - return None + ppath = ts.strftime("%Y/%m/%d/model/rtma/%H/rtma.t%Hz.awp2p5f000.grib2") + with archive_fetch(ppath) as fn: + if fn is None: + return None + tasks = { + "wind": [ + "10u", + "10v", + ], + "tmp": [ + "2t", + ], + "dwp": [ + "2d", + ], + } + res = [] + try: + grbs = pygrib.open(fn) + for task in tasks[kind]: + grb = grbs.select(shortName=task)[0] + res.append(grb2iemre(grb)) + return res[0] if len(res) == 1 else res + except Exception as exp: + LOG.info("%s exp:%s", fn, exp) def grid_wind(df, domain): @@ -306,7 +303,6 @@ def grid_hour(ts): dwpf = np.where(mask, tmpf, dwpf) write_grid(ts, "tmpk", masked_array(tmpf, data_units="degF").to("degK")) write_grid(ts, "dwpk", masked_array(dwpf, data_units="degF").to("degK")) - res = grid_skyc(df, domain) if res is None: LOG.warning("No obs for skyc at %s", ts) @@ -320,6 +316,8 @@ def write_grid(valid, vname, grid): This is isolated so that we don't 'lock' up our file while intensive work is done """ + if isinstance(grid, pint.Quantity): + grid = grid.m offset = iemre.hourly_offset(valid) with ncopen(iemre.get_hourly_ncname(valid.year), "a", timeout=300) as nc: LOG.info( From adacbf22fb7a163c5e3c0eca29fd49fa158ee868 Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 9 Feb 2024 12:07:08 -0600 Subject: [PATCH 3/4] address some lint --- scripts/iemre/hourly_analysis.py | 3 +++ scripts/iemre/init_dailyc.py | 3 ++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/scripts/iemre/hourly_analysis.py b/scripts/iemre/hourly_analysis.py index bbb9876ce4..270c388211 100644 --- a/scripts/iemre/hourly_analysis.py +++ b/scripts/iemre/hourly_analysis.py @@ -131,6 +131,9 @@ def use_rtma(ts, kind): return res[0] if len(res) == 1 else res except Exception as exp: LOG.info("%s exp:%s", fn, exp) + finally: + grbs.close() + return None def grid_wind(df, domain): diff --git a/scripts/iemre/init_dailyc.py b/scripts/iemre/init_dailyc.py index 6b2e45d58b..94c9152b4a 100644 --- a/scripts/iemre/init_dailyc.py +++ b/scripts/iemre/init_dailyc.py @@ -5,8 +5,9 @@ import geopandas as gpd import numpy as np from pyiem import iemre +from pyiem.database import get_dbconn from pyiem.grid.zs import CachingZonalStats -from pyiem.util import get_dbconn, logger, ncopen +from pyiem.util import logger, ncopen LOG = logger() From 7ca36f97f4b56b2a3806613b52d9d7e7c4dddef9 Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 9 Feb 2024 12:10:45 -0600 Subject: [PATCH 4/4] mnt: address lint --- scripts/dbutil/add_iem_data_entry.py | 3 +- scripts/dbutil/add_summary_entries.py | 3 +- scripts/dbutil/clean_afos.py | 3 +- scripts/dbutil/clean_mos.py | 3 +- scripts/dbutil/clean_telemetry.py | 3 +- scripts/dbutil/clean_unknown_hads.py | 3 +- scripts/iemre/daily_analysis.py | 44 +++++++++++++-------------- 7 files changed, 34 insertions(+), 28 deletions(-) diff --git a/scripts/dbutil/add_iem_data_entry.py b/scripts/dbutil/add_iem_data_entry.py index 9e832db260..2415dfb857 100644 --- a/scripts/dbutil/add_iem_data_entry.py +++ b/scripts/dbutil/add_iem_data_entry.py @@ -4,7 +4,8 @@ """ import datetime -from pyiem.util import get_dbconnc, logger +from pyiem.database import get_dbconnc +from pyiem.util import logger LOG = logger() diff --git a/scripts/dbutil/add_summary_entries.py b/scripts/dbutil/add_summary_entries.py index 2310d83b7d..6a26bbc42e 100644 --- a/scripts/dbutil/add_summary_entries.py +++ b/scripts/dbutil/add_summary_entries.py @@ -2,7 +2,8 @@ import datetime import sys -from pyiem.util import get_dbconn, logger +from pyiem.database import get_dbconn +from pyiem.util import logger LOG = logger() diff --git a/scripts/dbutil/clean_afos.py b/scripts/dbutil/clean_afos.py index 84c49afc1d..94841d03fb 100644 --- a/scripts/dbutil/clean_afos.py +++ b/scripts/dbutil/clean_afos.py @@ -3,8 +3,9 @@ called from RUN_2AM.sh """ +from pyiem.database import get_dbconn from pyiem.reference import state_names -from pyiem.util import get_dbconn, logger +from pyiem.util import logger LOG = logger() diff --git a/scripts/dbutil/clean_mos.py b/scripts/dbutil/clean_mos.py index b6f739e176..22c9ef2a49 100644 --- a/scripts/dbutil/clean_mos.py +++ b/scripts/dbutil/clean_mos.py @@ -3,7 +3,8 @@ Run from `RUN_2AM.sh` script. """ -from pyiem.util import get_dbconn, logger +from pyiem.database import get_dbconn +from pyiem.util import logger LOG = logger() diff --git a/scripts/dbutil/clean_telemetry.py b/scripts/dbutil/clean_telemetry.py index 2e5bfd40d8..a12de25f34 100644 --- a/scripts/dbutil/clean_telemetry.py +++ b/scripts/dbutil/clean_telemetry.py @@ -4,7 +4,8 @@ """ from datetime import datetime, timedelta -from pyiem.util import get_dbconnc, logger +from pyiem.database import get_dbconnc +from pyiem.util import logger LOG = logger() diff --git a/scripts/dbutil/clean_unknown_hads.py b/scripts/dbutil/clean_unknown_hads.py index e058e96a9d..66ec4116ab 100644 --- a/scripts/dbutil/clean_unknown_hads.py +++ b/scripts/dbutil/clean_unknown_hads.py @@ -5,7 +5,8 @@ Run from RUN_2AM.sh """ -from pyiem.util import get_dbconn, logger +from pyiem.database import get_dbconn +from pyiem.util import logger LOG = logger() diff --git a/scripts/iemre/daily_analysis.py b/scripts/iemre/daily_analysis.py index a50a6aaae1..aaf67561cb 100644 --- a/scripts/iemre/daily_analysis.py +++ b/scripts/iemre/daily_analysis.py @@ -133,28 +133,28 @@ def copy_iemre_hourly(ts, ds): ds[vname].values = res # One off - for vname in ["wind_speed"]: - hours = 0 - runningsum = None - for pair in pairs: - ncfn = iemre.get_hourly_ncname(pair[0].year) - tidx1 = iemre.hourly_offset(pair[0]) - tidx2 = iemre.hourly_offset(pair[1]) - with ncopen(ncfn, timeout=600) as nc: - uwnd = nc.variables["uwnd"] - vwnd = nc.variables["vwnd"] - for offset in range(tidx1, tidx2 + 1): - val = (uwnd[offset] ** 2 + vwnd[offset] ** 2) ** 0.5 - if val.mask.all(): - LOG.info("No wind for offset: %s", offset) - continue - if runningsum is None: - runningsum = val - else: - runningsum += val - hours += 1 - if hours > 0: - ds["wind_speed"].values = runningsum / hours + hours = 0 + runningsum = None + for pair in pairs: + ncfn = iemre.get_hourly_ncname(pair[0].year) + tidx1 = iemre.hourly_offset(pair[0]) + tidx2 = iemre.hourly_offset(pair[1]) + with ncopen(ncfn, timeout=600) as nc: + uwnd = nc.variables["uwnd"] + vwnd = nc.variables["vwnd"] + for offset in range(tidx1, tidx2 + 1): + val = (uwnd[offset] ** 2 + vwnd[offset] ** 2) ** 0.5 + if val.mask.all(): + LOG.info("No wind for offset: %s", offset) + continue + if runningsum is None: + runningsum = val + else: + runningsum += val + hours += 1 + if hours > 0: + ds["wind_speed"].values = runningsum / hours + # ----------------------------------------------------------------- for vname in ( "high_tmpk low_tmpk p01d high_soil4t avg_dwpk " "low_soil4t high_tmpk_12z low_tmpk_12z p01d_12z"