Weather model-based troposphere delay LUTs (#103)

* Add RAiDER to the requirements file

* Exposed troposphere to the CSLC interface

* Compute troposphere delay using weather model

* Save computed troposphere range LUT

* Add Docker spec file

* Format delay_type logic to avoid repetition

* reduce line of codes

* fix specifile
only forge - no default or .condarc channels

* correctly force forge w/o defauluts

* Correct zenith delay formula

* Remove unnecessary flags that forced conda forge channel

---------

Co-authored-by: Liang Yu <liangjyu@gmail.com>

requirements.txt

@@ -8,6 +8,7 @@ pandas
 pyproj
 pysolid
 pytest
+raider
 ruamel.yaml
 scipy
 yamale

src/compass/defaults/s1_cslc_geo.yaml

@@ -17,6 +17,8 @@ runconfig:
       dynamic_ancillary_file_group:
           # Digital elevation model
           dem_file:
+          # Troposphere weather model file
+          weather_model_file:
 
       static_ancillary_file_group:
           # burst database sqlite file
@@ -57,6 +59,11 @@ runconfig:
               range_spacing: 120
               # LUT spacing in azimuth direction in seconds
               azimuth_spacing: 0.028
+              # Troposphere delay using weather model
+              troposphere:
+                 # Type of troposphere delay. Any of 'dry', 'wet' or 'wet_dry' for
+                 # the sum of wet and dry delays
+                 delay_type: wet_dry
 
           rdr2geo:
               # Enable/disable computation of topo layers

src/compass/s1_geocode_slc.py

@@ -81,9 +81,11 @@ def run(cfg: GeoRunConfig):
         if cfg.lut_params.enabled:
             rg_lut, az_lut = cumulative_correction_luts(burst,
                                                         dem_path=cfg.dem,
+                                                        scratch_path=scratch_path,
+                                                        weather_model_path=cfg.weather_model_file,
                                                         rg_step=cfg.lut_params.range_spacing,
                                                         az_step=cfg.lut_params.azimuth_spacing,
-                                                        scratch_path=scratch_path)
+                                                        delay_type=cfg.tropo_params.delay_type)
         else:
             rg_lut = isce3.core.LUT2d()
             az_lut = isce3.core.LUT2d()
@@ -186,7 +188,9 @@ def run(cfg: GeoRunConfig):
 
             cslc_group = geo_burst_h5.require_group(f'{root_path}/CSLC')
             metadata_to_h5group(cslc_group, burst, cfg)
-            corrections_to_h5group(cslc_group, burst, cfg, rg_lut, az_lut, scratch_path)
+            corrections_to_h5group(cslc_group, burst, cfg, rg_lut, az_lut, scratch_path,
+                                   weather_model_path=cfg.weather_model_file,
+                                   delay_type=cfg.tropo_params.delay_type)
 
     dt = str(timedelta(seconds=time.time() - t_start)).split(".")[0]
     info_channel.log(f"{module_name} burst successfully ran in {dt} (hr:min:sec)")

src/compass/schemas/s1_cslc_geo.yaml

@@ -17,6 +17,8 @@ runconfig:
         dynamic_ancillary_file_group:
            # Digital Elevation Model.
            dem_file: str(required=False)
+           # Troposphere weather model
+           weather_model_file: str(required=False)
 
         static_ancillary_file_group:
            # burst database sqlite file
@@ -86,6 +88,13 @@ lut_options:
   range_spacing: num(min=0, required=False)
   # LUT spacing in azimuth direction in seconds
   azimuth_spacing: num(min=0, required=False)
+  # Troposphere delay using weather model
+  troposphere: include('troposphere_options', required=False)
+
+troposphere_options:
+  # Type of troposphere delay. Any of 'dry', 'wet' or 'wet_dry' for
+  # the sum of wet and dry delays
+  delay_type: enum('dry', 'wet', 'wet_dry', required=False)
 
 rdr2geo_options:
    # Enable/disable computation of topo layers

src/compass/utils/geo_runconfig.py

@@ -10,6 +10,7 @@
 
 from compass.utils.geo_grid import (generate_geogrids_from_db,
                                     generate_geogrids, geogrid_as_dict)
+from compass.utils.helpers import check_file_path
 from compass.utils.runconfig import (
     create_output_paths,
     runconfig_to_bursts,
@@ -64,6 +65,14 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str) -> GeoRunConfig:
         geocoding_dict = groups_cfg['processing']['geocoding']
         check_geocode_dict(geocoding_dict)
 
+        # Check troposphere weather model file if not None. This
+        # troposphere correction is applied only if this file is not None
+        weather_model_path = groups_cfg['dynamic_ancillary_file_group'][
+            'weather_model_file'
+        ]
+        if weather_model_path is not None:
+            check_file_path(weather_model_path)
+
         # Convert runconfig dict to SimpleNamespace
         sns = wrap_namespace(groups_cfg)
 
@@ -92,6 +101,10 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str) -> GeoRunConfig:
         return cls(cfg['runconfig']['name'], sns, bursts, empty_ref_dict,
                    user_plus_default_yaml_str, output_paths, geogrids)
 
+    @property
+    def weather_model_file(self) -> str:
+        return self.groups.dynamic_ancillary_file_group.weather_model_file
+
     @property
     def geocoding_params(self) -> dict:
         return self.groups.processing.geocoding
@@ -104,6 +117,10 @@ def rdr2geo_params(self) -> dict:
     def lut_params(self) -> dict:
         return self.groups.processing.correction_luts
 
+    @property
+    def tropo_params(self) -> dict:
+        return self.groups.processing.correction_luts.troposphere
+
     def as_dict(self):
         ''' Convert self to dict for write to YAML/JSON
 

src/compass/utils/h5_helpers.py

@@ -550,7 +550,8 @@ def poly1d_to_h5(group, poly1d_name, poly1d):
 
 
 def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut,
-                           scratch_path):
+                           scratch_path, weather_model_path=None,
+                           delay_type='dry'):
     '''
     Write azimuth, slant range, and EAP (if needed) correction LUT2ds to HDF5
 
@@ -568,6 +569,14 @@ def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut,
         LUT2d along azimuth direction
     scratch_path: str
         Path to the scratch directory
+    weather_model_path: str
+        Path to troposphere weather model in NetCDF4 format.
+        This is the only format supported by RAiDER. If None,
+        no weather model-based troposphere correction is applied
+        (default: None).
+    delay_type: str
+        Type of troposphere delay. Any between 'dry', or 'wet', or
+        'wet_dry' for the sum of wet and dry troposphere delays.
     '''
 
     # If enabled, save the correction LUTs
@@ -609,6 +618,18 @@ def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut,
                  f'Solid Earth tides (range) {desc}',
                  {'units': 'meters'}),
         ]
+        if weather_model_path is not None:
+            if 'wet' in delay_type:
+                correction_items.append(Meta('wet_los_troposphere_delay',
+                                             ds.GetRasterBand(5).ReadAsArray(),
+                                             f'Wet LOS troposphere delay {desc}',
+                                             {'units': 'meters'}))
+            if 'dry' in delay_type:
+                correction_items.append(Meta('dry_los_troposphere_delay',
+                                             ds.GetRasterBand(6).ReadAsArray(),
+                                             f'Dry LOS troposphere delay {desc}',
+                                             {'units': 'meters'}))
+
         for meta_item in correction_items:
             add_dataset_and_attrs(correction_group, meta_item)
 

src/compass/utils/lut.py

@@ -12,11 +12,16 @@
 from compass.utils.geometry_utils import enu2los, en2az
 from compass.utils.helpers import open_raster
 from compass.utils.helpers import write_raster
+from RAiDER.delay import tropo_delay
+from RAiDER.llreader import RasterRDR
+from RAiDER.losreader import Zenith
 
 
 def cumulative_correction_luts(burst, dem_path,
+                               scratch_path=None,
+                               weather_model_path=None,
                                rg_step=200, az_step=0.25,
-                               scratch_path=None):
+                               delay_type='dry'):
     '''
     Sum correction LUTs and returns cumulative correction LUT in slant range
     and azimuth directions
@@ -27,12 +32,20 @@ def cumulative_correction_luts(burst, dem_path,
         Sentinel-1 A/B burst SLC object
     dem_path: str
         Path to the DEM file
+    scratch_path: str
+        Path to the scratch directory
+    weather_model_path: str
+        Path to the weather model file in NetCDF4 format.
+        This file has been preprocessed by RAiDER and it is
+        the only file format supported by the package. If None,
+        no troposphere correction is performed.
     rg_step: float
         LUT spacing along slant range direction
     az_step: float
         LUT spacing along azimuth direction
-    scratch_path: str
-        Path to the scratch directory
+    delay_type: str
+        Type of troposphere delay. Any between 'dry', or 'wet', or
+        'wet_dry' for the sum of wet and dry troposphere delays.
 
     Returns
     -------
@@ -44,17 +57,26 @@ def cumulative_correction_luts(burst, dem_path,
         and slant range
     '''
     # Get individual LUTs
-    geometrical_steer_doppler, bistatic_delay, az_fm_mismatch, [tide_rg, _]= \
+    geometrical_steer_doppler, bistatic_delay, az_fm_mismatch, [tide_rg, _], \
+        [wet_los_tropo, dry_los_tropo] = \
         compute_geocoding_correction_luts(burst,
                                           dem_path=dem_path,
+                                          scratch_path=scratch_path,
+                                          weather_model_path=weather_model_path,
                                           rg_step=rg_step,
                                           az_step=az_step,
-                                          scratch_path=scratch_path)
+                                          )
 
     # Convert to geometrical doppler from range time (seconds) to range (m)
     geometry_doppler = geometrical_steer_doppler.data * isce3.core.speed_of_light * 0.5
     rg_lut_data = geometry_doppler + tide_rg
 
+    # Add troposphere delay to range LUT
+    if 'wet' in delay_type:
+        rg_lut_data += wet_los_tropo
+    if 'dry' in delay_type:
+        rg_lut_data += dry_los_tropo
+
     # Invert signs to correct for convention
     # TO DO: add azimuth SET to LUT
     az_lut_data = -(bistatic_delay.data + az_fm_mismatch.data)
@@ -80,32 +102,46 @@ def cumulative_correction_luts(burst, dem_path,
     descr = ['geometrical doppler', 'bistatic delay', 'azimuth FM rate mismatch',
              'slant range Solid Earth tides']
 
+    if 'wet' in delay_type:
+        data_list.append(wet_los_tropo)
+        descr.append('wet LOS troposphere')
+    if 'dry' in delay_type:
+        data_list.append(dry_los_tropo)
+        descr.append('dry LOS troposphere')
+
     write_raster(f'{output_path}/corrections', data_list, descr)
 
     return rg_lut, az_lut
 
 
 def compute_geocoding_correction_luts(burst, dem_path,
+                                      scratch_path=None,
+                                      weather_model_path=None,
                                       rg_step=200, az_step=0.25,
-                                      scratch_path=None):
+                                      ):
     '''
     Compute slant range and azimuth LUTs corrections
     to be applied during burst geocoding
+
     Parameters
     ----------
     burst: Sentinel1BurstSlc
         S1-A/B burst object
     dem_path: str
         Path to the DEM required for azimuth FM rate mismatch.
-    xstep: int
-        LUT spacing along x/slant range in meters
-    ystep: int
-        LUT spacing along y/azimuth in seconds
-
     scratch_path: str
         Path to the scratch directory.
         If `None`, `burst.az_fm_rate_mismatch_mitigation()` will
         create temporary directory internally.
+    weather_model_path: str
+        Path to troposphere weather model in NetCDF4 format.
+        This is the only format supported by RAiDER. If None,
+        no weather model-based troposphere correction is applied
+        (default: None).
+    rg_step: int
+        LUT spacing along slant range in meters
+    az_step: int
+        LUT spacing along azimuth in seconds
 
     Returns
     -------
@@ -127,11 +163,17 @@ def compute_geocoding_correction_luts(burst, dem_path,
         This correction needs to be added to the SLC tagged azimuth time to
         get the corrected azimuth times.
 
-    [rg_set, az_set]: list, np.ndarray
+    [rg_set, az_set]: list[np.ndarray]
         List of numpy.ndarray containing SET in slant range and azimuth directions
         in meters. These corrections need to be added to the slC tagged azimuth
         and slant range times.
+
+    [wet_los_tropo, dry_los_tropo]: list[np.ndarray]
+        List of numpy.ndarray containing the LOS wet and dry troposphere delays
+        computed from the file specified under 'weather_model_path'. These delays
+        need to be added to the slant range correction LUT2D.
     '''
+
     # Get DEM raster
     dem_raster = isce3.io.Raster(dem_path)
     epsg = dem_raster.get_epsg()
@@ -183,8 +225,30 @@ def compute_geocoding_correction_luts(burst, dem_path,
     rg_set = resize(rg_set_temp, out_shape, **kwargs)
     az_set = resize(az_set_temp, out_shape, **kwargs)
 
+    # Compute wet and dry troposphere delays using RAiDER
+    wet_los_tropo, dry_los_tropo = [np.zeros(out_shape) for _ in range(2)]
+
+    if weather_model_path is not None:
+        # Instantiate an "aoi" object to read lat/lon/height files
+        aoi = RasterRDR(rdr2geo_raster_paths[1], rdr2geo_raster_paths[0],
+                        rdr2geo_raster_paths[2])
+
+        # Instantiate the Zenith object. Note RAiDER LOS object requires
+        # the orbit files.
+        los = Zenith()
+
+        # Compute the troposphere delay along the Zenith
+        zen_wet, zen_dry = tropo_delay(burst.sensing_start,
+                                       weather_model_path,
+                                       aoi, los)
+
+        # RaiDER delay is one-way only. Get the LOS delay my multiplying
+        # by the incidence angle
+        wet_los_tropo = 2.0 * zen_wet / np.cos(np.deg2rad(inc_angle))
+        dry_los_tropo = 2.0 * zen_dry / np.cos(np.deg2rad(inc_angle))
+
     return geometrical_steering_doppler, bistatic_delay, az_fm_mismatch, [
-        rg_set, az_set]
+        rg_set, az_set], [wet_los_tropo, dry_los_tropo]
 
 
 def solid_earth_tides(burst, lat_radar_grid, lon_radar_grid, inc_angle,