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_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.2600022, 0.0358646 # 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/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('Pressure levels do not go high enough for HRRR.')
+
 
     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,17 @@
     -300,
     -500,
 ]
+
+## HRRR Model Levels
+# these are computed as the spatial (over CONUS) average of the HRRR pressure field
+LEVELS_50_HEIGHTS = [26143.662, 24907.1, 23820.475, 22836.285, 21934.553,
+                    21099.717, 20316.678, 19577.092, 18876.715, 18185.781,
+                    17453.205, 16682.344, 15936.782, 15257.467, 14641.991,
+                    14077.848, 13555.656, 13068.57, 12611.642, 12181.063,
+                    11773.499, 11377.56, 10973.201, 10551.642, 10113.527,
+                    9654.279, 9170.161, 8661.178, 8127.03, 7568.2935,
+                    6985.858, 6378.0234, 5743.416, 5095.357, 4469.7725,
+                    3899.744, 3392.3884, 2941.3154, 2540.2717, 2184.2336,
+                    1869.1932, 1591.8563, 1347.2079, 1128.5986, 935.72327,
+                    775.0872, 651.92017, 567.2376, 515.85925, 486.44925,
+                    200, 100, 0, -10, -20, -50, -100, -200, -500]

tools/RAiDER/models/weatherModel.py

@@ -326,6 +326,19 @@ def checkTime(self, time):
             raise RuntimeError(msg)
 
 
+    def setLevelType(self, levelType):
+        '''Set the level type to model levels or pressure levels'''
+        if levelType in 'ml pl nat prs'.split():
+            self._model_level_type = levelType
+        else:
+            raise RuntimeError(f'Level type {levelType} is not recognized')
+
+        if levelType in 'ml nat'.split():
+            self.__model_levels__()
+        else:
+            self.__pressure_levels__()
+
+
     def _convertmb2Pa(self, pres):
         '''
         Convert pressure in millibars to Pascals
@@ -860,7 +873,6 @@ def get_mapping(proj):
         return proj.to_wkt()
 
 
-
 def checkContainment_raw(path_wm_raw,
                         ll_bounds,
                         buffer_deg: float = 1e-5) -> bool: