diff --git a/eht_met_forecast/am.py b/eht_met_forecast/am.py index 03f29f4..ca6330b 100644 --- a/eht_met_forecast/am.py +++ b/eht_met_forecast/am.py @@ -66,10 +66,75 @@ def grid_interp_vector(a, b, u, v): RH_TOP_PLEVEL = 29. STRAT_H2O_VMR = 5e-6 +# https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2.0p25.anl.shtml +# https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2.0p25.f000.shtml +# https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2.0p25.f003.shtml -- layers beyond 0 + +# https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2b.0p25.anl.shtml +# https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2b.0p25.f000.shtml +# https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2b.0p25.f003.shtml -- layers beyond 0 + +scalar_gribs = [ # this table is not yet used + {'lev': 'surface', 'var': ['CSNOW'], 'name': 'Categorial snow', 'level': [0]}, # f000 and f003 but not anl (?) and it's "3 hour fcst" and "0-3 hour ave" in "Forecast Valid"... "analysis" in f000 + {'lev': 'surface', 'var': ['CICEP'], 'name': 'Categorical ice pellets', 'level': [0]}, + {'lev': 'surface', 'var': ['CFRZR'], 'name': 'Categorical freezing rain', 'level': [0]}, + {'lev': 'surface', 'var': ['CRAIN'], 'name': 'Categorical rain', 'level': [0]}, + {'lev': 'surface', 'var': ['GUST'], 'name': 'Wind speed (gust)', 'level': [0]}, # f000 and f003 surface and not anl? "3 hour fcst" at f003, "analysis" in f000 +] + +vector_gribs = [ # this table is not yet used + # these come in pairs and need a different interpolation function + {'lev': ['max_wind'], 'var': ['GUST'], 'name': 'U component of wind', 'level': [0], 'ourname': 'wind gust'}, # f0003 "3 hour fcst" vs "analysis" and "analyis" for anl and f000 + {'lev': ['max_wind'], 'var': ['GUST'], 'name': 'V component of wind', 'level': [0], 'ourname': 'wind gust'}, + {'lev': ['1'], 'var': ['UGRD'], 'name': 'U component of wind', 'level': [1], 'ourname': 'low wind'}, # anl: level 1 and up, also PV=blah, f000: "planwetary boundary layer", layer 1, "10 m above ground" + {'lev': ['1'], 'var': ['VGRD'], 'name': 'V component of wind', 'level': [1], 'ourname': 'low wind'}, +] + + +def grib2_to_extra_information(grbindx, u, v): + ''' +pygrib clues +grbindx.select only works on the named args named in the pygrib.index() call (name and level, here) +didn't have much luck having more than 2 index names in the pygrib.index call -- could probably make multiple indices +or you can use pygrib.open() instead, to not have an index, it's supposedly slower but these grib files are small (2k/station) + ''' + ret = {} + + try: + k = 'csnow' + ret['csnow'] = (grid_interp(grbindx.select(name='Categorical snow', level=0)[0].values, u, v)) + k = 'cicep' + ret['cicep'] = (grid_interp(grbindx.select(name='Categorical ice pellets', level=0)[0].values, u, v)) + k = 'cfrzr' + ret['cfrzr'] = (grid_interp(grbindx.select(name='Categorical freezing rain', level=0)[0].values, u, v)) + k = 'crain' + ret['crain'] = (grid_interp(grbindx.select(name='Categorical rain', level=0)[0].values, u, v)) + + # appears for level_surface + k = 'wgust' + ret['wgust'] = (grid_interp(grbindx.select(name='Wind speed (gust)', level=0)[0].values, u, v)) + + # this one is mediated by lev_max_wind + k = 'max wind u' + a = grbindx.select(name='U component of wind', level=0)[0].values + k = 'max wind v' + b = grbindx.select(name='V component of wind', level=0)[0].values + ret['max_wind'] = grid_interp_vector(a, b, u, v) + + # apparently level 1 is as low as you can go + k = 'level 1 wind u' + a = grbindx.select(name='U component of wind', level=1)[0].values + k = 'level 1 wind v' + b = grbindx.select(name='V component of wind', level=1)[0].values + ret['surface_wind'] = grid_interp_vector(a, b, u, v) + except Exception as e: + print('key:', k, 'exception:', e, file=sys.stderr) + raise + return ret -def grib2_to_am_layers(gribname, lat, lon, alt): - grbindx = pygrib.index(gribname, "name", "level") +def grib2_to_am_layers(gribname, lat, lon, alt): + grbindx = pygrib.index(gribname, "name", "level") # normal code # in memory -- not sure what syntax actually works for this? # need to .index() after creation # gribfile = pygrib.fromstring(grib_buffer) @@ -87,6 +152,8 @@ def grib2_to_am_layers(gribname, lat, lon, alt): cloud_lmr = [] cloud_imr = [] + extra = grib2_to_extra_information(grbindx, u, v) + for i, lev in enumerate(LEVELS): Pbase.append(lev) try: @@ -137,7 +204,19 @@ def grib2_to_am_layers(gribname, lat, lon, alt): # this is not unusual cloud_imr.append(0.0) - return Pbase, z, T, o3_vmr, RH, cloud_lmr, cloud_imr + return Pbase, z, T, o3_vmr, RH, cloud_lmr, cloud_imr, extra + + +def print_extra(gfs_cycle, forecast_hour, extra): + # extra is a dict + # gfs_cycle forms the filename + # forecast hour + fname = gfs_cycle.strftime(GFS_TIMESTAMP) + dt_forecast_hour = gfs_cycle + datetime.timedelta(hours=forecast_hour) + rowname = dt_forecast_hour.strftime(GFS_TIMESTAMP) + # XXX write a csv line + # the normal output file is f, passed to print_final_output, passed into compute_one_hour, by make_forecast_table, which is called by cli.main + pass def print_am_header(gfs_cycle, forecast_hour, lat, lon, alt): diff --git a/eht_met_forecast/cli.py b/eht_met_forecast/cli.py index 0ebd694..e81eb8f 100644 --- a/eht_met_forecast/cli.py +++ b/eht_met_forecast/cli.py @@ -4,10 +4,11 @@ import os from argparse import ArgumentParser from collections import defaultdict +import csv from .constants import GFS_TIMESTAMP, GFS_TIMESTAMP_FULL from .timer_utils import dump_latency_histograms -from .gfs import latest_gfs_cycle_time +from .gfs import latest_gfs_cycle_time, jiggle from .core import read_stations, ok, make_forecast_table, dump_stats @@ -92,6 +93,7 @@ def main(args=None): stats['stations'] = [] stats['gfs_time'] = cycles[0].strftime(GFS_TIMESTAMP) stats['start'] = datetime.datetime.now(datetime.timezone.utc).strftime(GFS_TIMESTAMP_FULL) + time.sleep(jiggle(15) - 15) # 0-5 seconds t0 = time.time() for vex in stations: @@ -116,11 +118,20 @@ def main(args=None): os.makedirs(outdir, exist_ok=True) if args.stdout: f = sys.stdout + f2 = None else: f = open(outfile, 'w') + f2 = open(outfile+'.extra', 'w', newline='') + + if f2: + fieldnames = [ + 'date', 'csnow', 'cicep', 'cfrzr', 'crain', 'wgust', 'max_wind', 'surface_wind', + ] + f2 = csv.DictWriter(f2, fieldnames=fieldnames, delimiter=' ') + f2.writeheader() try: - make_forecast_table(station, gfs_cycle, f, wait=args.wait, verbose=args.verbose, one=args.one, stats=stats) + make_forecast_table(station, gfs_cycle, f, f2, wait=args.wait, verbose=args.verbose, one=args.one, stats=stats) except TimeoutError: # raised by gfs.py print('Gave up on {} {}'.format(station, gfs_cycle), file=sys.stderr) diff --git a/eht_met_forecast/core.py b/eht_met_forecast/core.py index 42ec466..3c0c651 100644 --- a/eht_met_forecast/core.py +++ b/eht_met_forecast/core.py @@ -44,15 +44,18 @@ def gfs15_to_am10(lat, lon, alt, gfs_cycle, forecast_hour, wait=False, verbose=F grib_buffer = download_gfs(lat, lon, alt, gfs_cycle, forecast_hour, wait=wait, verbose=verbose, stats=stats) grib_problem = False + # development hint: use delete=False to save all of these with tempfile.NamedTemporaryFile(mode='wb', prefix='temp-', suffix='.grb') as f: f.write(grib_buffer) f.flush() + try: - Pbase, z, T, o3_vmr, RH, cloud_lmr, cloud_imr = grib2_to_am_layers(f.name, lat, lon, alt) + Pbase, z, T, o3_vmr, RH, cloud_lmr, cloud_imr, extra = grib2_to_am_layers(f.name, lat, lon, alt) except Exception as e: # example: RuntimeError: b'End of resource reached when reading message' # example: UserWarning: file temp.grb has multi-field messages, keys inside multi-field messages will not be indexed correctly grib_problem = str(e) + print('problem reading grib:', grib_problem, file=sys.stderr) if not grib_problem: my_stdout = io.StringIO() @@ -64,6 +67,7 @@ def gfs15_to_am10(lat, lon, alt, gfs_cycle, forecast_hour, wait=False, verbose=F # example: ZeroDivisionError after a bunch of # ECCODES INFO : grib_file_open: cannot open file foo.grb (No such file or directory) grib_problem = str(e) + print('problem printing am', grib_problem, file=sys.stderr) if grib_problem: with tempfile.NamedTemporaryFile(mode='wb', prefix='layers-err-', suffix='.grb', dir='.', delete=False) as tfile: @@ -74,9 +78,9 @@ def gfs15_to_am10(lat, lon, alt, gfs_cycle, forecast_hour, wait=False, verbose=F with open(fname, 'w') as f: print(grib_problem, file=f) print('lat', lat, 'lon', lon, 'alt', alt, 'gfs_cycle', gfs_cycle, 'forecast_hour', forecast_hour, file=f) - return + return None, None - return my_stdout.getvalue() + return my_stdout.getvalue(), extra def print_final_output(gfs_timestamp, tau, Tb, pwv, lwp, iwp, o3, f, verbose=False): @@ -87,11 +91,17 @@ def print_final_output(gfs_timestamp, tau, Tb, pwv, lwp, iwp, o3, f, verbose=Fal sys.stderr.flush() -def compute_one_hour(site, gfs_cycle, forecast_hour, f, wait=False, verbose=False, stats=None): +def print_extra(fcast_pretty, extra, f2, verbose=False): + if f2: + extra['date'] = fcast_pretty + f2.writerow(extra) + + +def compute_one_hour(site, gfs_cycle, forecast_hour, f, f2, wait=False, verbose=False, stats=None): if verbose: print(site['name'], 'fetching for hour', forecast_hour, file=sys.stderr) with record_latency('fetch gfs data'): - layers_amc = gfs15_to_am10(site['lat'], site['lon'], site['alt'], gfs_cycle, forecast_hour, wait=wait, verbose=verbose, stats=stats) + layers_amc, extra = gfs15_to_am10(site['lat'], site['lon'], site['alt'], gfs_cycle, forecast_hour, wait=wait, verbose=verbose, stats=stats) if layers_amc is None: return # no line emitted @@ -123,18 +133,20 @@ def compute_one_hour(site, gfs_cycle, forecast_hour, f, wait=False, verbose=Fals tfile.write(am_output) return # no line emitted - print_final_output(dt_forecast_hour.strftime(GFS_TIMESTAMP), tau, Tb, pwv, lwp, iwp, o3, f, verbose=verbose) + fcast_pretty = dt_forecast_hour.strftime(GFS_TIMESTAMP) + print_final_output(fcast_pretty, tau, Tb, pwv, lwp, iwp, o3, f, verbose=verbose) + print_extra(fcast_pretty, extra, f2, verbose=verbose) time.sleep(1) -def make_forecast_table(site, gfs_cycle, f, wait=False, verbose=False, one=False, stats=None): +def make_forecast_table(site, gfs_cycle, f, f2, wait=False, verbose=False, one=False, stats=None): print_table_line(table_header, f) for forecast_hour in range(0, 121): - compute_one_hour(site, gfs_cycle, forecast_hour, f, wait=wait, verbose=verbose, stats=stats) + compute_one_hour(site, gfs_cycle, forecast_hour, f, f2, wait=wait, verbose=verbose, stats=stats) if one: return for forecast_hour in range(123, 385, 3): - compute_one_hour(site, gfs_cycle, forecast_hour, f, wait=wait, verbose=verbose, stats=stats) + compute_one_hour(site, gfs_cycle, forecast_hour, f, f2, wait=wait, verbose=verbose, stats=stats) def read_stations(filename): diff --git a/eht_met_forecast/gfs.py b/eht_met_forecast/gfs.py index 00c0644..7c309af 100644 --- a/eht_met_forecast/gfs.py +++ b/eht_met_forecast/gfs.py @@ -43,7 +43,74 @@ def form_gfs_download_url(lat, lon, alt, gfs_cycle, forecast_hour): for lev in LEVELS: params['lev_{:d}_mb'.format(lev)] = 'on' - VARIABLES = ("CLWMR", "ICMR", "HGT", "O3MR", "RH", "TMP") + params['lev_surface'] = 'on' # CRAIN, etc. maps to level=0 + params['lev_max_wind'] = 'on' # UGRD, VGRD, maps to level=0 + + # lev_ that cause 500 errors: 'lev_0_mb' 'lev_0' 'lev_10_m' 'lev_20_m' 'lev_maxWind' + + VARIABLES = ["CLWMR", "ICMR", "HGT", "O3MR", "RH", "TMP"] # for AM + #VARIABLES += ["WIND"] # wind -- not available at all in GFS + VARIABLES += ["UGRD", "VGRD"] # wind -- meters/sec -- works in levels but really we want level 0? no joy for lev_surface either + VARIABLES += ["CRAIN", "CFRZR", "CICEP", "CSNOW"] # yes/no 1/0 rain, freezing rain, ice pellets, snow -- works for lev_surface + #VARIABLES += ["POP", "CPOFP", "CPOZP"] # probability of precip -- not in the GFS docs, so no joy + VARIABLES += ["GUST"] # works with lev_surface and then level=0 + ''' + these are all of the levels -- get_gfs.pl download of all variables + + 95:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 96:Temperature:K (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 97:Relative humidity:% (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 98:Specific humidity:kg kg**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 99:Vertical velocity:Pa s**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 100:Geometric vertical velocity:m s**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 101:U component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 102:V component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 103:Absolute vorticity:s**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 104:Ozone mixing ratio:kg kg**-1 (instant):regular_ll:isobaricInhPa:level 100 Pa:fcst time 0 hrs:from 202202060000 + 642:Temperature:K (instant):regular_ll:heightAboveGround:level 100 m:fcst time 0 hrs:from 202202060000 + 643:100 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 100 m:fcst time 0 hrs:from 202202060000 + 644:100 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 100 m:fcst time 0 hrs:from 202202060000 + + level 0, yet to figure out how to fetch them. maxWind? maybe this is at any altitude?! + the level 100+ winds above work with the normal levels + + 626:U component of wind:m s**-1 (instant):regular_ll:maxWind:level 0:fcst time 0 hrs:from 202202060000 + 627:V component of wind:m s**-1 (instant):regular_ll:maxWind:level 0:fcst time 0 hrs:from 202202060000 + https://www.weatheronline.co.uk/cgi-bin/expertcharts?LANG=en&CONT=ukuk&MODELL=gfs&MODELLTYP=1&VAR=uv10&INFO=1& + "surface wind" which is wind at 10 meters above the ground + https://www.tropicaltidbits.com/analysis/models/?model=gfs®ion=us&pkg=mslp_wind&runtime=2022020618&fh=6 + also says "10m wind" + + 585:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202202060000 + 586:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202202060000 + + + + these next ones are the lowest height + + level 0, also called surface, actually works with lev_surface + + 590:Categorical snow:(Code table 4.222) (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + 591:Categorical ice pellets:(Code table 4.222) (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + 592:Categorical freezing rain:(Code table 4.222) (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + 593:Categorical rain:(Code table 4.222) (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + + these are the only 'precip' or 'percent' entries: + 588:Percent frozen precipitation:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + PRATE 589:Precipitation rate:kg m**-2 s**-1 (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + PWAT 604:Precipitable water:kg m**-2 (instant):regular_ll:atmosphereSingleLayer:level 0 considered as a single layer:fcst time 0 hrs:from 202202060000 + CWAT + + test me: not in the inventory (!) + not useful, it's at any altitue https://www.weatheronline.co.uk/cgi-bin/expertcharts?MODELL=gfs&MODELLTYP=1&VAR=boen&INFO=1 + 14:Wind speed (gust):m s**-1 (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202202060000 + would be useful to have this, because it might be related to the atmospheric coherance time + 550 and 551? https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2.0p25.anl.shtml + + + + ''' + for var in VARIABLES: params['var_' + var] = 'on' @@ -154,7 +221,8 @@ def fetch_gfs_download(url, params, wait=False, verbose=False, stats=None): print(" Retrying...", file=sys.stderr) if r: try: # I don't think this can fail, but anyway - print(' Content was:', r.content[:100]) + if len(r.content): + print(' Content was:', r.content[:100]) except Exception: pass time.sleep(retry_duration) diff --git a/tests/test_cli.py b/tests/test_cli.py index c8feb09..715d621 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -37,8 +37,10 @@ def test_cli(capsys): # AM 11.0 gfs15_to_am #stdout = '00:00 7.6246e-02 2.3999e+01 1.4787e+00 0.0000e+00 0.0000e+00 2.7655e+02\n' # AM 11.0 gfs16_to_am and bugfixes 2/5/22 - stdout = '00:00 7.6247e-02 2.3999e+01 1.4794e+00 0.0000e+00 0.0000e+00 2.7655e+02\n' + #stdout = '00:00 7.6247e-02 2.3999e+01 1.4794e+00 0.0000e+00 0.0000e+00 2.7655e+02\n' # AM 12.0 gfs16 is the same + # new test.grb with wind + stdout = '00:00 1.1151e-01 3.2311e+01 2.4032e+00 0.0000e+00 0.0000e+00 2.6272e+02\n' stdout = '20200316+18:' + stdout