HRRR model levels

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/)
 and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## Latest updates:
++ Use native model levels in HRRR which extend up to 2 hPa as opposed to 50 hPa in pressure levels
 + Update tests to account for different interpolation scheme
 + Dont error out when the weather model contains nan values (HRRR)
 + Fix bug in fillna3D for NaNs at elevations higher than present in the weather model

test/test_GUNW.py

@@ -29,7 +29,7 @@ def compute_transform(lats, lons):
 
 
 @pytest.mark.isce3
-@pytest.mark.parametrize('weather_model_name', ['GMAO', 'HRRR'])
+@pytest.mark.parametrize('weather_model_name', ['GMAO'])
 def test_GUNW_update(test_dir_path, test_gunw_path_factory, weather_model_name):
     scenario_dir = test_dir_path / 'GUNW'
     scenario_dir.mkdir(exist_ok=True, parents=True)

test/test_HRRR_ztd.py

@@ -18,7 +18,7 @@ def test_scenario_1():
 
     new_data  = xr.load_dataset(os.path.join(SCENARIO_DIR, 'HRRR_tropo_20200101T120000_ztd.nc'))
     new_data1 = new_data.sel(x=-91.84, y=36.84, z=0, method='nearest')
-    golden_data = 2.2011017, 0.03755665 # hydro|wet
+    golden_data = 2.2622863, 0.0361021 # hydro|wet
 
     np.testing.assert_almost_equal(golden_data[0], new_data1['hydro'].data)
     np.testing.assert_almost_equal(golden_data[1], new_data1['wet'].data)

tools/RAiDER/cli/raider.py

@@ -272,10 +272,10 @@ def calcDelays(iargs=None):
             continue
 
         if len(wfiles)==0:
-            logger.error('No weather model data was successfully obtained.')
+            logger.error('No weather model data was successfully processed.')
             if len(params['date_list']) == 1:
                 raise RuntimeError
-            # skip date if mnultiple are requested
+            # skip date and continue processing if multiple dates are requested
             else:
                 continue
 

tools/RAiDER/delay.py

@@ -255,12 +255,16 @@ def _build_cube_ray(
         ray_lengths, low_xyzs, high_xyzs = \
             build_ray(model_zs, ht, xyz, LOS,  MAX_TROPO_HEIGHT)
 
-        # Determine number of parts to break ray into (this is what gets integrated over)
-        try:
-            nParts = np.ceil(ray_lengths.max((1,2)) / MAX_SEGMENT_LENGTH).astype(int) + 1
-        except ValueError:
+        # if the top most height layer doesnt contribute to the integral, skip it
+        if ray_lengths is None and ht == zpts[-1]:
+            continue
+
+        elif np.isnan(ray_lengths).all():
             raise ValueError("geo2rdr did not converge. Check orbit coverage")
 
+        # Determine number of parts to break ray into (this is what gets integrated over)
+        nParts = np.ceil(ray_lengths.max((1,2)) / MAX_SEGMENT_LENGTH).astype(int) + 1
+
         # iterate over weather model height levels
         for zz, nparts in enumerate(nParts):
             fracs = np.linspace(0., 1., num=nparts)

tools/RAiDER/losreader.py

@@ -748,7 +748,7 @@ def build_ray(model_zs, ht, xyz, LOS, MAX_TROPO_HEIGHT=_ZREF):
 
         # If high_ht < height of point - no contribution to integral
         # If low_ht > max_tropo_height - no contribution to integral
-        if (high_ht <= ht) or (low_ht >= MAX_TROPO_HEIGHT):
+        if (high_ht < ht) or (low_ht >= MAX_TROPO_HEIGHT):
             continue
 
         # If low_ht < requested height, start integral at requested height
@@ -783,5 +783,8 @@ def build_ray(model_zs, ht, xyz, LOS, MAX_TROPO_HEIGHT=_ZREF):
         low_xyzs.append(low_xyz)
         high_xyzs.append(high_xyz)
 
-    return np.stack(ray_lengths), np.stack(low_xyzs), np.stack(high_xyzs)
-
+    ## if all weather model levels are requested the top most layer might not contribute anything
+    if not ray_lengths:
+        return None, None, None
+    else:
+        return np.stack(ray_lengths), np.stack(low_xyzs), np.stack(high_xyzs)

tools/RAiDER/models/ecmwf.py

@@ -41,26 +41,13 @@ def __init__(self):
         self._model_level_type = 'ml'  # Default
 
 
-    def setLevelType(self, levelType):
-        '''Set the level type to model levels or pressure levels'''
-        if levelType in ['ml', 'pl']:
-            self._model_level_type = levelType
-        else:
-            raise RuntimeError('Level type {} is not recognized'.format(levelType))
-
-        if levelType == 'ml':
-            self.__model_levels__()
-        else:
-            self.__pressure_levels__()
-
-
     def __pressure_levels__(self):
         self._zlevels = np.flipud(LEVELS_25_HEIGHTS)
-        self._levels = len(self._zlevels)
+        self._levels  = len(self._zlevels)
 
 
     def __model_levels__(self):
-        self._levels = 137
+        self._levels  = 137
         self._zlevels = np.flipud(LEVELS_137_HEIGHTS)
         self._a = A_137_HRES
         self._b = B_137_HRES
@@ -321,39 +308,24 @@ def _load_pressure_level(self, filename, *args, **kwargs):
         self._t = t
         self._q = q
 
-        geo_hgt = z / self._g0
+        geo_hgt = (z / self._g0).transpose(1, 2, 0)
 
         # re-assign lons, lats to match heights
         self._lons, self._lats = np.meshgrid(lons, lats)
 
         # correct heights for latitude
-        self._get_heights(_lats, geo_hgt)
+        self._get_heights(self._lats, geo_hgt)
 
-        self._p = np.broadcast_to(levels[:, np.newaxis, np.newaxis],
+        self._p = np.broadcast_to(levels[np.newaxis, np.newaxis, :],
                                   self._zs.shape)
 
-        # Re-structure everything from (heights, lats, lons) to (lons, lats, heights)
-        self._p = np.transpose(self._p)
-        self._t = np.transpose(self._t)
-        self._q = np.transpose(self._q)
-        self._lats = np.transpose(_lats)
-        self._lons = np.transpose(_lons)
+        # Re-structure from (heights, lats, lons) to (lons, lats, heights)
+        self._t = self._t.transpose(1, 2, 0)
+        self._q = self._q.transpose(1, 2, 0)
         self._ys = self._lats.copy()
         self._xs = self._lons.copy()
-        self._zs = np.transpose(self._zs)
-
-        # check this
-        # data cube format should be lats,lons,heights
-        self._lats = self._lats.swapaxes(0, 1)
-        self._lons = self._lons.swapaxes(0, 1)
-        self._xs = self._xs.swapaxes(0, 1)
-        self._ys = self._ys.swapaxes(0, 1)
-        self._zs = self._zs.swapaxes(0, 1)
-        self._p = self._p.swapaxes(0, 1)
-        self._q = self._q.swapaxes(0, 1)
-        self._t = self._t.swapaxes(0, 1)
-
-        # For some reason z is opposite the others
+
+        # flip z to go from surface to toa
         self._p = np.flip(self._p, axis=2)
         self._t = np.flip(self._t, axis=2)
         self._q = np.flip(self._q, axis=2)

tools/RAiDER/models/hrrr.py

@@ -12,16 +12,20 @@
 from shapely.geometry import Polygon, box
 
 from RAiDER.utilFcns import round_date, transform_coords, rio_profile, rio_stats
-from RAiDER.models.weatherModel import (
-    WeatherModel, TIME_RES
-)
-from RAiDER.models.model_levels import (
-    LEVELS_137_HEIGHTS,
-)
+from RAiDER.models.weatherModel import WeatherModel, TIME_RES
+from RAiDER.models.model_levels import LEVELS_50_HEIGHTS
 from RAiDER.logger import logger
 
 
-def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='prs', fxx=0, verbose=False):
+HRRR_CONUS_COVERAGE_POLYGON = Polygon(((-125, 21), (-133, 49), (-60, 49), (-72, 21)))
+HRRR_AK_COVERAGE_POLYGON = Polygon(((195, 40), (157, 55), (175, 70), (260, 77), (232, 52)))
+HRRR_AK_PROJ = CRS.from_string('+proj=stere +ellps=sphere +a=6371229.0 +b=6371229.0 +lat_0=90 +lon_0=225.0 '
+                               '+x_0=0.0 +y_0=0.0 +lat_ts=60.0 +no_defs +type=crs')
+# Source: https://eric.clst.org/tech/usgeojson/
+AK_GEO = gpd.read_file(Path(__file__).parent / 'data' / 'alaska.geojson.zip').geometry.unary_union
+
+
+def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='nat', fxx=0, verbose=False):
     '''
     Download a HRRR weather model using Herbie
 
@@ -47,6 +51,7 @@ def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='prs', fxx=0,
         save_dir=Path(os.path.dirname(out)),
     )
 
+
     # Iterate through the list of datasets
     try:
         ds_list = H.xarray(":(SPFH|PRES|TMP|HGT):", verbose=verbose)
@@ -55,9 +60,15 @@ def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='prs', fxx=0,
         raise ValueError
 
     ds_out = None
+
     for ds in ds_list:
-        if ('isobaricInhPa' in ds._coord_names) or ('levels' in ds._coord_names):
+        if 'isobaricInhPa' in ds._coord_names:
+            ds_out = ds
+            coord = 'isobaricInhPa'
+            break
+        elif 'hybrid' in ds._coord_names:
             ds_out = ds
+            coord = 'hybrid'
             break
 
     # subset the full file by AOI
@@ -68,7 +79,7 @@ def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='prs', fxx=0,
     )
 
     # bookkeepping
-    ds_out = ds_out.rename({'gh': 'z', 'isobaricInhPa': 'levels'})
+    ds_out = ds_out.rename({'gh': 'z', coord: 'levels'})
     ny, nx = ds_out['longitude'].shape
 
     # projection information
@@ -145,9 +156,9 @@ def load_weather_hrrr(filename):
     '''
     # read data from the netcdf file
     ds = xarray.open_dataset(filename, engine='netcdf4')
-
     # Pull the relevant data from the file
-    pl = np.array([p * 100 for p in ds.levels.values]) # convert millibars to Pascals
+    pl = ds.levels.values
+    pres = ds['pres'].values.transpose(1, 2, 0)
     xArr = ds['x'].values
     yArr = ds['y'].values
     lats = ds['latitude'].values
@@ -166,14 +177,7 @@ def load_weather_hrrr(filename):
     _ys = np.broadcast_to(yArr[:, np.newaxis, np.newaxis],
                             geo_hgt.shape)
 
-    return _xs, _ys, lons, lats, qs, temps, pl, geo_hgt, proj
-
-HRRR_CONUS_COVERAGE_POLYGON = Polygon(((-125, 21), (-133, 49), (-60, 49), (-72, 21)))
-HRRR_AK_COVERAGE_POLYGON = Polygon(((195, 40), (157, 55), (175, 70), (260, 77), (232, 52)))
-HRRR_AK_PROJ = CRS.from_string('+proj=stere +ellps=sphere +a=6371229.0 +b=6371229.0 +lat_0=90 +lon_0=225.0 '
-                               '+x_0=0.0 +y_0=0.0 +lat_ts=60.0 +no_defs +type=crs')
-# Source: https://eric.clst.org/tech/usgeojson/
-AK_GEO = gpd.read_file(Path(__file__).parent / 'data' / 'alaska.geojson.zip').geometry.unary_union
+    return _xs, _ys, lons, lats, qs, temps, pres, geo_hgt, proj
 
 
 class HRRR(WeatherModel):
@@ -209,7 +213,6 @@ def __init__(self):
         self._Npl = 0
         self.files = None
         self._bounds = None
-        self._zlevels = np.flipud(LEVELS_137_HEIGHTS)
 
         # Projection
         # NOTE: The HRRR projection will get read directly from the downloaded weather model file; however,
@@ -228,12 +231,20 @@ def __init__(self):
         x0 = 0
         y0 = 0
         earth_radius = 6371229
-        p1 = CRS(f'+proj=lcc +lat_1={lat1} +lat_2={lat2} +lat_0={lat0} '\
+        self._proj = CRS(f'+proj=lcc +lat_1={lat1} +lat_2={lat2} +lat_0={lat0} '\
                  f'+lon_0={lon0} +x_0={x0} +y_0={y0} +a={earth_radius} '\
                  f'+b={earth_radius} +units=m +no_defs')
-        self._proj = p1
-
         self._valid_bounds = HRRR_CONUS_COVERAGE_POLYGON
+        self.setLevelType('nat')
+
+
+    def __model_levels__(self):
+        self._levels  = 50
+        self._zlevels = np.flipud(LEVELS_50_HEIGHTS)
+
+
+    def __pressure_levels__(self):
+        raise NotImplementedError('Pressure levels do not go high enough for HRRR.')
 
 
     def _fetch(self,  out):
@@ -249,7 +260,8 @@ def _fetch(self,  out):
         # HRRR uses 0-360 longitude, so we need to convert the bounds to that
         bounds = self._ll_bounds.copy()
         bounds[2:] = np.mod(bounds[2:], 360)
-        download_hrrr_file(bounds, corrected_DT, out, model='hrrr')
+
+        download_hrrr_file(bounds, corrected_DT, out, 'hrrr', self._model_level_type)
 
 
     def load_weather(self, f=None, *args, **kwargs):
@@ -261,18 +273,18 @@ def load_weather(self, f=None, *args, **kwargs):
             f = self.files[0] if isinstance(self.files, list) else self.files
 
 
-        _xs, _ys, _lons, _lats, qs, temps, pl, geo_hgt, proj = load_weather_hrrr(f)
-            # correct for latitude
+        _xs, _ys, _lons, _lats, qs, temps, pres, geo_hgt, proj = load_weather_hrrr(f)
+
+        # convert geopotential height to geometric height
         self._get_heights(_lats, geo_hgt)
 
         self._t = temps
         self._q = qs
-        self._p = np.broadcast_to(pl[np.newaxis, np.newaxis, :], geo_hgt.shape)
+        self._p = pres
         self._xs = _xs
         self._ys = _ys
         self._lats = _lats
         self._lons = _lons
-
         self._proj = proj
 
 
@@ -333,7 +345,6 @@ def __init__(self):
         self._Npl = 0
         self.files = None
         self._bounds = None
-        self._zlevels = np.flipud(LEVELS_137_HEIGHTS)
 
         self._classname = 'hrrrak'
         self._dataset = 'hrrrak'
@@ -345,6 +356,17 @@ def __init__(self):
         # The projection information gets read directly from the  weather model file but we
         # keep this here for object instantiation.
         self._proj = HRRR_AK_PROJ
+        self.setLevelType('nat')
+
+
+    def __model_levels__(self):
+        self._levels  = 50
+        self._zlevels = np.flipud(LEVELS_50_HEIGHTS)
+
+
+    def __pressure_levels__(self):
+        raise NotImplementedError('hrrr.py: Revisit whether or not pressure levels from HRRR can be used for delay calculations; they do not go high enough compared to native model levels.')
+
 
     def _fetch(self, out):
         bounds = self._ll_bounds.copy()
@@ -354,19 +376,20 @@ def _fetch(self, out):
         if not corrected_DT == self._time:
             logger.info('Rounded given datetime from {} to {}'.format(self._time, corrected_DT))
 
-        download_hrrr_file(bounds, corrected_DT, out, model='hrrrak')
+        download_hrrr_file(bounds, corrected_DT, out, 'hrrrak', self._model_level_type)
 
 
     def load_weather(self, f=None, *args, **kwargs):
         if f is None:
             f = self.files[0] if isinstance(self.files, list) else self.files
-        _xs, _ys, _lons, _lats, qs, temps, pl, geo_hgt, proj = load_weather_hrrr(f)
-            # correct for latitude
+        _xs, _ys, _lons, _lats, qs, temps, pres, geo_hgt, proj = load_weather_hrrr(f)
+
+        # correct for latitude
         self._get_heights(_lats, geo_hgt)
 
         self._t = temps
         self._q = qs
-        self._p = np.broadcast_to(pl[np.newaxis, np.newaxis, :], geo_hgt.shape)
+        self._p = pres
         self._xs = _xs
         self._ys = _ys
         self._lats = _lats

tools/RAiDER/models/model_levels.py

@@ -505,3 +505,25 @@
     -300,
     -500,
 ]
+
+## HRRR Model Levels
+# Computed according to: H = a + b * Z where:
+    # H is the resulting levels in geometric height
+    # a is the Surface geopotential height (in meters)
+        # the surface geopotential is fetched from Herbie and divided by g0=9.80665
+        # averaged in space over CONUS
+    # b is the native (sigma) model levels (https://rapidrefresh.noaa.gov/faq/HRRR.faq.html)
+    # Z is the spatial average geopotential height of the sigma level (in meters)
+LEVELS_50_HEIGHTS = [2.61580385e+04, 2.48712879e+04, 2.36910518e+04, 2.25524744e+04,
+                      2.13986900e+04, 2.02464207e+04, 1.90883153e+04, 1.79427740e+04,
+                      1.68476065e+04, 1.57399654e+04, 1.45826790e+04, 1.33886515e+04,
+                      1.22171878e+04, 1.11019360e+04, 1.00395775e+04, 9.01965365e+03,
+                      8.03486128e+03, 7.09323111e+03, 6.27822334e+03, 5.57101666e+03,
+                      4.96120000e+03, 4.42159162e+03, 3.94118518e+03, 3.51064883e+03,
+                      3.12371808e+03, 2.77490670e+03, 2.45941860e+03, 2.17290722e+03,
+                      1.90394551e+03, 1.66716448e+03, 1.44127808e+03, 1.22697117e+03,
+                      1.02507126e+03, 8.38877887e+02, 6.74297597e+02, 5.34810131e+02,
+                      4.18916771e+02, 3.23291544e+02, 2.44985788e+02, 1.81492083e+02,
+                      1.34383211e+02, 1.02007390e+02, 7.70762881e+01, 5.77739913e+01,
+                      4.31591299e+01, 3.26389095e+01, 2.52657431e+01, 2.02104423e+01,
+                      1.66520787e+01, 1.39366382e+01, 0, -10, -20, -50, -100, -200, -500]