mock weather model files

test/conftest.py

@@ -61,14 +61,39 @@ def gunw_azimuth_test():
 
 
 @pytest.fixture(scope='session')
-def orbit_dict_for_azimuth_test():
+def orbit_dict_for_azimuth_time_test():
     test_data = TEST_DIR / 'gunw_azimuth_test_data'
     return {'reference': test_data / 'S1B_OPER_AUX_POEORB_OPOD_20210812T111941_V20210722T225942_20210724T005942.EOF',
             'secondary': test_data / 'S1B_OPER_AUX_POEORB_OPOD_20210731T111940_V20210710T225942_20210712T005942.EOF'}
 
 
 @pytest.fixture(scope='session')
-def slc_id_dict_for_azimuth_test():
+def slc_id_dict_for_azimuth_time_test():
     test_data = TEST_DIR / 'gunw_azimuth_test_data'
     return {'reference': test_data / 'S1B_IW_SLC__1SDV_20210723T014947_20210723T015014_027915_0354B4_B3A9',
             'secondary': test_data / 'S1B_IW_SLC__1SDV_20210711T014947_20210711T015013_027740_034F80_D404'}
+
+
+@pytest.fixture(scope='session')
+def weather_model_dict_for_azimuth_time_test():
+    """The order is important; will be closest to InSAR acq time so goes 2, 1, 3 AM."""
+    test_data = TEST_DIR / 'gunw_azimuth_test_data' / 'weather_files'
+    return {'HRRR': [test_data / 'HRRR_2021_07_23_T02_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_23_T01_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_23_T03_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_11_T02_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_11_T01_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_11_T03_00_00_33N_36N_120W_115W.nc'
+                     ]}
+
+
+@pytest.fixture(scope='session')
+def weather_model_dict_for_center_time_test():
+    """Order is important here; will be in chronological order with respect to closest date times"""
+    test_data = TEST_DIR / 'gunw_azimuth_test_data' / 'weather_files'
+    return {'HRRR': [test_data / 'HRRR_2021_07_23_T01_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_23_T02_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_11_T01_00_00_33N_36N_120W_115W.nc',
+                     test_data / 'HRRR_2021_07_11_T02_00_00_33N_36N_120W_115W.nc',
+                     ]
+            }

test/test_GUNW.py

@@ -4,13 +4,12 @@
 import shutil
 import subprocess
 import unittest
+from pathlib import Path
 
+import hyp3lib
 import jsonschema
-<<<<<<< HEAD
 import numpy as np
-=======
 import pandas as pd
->>>>>>> dev
 import pytest
 import rasterio as rio
 import xarray as xr
@@ -24,9 +23,9 @@
 
 def compute_transform(lats, lons):
     """ Hand roll an affine transform from lat/lon coords """
-    a = lons[1] - lons[0] # lon spacing
+    a = lons[1] - lons[0]  # lon spacing
     b = 0
-    c = lons[0] - a/2 # lon start, adjusted by half a grid cell
+    c = lons[0] - a/2  # lon start, adjusted by half a grid cell
     d = 0
     e = lats[1] - lats[0]
     f = lats[0] - e/2
@@ -55,7 +54,7 @@ def test_GUNW_update(test_dir_path, test_gunw_path_factory, weather_model_name):
         for v in 'troposphereWet troposphereHydrostatic'.split():
             da = ds[v]
             lats, lons = da.latitudeMeta.to_numpy(), da.longitudeMeta.to_numpy()
-            transform  = compute_transform(lats, lons)
+            transform = compute_transform(lats, lons)
             assert da.rio.transform().almost_equals(transform), 'Affine Transform incorrect'
 
         crs = rio.crs.CRS.from_wkt(ds['crs'].crs_wkt)
@@ -124,36 +123,123 @@ def test_GUNW_metadata_update(test_gunw_json_path, test_gunw_json_schema_path, t
     ]
 
 
-@pytest.mark.parametrize('model', ['HRRR'])
-def test_azimuth_timing_against_interpolation(model, tmp_path, gunw_azimuth_test):
+@pytest.mark.parametrize('weather_model_name', ['HRRR'])
+def test_azimuth_timing_against_interpolation(weather_model_name: str,
+                                              tmp_path: Path,
+                                              gunw_azimuth_test: Path,
+                                              orbit_dict_for_azimuth_time_test: dict[str],
+                                              weather_model_dict_for_azimuth_time_test,
+                                              weather_model_dict_for_center_time_test,
+                                              mocker):
     """This test shows that the azimuth timing interpolation does not deviate from
     the center time by more than 1 mm for the HRRR model. This is expected since the model times are
     6 hours apart and a the azimuth time is changing the interpolation weights for a given pixel at the order
-    of seconds and thus these two approaches are quite similar."""
+    of seconds and thus these two approaches are quite similar.
+
+    Effectively, this mocks the following CL script
+
+    ```
+    cmd = f'raider.py ++process calcDelaysGUNW -f {out_path_0} -m {weather_model_name} -interp center_time'
+    subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+
+    cmd = f'raider.py ++process calcDelaysGUNW -f {out_path_1} -m {weather_model_name} -interp azimuth_time_grid'
+    subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+    ```
+
+    Getting weather model file names requires patience. For HRRR and center time, here are the calls:
+
+    ```
+    wfile_0 = prepareWeatherModel(model, datetime.datetime(2021, 7, 23, 1, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_1 = prepareWeatherModel(model, datetime.datetime(2021, 7, 23, 2, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_2 = prepareWeatherModel(model, datetime.datetime(2021, 7, 11, 1, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_3 = prepareWeatherModel(model, datetime.datetime(2021, 7, 11, 2, 0), [33.45, 35.45, -119.15, -115.95])
+    ```
+
+    And similarly for the azimuth time:
+
+    ```
+    wfile_0 = prepareWeatherModel(model, datetime.datetime(2021, 7, 23, 2, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_1 = prepareWeatherModel(model, datetime.datetime(2021, 7, 23, 1, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_2 = prepareWeatherModel(model, datetime.datetime(2021, 7, 23, 3, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_3 = prepareWeatherModel(model, datetime.datetime(2021, 7, 11, 2, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_4 = prepareWeatherModel(model, datetime.datetime(2021, 7, 11, 1, 0), [33.45, 35.45, -119.15, -115.95])
+    wfile_5 = prepareWeatherModel(model, datetime.datetime(2021, 7, 11, 3, 0), [33.45, 35.45, -119.15, -115.95])
+    ```
+
+    Note for azimuth time they are acquired in order of proximity to acq as opposed to chronological.
+    """
 
     out_0 = gunw_azimuth_test.name.replace('.nc', '__ct_interp.nc')
     out_1 = gunw_azimuth_test.name.replace('.nc', '__az_interp.nc')
 
     out_path_0 = shutil.copy(gunw_azimuth_test, tmp_path / out_0)
     out_path_1 = shutil.copy(gunw_azimuth_test, tmp_path / out_1)
 
-    cmd = f'raider.py ++process calcDelaysGUNW -f {out_path_0} -m HRRR -interp center_time'
-    subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
-
-    cmd = f'raider.py ++process calcDelaysGUNW -f {out_path_1} -m HRRR -interp azimuth_time_grid'
-    subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+    # In the loop here: https://github.com/dbekaert/RAiDER/blob/
+    # f77af9ce2d3875b00730603305c0e92d6c83adc2/tools/RAiDER/cli/raider.py#L233-L237
+    # Which reads the datelist from the YAML
+    # We note that reference scene is processed *first* and then secondary
+    # as yaml is created using these:
+    # https://github.com/dbekaert/RAiDER/blob/
+    # f77af9ce2d3875b00730603305c0e92d6c83adc2/tools/RAiDER/aria/prepFromGUNW.py#L151-L200
+
+    mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile',
+                 # Hyp3 Lib returns 2 values
+                 side_effect=[
+                               # For center time
+                               (orbit_dict_for_azimuth_time_test['reference'], ''),
+                               (orbit_dict_for_azimuth_time_test['secondary'], ''),
+                               # For azimuth time
+                               (orbit_dict_for_azimuth_time_test['reference'], ''),
+                               (orbit_dict_for_azimuth_time_test['secondary'], ''),
+                               ]
+                 )
+    mocker.patch('RAiDER.s1_azimuth_timing.get_slc_id_from_point_and_time',
+                 return_value=[
+                               # Center_time
+                               'reference_slc_id', 'secondary_slc_id',
+                               # Azimuth time
+                               'reference_slc_id', 'secondary_slc_id',
+                               ])
+
+    side_effect = (weather_model_dict_for_center_time_test[weather_model_name] +
+                   weather_model_dict_for_azimuth_time_test[weather_model_name])
+    # RAiDER needs strings for paths
+    side_effect = list(map(str, side_effect))
+    mocker.patch('RAiDER.processWM.prepareWeatherModel',
+                 side_effect=side_effect)
+    iargs_0 = [
+               '--file', str(out_path_0),
+               '--weather-model', weather_model_name,
+               '-interp', 'center_time'
+               ]
+    calcDelaysGUNW(iargs_0)
+
+    iargs_1 = [
+               '--file', str(out_path_1),
+               '--weather-model', weather_model_name,
+               '-interp', 'azimuth_time_grid'
+               ]
+    calcDelaysGUNW(iargs_1)
+
+    # Calls 6 times for azimuth time and 4 times for center time
+    assert RAiDER.processWM.prepareWeatherModel.call_count == 10
+    # Only calls for azimuth timing for reference and secondary
+    assert hyp3lib.get_orb.downloadSentinelOrbitFile.call_count == 2
 
     for ifg_type in ['reference', 'secondary']:
         for var in ['troposphereHydrostatic', 'troposphereWet']:
-            group = f'science/grids/corrections/external/troposphere/{model}/{ifg_type}'
-            with xr.open_dataset(out_0, group=group) as ds:
+            group = f'science/grids/corrections/external/troposphere/{weather_model_name}/{ifg_type}'
+            with xr.open_dataset(out_path_0, group=group) as ds:
                 da_0 = ds[var]
-            with xr.open_dataset(out_1, group=group) as ds:
+            with xr.open_dataset(out_path_1, group=group) as ds:
                 da_1 = ds[var]
             # diff * wavelength / (4 pi) transforms to meters; then x 1000 to mm
-            diff_mm = (da_1 - da_0).data * 0.055465761572122574 / (4 * np.pi) * 1_000
+            abs_diff_mm = np.abs((da_1 - da_0).data) * 0.055465761572122574 / (4 * np.pi) * 1_000
             # Differences in mm are bounded by 1
-            assert np.all(diff_mm < 1)
+            assert np.nanmax(abs_diff_mm) < 1
+
+
 @pytest.mark.parametrize('weather_model_name', ['GMAO', 'HRRR', 'HRES', 'ERA5', 'ERA5'])
 def test_check_weather_model_availability(test_gunw_path_factory, weather_model_name, mocker):
     # Should be True for all weather models

test/test_s1_time_grid.py

@@ -58,8 +58,8 @@ def test_get_slc_id():
 @pytest.mark.parametrize('ifg_type', ['reference', 'secondary'],)
 def test_s1_timing_array_wrt_slc_center_time(gunw_azimuth_test: Path,
                                              ifg_type: str,
-                                             orbit_dict_for_azimuth_test: dict,
-                                             slc_id_dict_for_azimuth_test: dict,
+                                             orbit_dict_for_azimuth_time_test: dict,
+                                             slc_id_dict_for_azimuth_time_test: dict,
                                              mocker):
     """Make sure the SLC start time is within reasonable amount of grid. The flow chart is:
 
@@ -88,9 +88,9 @@ def test_s1_timing_array_wrt_slc_center_time(gunw_azimuth_test: Path,
     # Azimuth time grid
     mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile',
                  # Hyp3 Lib returns 2 values
-                 return_value=(orbit_dict_for_azimuth_test[ifg_type], ''))
+                 return_value=(orbit_dict_for_azimuth_time_test[ifg_type], ''))
     mocker.patch('RAiDER.s1_azimuth_timing._asf_query',
-                 return_value=(slc_id_dict_for_azimuth_test[ifg_type], ''))
+                 return_value=(slc_id_dict_for_azimuth_time_test[ifg_type], ''))
     time_grid = get_s1_azimuth_time_grid(lon, lat, hgt, slc_start_time)
 
     abs_diff = np.abs(time_grid - np.datetime64(slc_start_time)) / np.timedelta64(1, 's')
@@ -106,8 +106,8 @@ def test_s1_timing_array_wrt_slc_center_time(gunw_azimuth_test: Path,
 @pytest.mark.parametrize('ifg_type', ['reference', 'secondary'])
 def test_s1_timing_array_wrt_variance(gunw_azimuth_test: Path,
                                       ifg_type: str,
-                                      orbit_dict_for_azimuth_test: dict,
-                                      slc_id_dict_for_azimuth_test: dict,
+                                      orbit_dict_for_azimuth_time_test: dict,
+                                      slc_id_dict_for_azimuth_time_test: dict,
                                       mocker):
     """Make sure along the hgt dimension of grid there is very small deviations
     """
@@ -132,9 +132,9 @@ def test_s1_timing_array_wrt_variance(gunw_azimuth_test: Path,
     # Azimuth time grid
     mocker.patch('hyp3lib.get_orb.downloadSentinelOrbitFile',
                  # Hyp3 Lib returns 2 values
-                 return_value=(orbit_dict_for_azimuth_test[ifg_type], ''))
+                 return_value=(orbit_dict_for_azimuth_time_test[ifg_type], ''))
     mocker.patch('RAiDER.s1_azimuth_timing._asf_query',
-                 return_value=(slc_id_dict_for_azimuth_test[ifg_type], ''))
+                 return_value=(slc_id_dict_for_azimuth_time_test[ifg_type], ''))
     X = get_s1_azimuth_time_grid(lon, lat, hgt, slc_start_time)
 
     Z = (X - X.min()) / np.timedelta64(1, 's')

tools/RAiDER/cli/raider.py

@@ -135,9 +135,9 @@ def drop_nans(d):
 def calcDelays(iargs=None):
     """ Parse command line arguments using argparse. """
     import RAiDER
+    import RAiDER.processWM
     from RAiDER.delay import tropo_delay
     from RAiDER.checkArgs import checkArgs
-    from RAiDER.processWM import prepareWeatherModel
     from RAiDER.utilFcns import writeDelays, get_nearest_wmtimes
     examples = 'Examples of use:' \
         '\n\t raider.py customTemplatefile.cfg' \
@@ -263,7 +263,10 @@ def calcDelays(iargs=None):
         wfiles = []
         for tt in times:
             try:
-                wfile = prepareWeatherModel(model, tt, aoi.bounds(), makePlots=params['verbose'])
+                wfile = RAiDER.processWM.prepareWeatherModel(model,
+                                                             tt,
+                                                             aoi.bounds(),
+                                                             makePlots=params['verbose'])
                 wfiles.append(wfile)
 
             # catch when requested datetime fails
@@ -541,6 +544,7 @@ def calcDelaysGUNW(iargs: list[str] = None):
 
     p.add_argument(
         '-interp', '--interpolate-time', default='center_time', type=str,
+        choices=['none', 'center_time', 'azimuth_time_grid'],
         help=('How to interpolate across model time steps. Possible options are: '
               '[\'none\', \'center_time\', \'azimuth_time_grid\'] '
               'None: means nearest model time; center_time: linearly across center time; '