Merge pull request #719 from jlmaurer/fix_hrrr

Use known CRS for HRRR-AK instead of reading CRS from file

Author: jlmaurer

CHANGELOG.md

@@ -15,6 +15,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 ### Fixed
 * [721](https://github.com/dbekaert/RAiDER/pull/721) - Fixed bug mixing station_delay_mean and station_delay_median
 * [715](https://github.com/dbekaert/RAiDER/pull/715) - Fixed the coverage test Github action and a timing issue with raiderCombine
+* [719](https://github.com/dbekaert/RAiDER/pull/719) - Fixed the automatic switch to HRRR-AK based on bounding box when generic HRRR is specified as the weather model
 * [720](https://github.com/dbekaert/RAiDER/pull/720) - Bug-fix to properly pass expected `user_title` argument through raiderStats function.
 
 ## [0.5.4]

test/gunw_test_data/weather_files/GMAO_2020_01_24T13_52_44_timeInterp_32N_36N_121W_114W.nc

@@ -22,9 +22,10 @@
 from RAiDER.aria.prepFromGUNW import (
     check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation,
     check_weather_model_availability,_get_acq_time_from_gunw_id,
-    get_slc_ids_from_gunw,get_acq_time_from_slc_id
+    get_slc_ids_from_gunw,get_acq_time_from_slc_id,identify_which_hrrr
 )
 from RAiDER.cli.raider import calcDelaysGUNW
+from RAiDER.models.hrrr import HRRR, HRRRAK
 from RAiDER.models.customExceptions import (
      NoWeatherModelData, WrongNumberOfFiles,
 ) 
@@ -632,3 +633,72 @@ def test_get_acq_time_invalid_slc_id():
   invalid_slc_id = "test/gunw_azimuth_test_data/S1B_OPER_AUX_POEORB_OPOD_20210731T111940_V20210710T225942_20210712T005942.EOF"
   with pytest.raises(ValueError):
     get_acq_time_from_slc_id(invalid_slc_id)
+
+
+def test_identify_which_hrrr_1():
+    """Tests if function identifies the correct HRRR file"""
+    gunw_id = Path("test/gunw_azimuth_test_data/S1-GUNW-A-R-064-tops-20210723_20210711-015000-00119W_00033N-PP-6267-v2_0_6.nc")
+    result = identify_which_hrrr(gunw_id)
+    assert result == "HRRR"
+
+
+def test_identify_which_hrrr_2():
+    """Tests if function identifies the correct HRRR file"""
+    gunw_id = Path("test/gunw_test_data/S1-GUNW-D-R-059-tops-20230320_20220418-180300-00179W_00051N-PP-c92e-v2_0_6.nc")
+    result = identify_which_hrrr(gunw_id)
+    assert result == "HRRRAK"
+
+
+def test_cast_to_hrrrak_1():
+    """Tests if function casts the HRRR file to HRRRAK"""
+    ak_bounds = [51.0, 71.0, -175., -130.0]
+    conus_bounds = [34.0,35.0, -91,  -90.0]
+    model = HRRR()
+    model.checkValidBounds(conus_bounds)
+    model.checkValidBounds(ak_bounds)
+    assert model._Name == "HRRR-AK"
+
+
+def test_cast_to_hrrrak_2():
+    """Tests if function casts the HRRR file to HRRRAK"""
+    ak_bounds = [51.0, 71.0, -175., -130.0]
+    model = HRRRAK()
+    model.checkValidBounds(ak_bounds)
+    assert model._Name == "HRRR-AK"
+
+
+def test_cast_to_hrrrak_2b():
+    """Tests if function casts the HRRR file to HRRRAK"""
+    ak_bounds = [60.0, 65.0, -150., -120.0]
+    model = HRRRAK()
+    model.checkValidBounds(ak_bounds)
+    assert model._Name == "HRRR-AK"
+
+
+def test_cast_to_hrrrak_3():
+    """Tests if function casts the HRRR file to HRRRAK"""
+    conus_bounds = [34.0,35.0, -91,  -90.0]
+    model = HRRR()
+    model.checkValidBounds(conus_bounds)
+    assert model._Name == "HRRR"
+
+
+def test_cast_to_hrrrak_4():
+    """Tests if function casts the HRRR file to HRRRAK"""
+    europe_bounds = [-1, 1, -1, 1]
+    model = HRRR()
+    with pytest.raises(ValueError):
+        model.checkValidBounds(europe_bounds)
+
+
+def test_identify_which_hrrr_invalid():
+    """Tests if function raises error for an invalid gunw_id format"""
+    invalid_gunw_id = "dummy.nc"
+    with pytest.raises(NoWeatherModelData):
+        identify_which_hrrr(invalid_gunw_id)
+
+
+def test_check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation_again():
+    """Tests if function raises error for an invalid gunw_id format"""
+    gunw_id = "S1-GUNW-D-R-044-tops-20240418_20240406-171649-00163W_00069N-PP-af6b-v3_0_1.nc"
+    assert check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id, 'hrrrak') is True

test/gunw_test_data/weather_files/GMAO_2020_01_30T13_52_44_timeInterp_32N_36N_121W_114W.nc

@@ -23,7 +23,9 @@
 from RAiDER.models.gmao import GMAO
 from RAiDER.models.merra2 import MERRA2
 from RAiDER.models.ncmr import NCMR
-from RAiDER.models.customExceptions import DatetimeOutsideRange
+from RAiDER.models.customExceptions import (
+    DatetimeOutsideRange, NoWeatherModelData,
+)
 
 
 _LON0 = 0
@@ -435,7 +437,7 @@ def test_get_bounds_indices_3() -> None:
     l = np.arange(-20, 20)
     l2 = (((np.arange(160, 200) + 180) % 360) - 180)
     lats, lons = np.meshgrid(l, l2)
-    with pytest.raises(RuntimeError):
+    with pytest.raises(NoWeatherModelData):
         get_bounds_indices(snwe, lats, lons)
 
 

test/test_GUNW.py

@@ -20,6 +20,7 @@
 
 from RAiDER.logger import logger
 from RAiDER.models import credentials
+from RAiDER.models.customExceptions import NoWeatherModelData
 from RAiDER.models.hrrr import AK_GEO, HRRR_CONUS_COVERAGE_POLYGON, check_hrrr_dataset_availability
 from RAiDER.s1_azimuth_timing import get_times_for_azimuth_interpolation
 from RAiDER.s1_orbits import get_orbits_from_slc_ids
@@ -50,7 +51,7 @@ def _get_acq_time_from_gunw_id(gunw_id: str, reference_or_secondary: str) -> dt.
     return cen_acq_time
 
 
-def check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id: str) -> bool:
+def check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id: str, weather_model_name: str='hrrr') -> bool:
     """
     Determine if all the times for azimuth interpolation are available using
     Herbie. Note that not all 1 hour times are available within the said date
@@ -71,11 +72,11 @@ def check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id: st
     ref_acq_time = _get_acq_time_from_gunw_id(gunw_id, 'reference')
     sec_acq_time = _get_acq_time_from_gunw_id(gunw_id, 'secondary')
 
-    model_step_hours = 1
+    model_step_hours = [1 if weather_model_name == 'hrrr' else 3][0]
     ref_times_for_interp = get_times_for_azimuth_interpolation(ref_acq_time, model_step_hours)
     sec_times_for_interp = get_times_for_azimuth_interpolation(sec_acq_time, model_step_hours)
-    ref_dataset_availability = list(map(check_hrrr_dataset_availability, ref_times_for_interp))
-    sec_dataset_availability = list(map(check_hrrr_dataset_availability, sec_times_for_interp))
+    ref_dataset_availability = list(map(check_hrrr_dataset_availability, ref_times_for_interp, [weather_model_name]*len(ref_times_for_interp)))
+    sec_dataset_availability = list(map(check_hrrr_dataset_availability, sec_times_for_interp, [weather_model_name]*len(sec_times_for_interp)))
 
     return all(ref_dataset_availability) and all(sec_dataset_availability)
 
@@ -125,14 +126,9 @@ def check_weather_model_availability(gunw_path: Path, weather_model_name: str) -
     sec_ts = get_acq_time_from_slc_id(sec_slc_ids[0]).replace(tzinfo=dt.timezone(offset=dt.timedelta()))
 
     if weather_model_name == 'HRRR':
-        group = '/science/grids/data/'
-        with xr.open_dataset(gunw_path, group=f'{group}') as ds:
-            gunw_poly = box(*ds.rio.bounds())
-        if HRRR_CONUS_COVERAGE_POLYGON.intersects(gunw_poly):
-            pass
-        elif AK_GEO.intersects(gunw_poly):
-            weather_model_name = 'HRRRAK'
-        else:
+        try:
+            weather_model_name = identify_which_hrrr(gunw_path)
+        except NoWeatherModelData:
             return False
 
     # source: https://stackoverflow.com/a/7668273
@@ -387,3 +383,26 @@ def main(args: CalcDelaysArgs) -> tuple[Path, float]:
     path_cfg = Path(f'GUNW_{GUNWObj.name}.yaml')
     write_yaml(raider_cfg, path_cfg)
     return path_cfg, GUNWObj.wavelength
+
+
+def identify_which_hrrr(gunw_path: Path) -> str:
+    group = '/science/grids/data/'
+    try:
+        with xr.open_dataset(gunw_path, group=f'{group}') as ds:
+            gunw_poly = box(*ds.rio.bounds())
+        if HRRR_CONUS_COVERAGE_POLYGON.intersects(gunw_poly):
+            weather_model_name = 'HRRR'
+        elif AK_GEO.intersects(gunw_poly):
+            weather_model_name = 'HRRRAK'
+        else:
+            raise NoWeatherModelData(
+                f'GUNW {gunw_path} does not intersect with any HRRR coverage area. '
+                'Please use a different weather model.'
+            )
+    except FileNotFoundError:
+        raise NoWeatherModelData(
+            f'''GUNW {gunw_path} does not exist or is not a valid HRRR file.
+            Please check the file path.
+            '''
+        )
+    return weather_model_name

test/test_weather_model.py

@@ -30,6 +30,7 @@
     RuntimeGroup,
     TimeGroup,
 )
+from RAiDER.aria.prepFromGUNW import identify_which_hrrr
 from RAiDER.cli.validators import DateListAction, date_type
 from RAiDER.gnss.types import RAiDERCombineArgs
 from RAiDER.logger import logger, logging
@@ -334,7 +335,7 @@ def calcDelays(iargs: Optional[Sequence[str]]=None) -> list[Path]:
 
         if len(wfiles) == 0:
             logger.error('No weather model data was successfully processed.')
-            raise NoWeatherModelData()
+            raise NoWeatherModelData('Weather model processing failed for all times')
 
         # Get the weather model file
         weather_model_file = getWeatherFile(wfiles, times, t, model._Name, interp_method)
@@ -604,7 +605,8 @@ def calcDelaysGUNW(iargs: Optional[list[str]] = None) -> Optional[xr.Dataset]:
         args.interpolate_time == 'azimuth_time_grid'
     ):
         gunw_id = args.file.name.replace('.nc', '')
-        if not RAiDER.aria.prepFromGUNW.check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id):
+        weather_model_name = identify_which_hrrr(args.file)
+        if not RAiDER.aria.prepFromGUNW.check_hrrr_dataset_availablity_for_s1_azimuth_time_interpolation(gunw_id, weather_model_name):
             raise NoWeatherModelData('The required HRRR data for time-grid interpolation is not available')
 
     if args.file is None:

tools/RAiDER/aria/prepFromGUNW.py

@@ -8,8 +8,10 @@
 from herbie import Herbie
 from pyproj import CRS, Transformer
 from shapely.geometry import Polygon, box
+from typing import Optional, Union, List, Tuple
 
 from RAiDER.logger import logger
+from RAiDER.models.customExceptions import NoWeatherModelData
 from RAiDER.models.model_levels import LEVELS_50_HEIGHTS
 from RAiDER.models.weatherModel import TIME_RES, WeatherModel
 from RAiDER.utilFcns import round_date
@@ -25,11 +27,11 @@
 AK_GEO = gpd.read_file(Path(__file__).parent / 'data' / 'alaska.geojson.zip').geometry.unary_union
 
 
-def check_hrrr_dataset_availability(datetime: dt.datetime) -> bool:
+def check_hrrr_dataset_availability(datetime: dt.datetime, model='hrrr') -> bool:
     """Note a file could still be missing within the models valid range."""
     herbie = Herbie(
         datetime,
-        model='hrrr',
+        model=model,
         product='nat',
         fxx=0,
     )
@@ -87,11 +89,17 @@ def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='nat', fxx=0,
         raise RuntimeError('Herbie did not obtain an HRRR dataset with the expected layers and coordinates')
 
     # subset the full file by AOI
-    x_min, x_max, y_min, y_max = get_bounds_indices(
-        ll_bounds,
-        ds_out.latitude.to_numpy(),
-        ds_out.longitude.to_numpy(),
-    )
+    try:
+        x_min, x_max, y_min, y_max = get_bounds_indices(
+            ll_bounds,
+            ds_out.latitude.to_numpy(),
+            ds_out.longitude.to_numpy(),
+        )
+    except NoWeatherModelData as e:
+        logger.error(e)
+        logger.error('lat/lon bounds: %s', ll_bounds)
+        logger.error('Weather model is {}'.format(model))
+        raise
 
     # bookkeepping
     ds_out = ds_out.rename({'gh': 'z', coord: 'levels'})
@@ -104,8 +112,15 @@ def download_hrrr_file(ll_bounds, DATE, out, model='hrrr', product='nat', fxx=0,
         ds_out[var].attrs['grid_mapping'] = 'proj'
 
     # pull the grid information
-    proj = CRS.from_cf(ds_out['proj'].attrs)
+    if model=='hrrr':
+        proj = CRS.from_cf(ds_out['proj'].attrs)
+    elif model=='hrrrak':
+        proj = HRRR_AK_PROJ
+    else:
+        raise ValueError('Model not recognized. Please use either hrrr or hrrrak')
+
     t = Transformer.from_crs(4326, proj, always_xy=True)
+
     xl, yl = t.transform(ds_out['longitude'].values, ds_out['latitude'].values)
     W, E, S, N = np.nanmin(xl), np.nanmax(xl), np.nanmin(yl), np.nanmax(yl)
 
@@ -139,7 +154,7 @@ def get_bounds_indices(SNWE, lats, lons):
         W, E = np.mod([W, E], 360)
         m1 = (S <= lats) & (N >= lats) & (W <= lons) & (E >= lons)
         if np.sum(m1) == 0:
-            raise RuntimeError('Area of Interest has no overlap with the HRRR model available extent')
+            raise NoWeatherModelData('Area of Interest has no overlap with the HRRR model available extent')
 
     # Y extent
     shp = lats.shape
@@ -266,6 +281,23 @@ def _fetch(self, out) -> None:
 
         download_hrrr_file(bounds, corrected_DT, out, 'hrrr', self._model_level_type)
 
+    def _cast_to_hrrrak(self) -> None:
+        """
+        Update the weather model to HRRR-AK if the bounding box is in Alaska.
+        """
+        self.__class__ = HRRRAK
+        self._dataset = 'hrrrak'
+        self._valid_bounds = HRRR_AK_COVERAGE_POLYGON
+        self._proj = HRRR_AK_PROJ
+        self._Name = 'HRRR-AK'
+        self._time_res = TIME_RES['HRRR-AK']
+        self._valid_range = (
+            dt.datetime(2018, 7, 13).replace(tzinfo=dt.timezone(offset=dt.timedelta())),
+            dt.datetime.now(dt.timezone.utc),
+        )
+        self.setLevelType('nat')
+        
+
     def load_weather(self, f=None, *args, **kwargs) -> None:
         """
         Load a weather model into a python weatherModel object, from self.files if no
@@ -294,7 +326,7 @@ def checkValidBounds(self, ll_bounds: np.ndarray) -> None:
         (i.e., intersects with the model domain at all).
 
         Args:
-        ll_bounds : np.ndarray
+        ll_bounds : np.ndarray SNWE bounding box
         """
         S, N, W, E = ll_bounds
         aoi = box(W, S, E, N)
@@ -306,13 +338,18 @@ def checkValidBounds(self, ll_bounds: np.ndarray) -> None:
             logger.critical('The HRRR weather model extent does not completely cover your AOI!')
 
         else:
+            logger.info('The HRRR weather model extent does not include your AOI!')
+            logger.info('Checking the HRRR-AK model.')
             Mod = HRRRAK()
             # valid bounds are in 0->360 to account for dateline crossing
             W, E = np.mod([W, E], 360)
             aoi = box(W, S, E, N)
             if Mod._valid_bounds.contains(aoi):
-                pass
+                self._cast_to_hrrrak()
+                logger.info('Casting self to the HRRR-AK weather model.')
             elif aoi.intersects(Mod._valid_bounds):
+                self._cast_to_hrrrak()
+                logger.info('Casting self to the HRRR-AK weather model.')
                 logger.critical('The HRRR-AK weather model extent does not completely cover your AOI!')
 
             else:
@@ -391,3 +428,4 @@ def load_weather(self, f=None, *args, **kwargs) -> None:
         self._lats = _lats
         self._lons = _lons
         self._proj = proj
+

tools/RAiDER/cli/raider.py

@@ -465,7 +465,9 @@ def checkValidBounds(
         """
         S, N, W, E = ll_bounds
         if not box(W, S, E, N).intersects(self._valid_bounds):
-            raise ValueError(f'The requested location is unavailable for {self._Name}')
+            #check the longitude wrapping
+            if not box(360+W, S, 360+E, N).intersects(self._valid_bounds):
+                raise ValueError(f'The requested location is unavailable for {self._Name}')
 
     def checkContainment(self, ll_bounds: Union[List, Tuple,np.ndarray], buffer_deg: float = 1e-5) -> bool:
         """

tools/RAiDER/models/hrrr.py

@@ -89,7 +89,7 @@ def prepareWeatherModel(
         logger.warning('The processed weather model file already exists, so I will use that.')
 
         containment = weather_model.checkContainment(ll_bounds)
-        if not containment and weather_model.Model() not in 'HRRR'.split():
+        if not containment and weather_model.Model() not in 'HRRR HRRRAK'.split():
             raise ExistingWeatherModelTooSmall
 
         return f