WIP wind

Smithsonian · Feb 21, 2022 · bba5fa6 · bba5fa6
1 parent 73d9ab0
commit bba5fa6
Show file tree

Hide file tree

Showing 5 changed files with 189 additions and 17 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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&region=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
@@ -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