From 65100af3f6b04fa13596c0dbfa6634908eef8871 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 14:31:31 -0700 Subject: [PATCH 01/16] Expose boolean flags to turn on/off individual corrections --- src/compass/defaults/s1_cslc_geo.yaml | 22 +++++++++++++++------- src/compass/schemas/s1_cslc_geo.yaml | 22 ++++++++++++++++++---- 2 files changed, 33 insertions(+), 11 deletions(-) diff --git a/src/compass/defaults/s1_cslc_geo.yaml b/src/compass/defaults/s1_cslc_geo.yaml index 3b8fb0c5..b3aa5087 100644 --- a/src/compass/defaults/s1_cslc_geo.yaml +++ b/src/compass/defaults/s1_cslc_geo.yaml @@ -55,18 +55,26 @@ runconfig: numiter: 25 correction_luts: - # Boolean flag to activate/deactivate model-based - # corrections while geocoding the burst - enabled: True # LUT spacing in range direction in meters range_spacing: 120 # LUT spacing in azimuth direction in seconds azimuth_spacing: 0.028 + # Enable/disable geometry steering doppler correction + geometry_steering_doppler: True + # Enable/disable bistatic delay correction + bistatic_delay: True + # Enable/disable azimuth FM rate mismatch correction + azimuth_fm_rate: True + # Enable/disable Solid Earth tides correction + solid_earth_tide: True + # Enable/disable ionosphere TEC correction + tec_ionosphere: True + # Enable/disable static troposphere correction + static_troposphere: True # 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 + weather_model_troposphere: + enabled: False + delay_type: wet_dry rdr2geo: # Enable/disable computation of topo layers diff --git a/src/compass/schemas/s1_cslc_geo.yaml b/src/compass/schemas/s1_cslc_geo.yaml index fc5d8603..12430b5b 100644 --- a/src/compass/schemas/s1_cslc_geo.yaml +++ b/src/compass/schemas/s1_cslc_geo.yaml @@ -85,17 +85,31 @@ geo2rdr_options: lines_per_block: int(min=1, required=False) lut_options: - # Boolean flag to activate/deactivate model-based - # corrections while geocoding the burst - enabled: bool(required=False) # LUT spacing in range direction in meters range_spacing: num(min=0, required=False) # LUT spacing in azimuth direction in seconds azimuth_spacing: num(min=0, required=False) + # Enable/disable geometry steering doppler correction + geometry_steering_doppler: include('enable_options', required=False) + # Enable/disable bistatic delay correction + bistatic_delay: include('enable_options', required=False) + # Enable/disable azimuth FM rate mismatch correction + azimuth_fm_rate: include('enable_options', required=False) + # Enable/disable Solid Earth tides correction + solid_earth_tide: include('enable_options', required=False) + # Enable/disable ionosphere TEC correction + tec_ionosphere: include('enable_options', required=False) + # Enable/disable static troposphere correction + static_troposphere: include('enable_options', required=False) # Troposphere delay using weather model - troposphere: include('troposphere_options', required=False) + weather_model_troposphere: include('troposphere_options', required=False) + +enable_options: + # Boolean flag to enable a LUT correction + enabled: bool(required=False) troposphere_options: + enabled: bool(required=False) # 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) From b333166f022497a891daf48e0b14646e9be27dc5 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 18:10:55 -0700 Subject: [PATCH 02/16] Correct typos in schema entry --- src/compass/defaults/s1_cslc_geo.yaml | 4 ++-- src/compass/schemas/s1_cslc_geo.yaml | 18 +++++++----------- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/src/compass/defaults/s1_cslc_geo.yaml b/src/compass/defaults/s1_cslc_geo.yaml index b3aa5087..2644c276 100644 --- a/src/compass/defaults/s1_cslc_geo.yaml +++ b/src/compass/defaults/s1_cslc_geo.yaml @@ -66,9 +66,9 @@ runconfig: # Enable/disable azimuth FM rate mismatch correction azimuth_fm_rate: True # Enable/disable Solid Earth tides correction - solid_earth_tide: True + solid_earth_tides: True # Enable/disable ionosphere TEC correction - tec_ionosphere: True + ionosphere_tec: True # Enable/disable static troposphere correction static_troposphere: True # Troposphere delay using weather model diff --git a/src/compass/schemas/s1_cslc_geo.yaml b/src/compass/schemas/s1_cslc_geo.yaml index 12430b5b..9f41cf3f 100644 --- a/src/compass/schemas/s1_cslc_geo.yaml +++ b/src/compass/schemas/s1_cslc_geo.yaml @@ -90,23 +90,19 @@ lut_options: # LUT spacing in azimuth direction in seconds azimuth_spacing: num(min=0, required=False) # Enable/disable geometry steering doppler correction - geometry_steering_doppler: include('enable_options', required=False) + geometry_steering_doppler: bool(required=False) # Enable/disable bistatic delay correction - bistatic_delay: include('enable_options', required=False) + bistatic_delay: bool(required=False) # Enable/disable azimuth FM rate mismatch correction - azimuth_fm_rate: include('enable_options', required=False) + azimuth_fm_rate: bool(required=False) # Enable/disable Solid Earth tides correction - solid_earth_tide: include('enable_options', required=False) + solid_earth_tides: bool(required=False) # Enable/disable ionosphere TEC correction - tec_ionosphere: include('enable_options', required=False) + ionosphere_tec: bool(required=False) # Enable/disable static troposphere correction - static_troposphere: include('enable_options', required=False) + static_troposphere: bool(required=False) # Troposphere delay using weather model - weather_model_troposphere: include('troposphere_options', required=False) - -enable_options: - # Boolean flag to enable a LUT correction - enabled: bool(required=False) + weather_model_troposphere: bool(required=False) troposphere_options: enabled: bool(required=False) From 930cb03272eb0c25b4665b7a3623533f20424180 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 18:17:14 -0700 Subject: [PATCH 03/16] Selectively compute LUT corrections --- src/compass/utils/lut.py | 340 ++++++++++++++------------------------- 1 file changed, 119 insertions(+), 221 deletions(-) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index aeed9c95..294390be 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -18,236 +18,120 @@ from RAiDER.losreader import Zenith -def cumulative_correction_luts(burst, dem_path, tec_path, - scratch_path=None, - weather_model_path=None, - rg_step=200, az_step=0.25, - delay_type='dry'): - ''' - Sum correction LUTs and returns cumulative correction LUT in slant range - and azimuth directions - - Parameters - ---------- - burst: Sentinel1BurstSlc - Sentinel-1 A/B burst SLC object - dem_path: str - Path to the DEM file - tec_path: str - Path to the TEC file in IONEX format - 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 - 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 - ------- - rg_lut: isce3.core.LUT2d - Sum of slant range correction LUTs in meters as a function of azimuth - time and slant range - az_lut: isce3.core.LUT2d - Sum of azimuth correction LUTs in seconds as a function of azimuth time - and slant range - ''' - # Get individual LUTs - geometrical_steer_doppler, bistatic_delay, az_fm_mismatch, [tide_rg, _], \ - los_ionosphere, [wet_los_tropo, dry_los_tropo], los_static_tropo = \ - compute_geocoding_correction_luts(burst, - dem_path=dem_path, - tec_path=tec_path, - scratch_path=scratch_path, - weather_model_path=weather_model_path, - rg_step=rg_step, - az_step=az_step) - - # 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 + los_ionosphere + los_static_tropo - - # 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) - - rg_lut = isce3.core.LUT2d(bistatic_delay.x_start, - bistatic_delay.y_start, - bistatic_delay.x_spacing, - bistatic_delay.y_spacing, - rg_lut_data) - az_lut = isce3.core.LUT2d(bistatic_delay.x_start, - bistatic_delay.y_start, - bistatic_delay.x_spacing, - bistatic_delay.y_spacing, - az_lut_data) - - # Save corrections on disk. In this way, we should avoid running - # the corrections again when allocating data inside the HDF5 product - # Create a directory in the scratch path to save corrections - output_path = f'{scratch_path}/corrections' - os.makedirs(output_path, exist_ok=True) - data_list = [geometry_doppler, bistatic_delay.data, az_fm_mismatch.data, - tide_rg, los_ionosphere] - descr = ['slant range geometrical doppler', 'azimuth bistatic delay', 'azimuth FM rate mismatch', - 'slant range Solid Earth tides', 'line-of-sight ionospheric delay'] - - if weather_model_path is not None: - 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, tec_path, - scratch_path=None, - weather_model_path=None, - rg_step=200, az_step=0.25,): - ''' - 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. - tec_path: str - Path to the TEC file for ionosphere correction - 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 - ------- - geometrical_steering_doppler: isce3.core.LUT2d: - LUT2D object of total doppler (geometrical doppler + steering doppler) - in seconds as the function of the azimuth time and slant range. - This correction needs to be added to the SLC tagged range time to - get the corrected range times. - - bistatic_delay: isce3.core.LUT2d: - LUT2D object of bistatic delay correction in seconds as a function - of the azimuth time and slant range. - This correction needs to be added to the SLC tagged azimuth time to - get the corrected azimuth times. - - az_fm_mismatch: isce3.core.LUT2d: - LUT2D object of azimuth FM rate mismatch mitigation, - in seconds as the function of the azimuth time and slant range. - 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] - 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. - - ionosphere: np.ndarray - numpy.ndarray for ionosphere delay in line-of-sight direction in meters. - This correction needs to be added to the SLC tagged range time to - get the corrected 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 +def correction_luts(burst, lut_par, dem_path, tec_path, + scratch_path=None, + weather_model_path=None): + # Dem info dem_raster = isce3.io.Raster(dem_path) epsg = dem_raster.get_epsg() proj = isce3.core.make_projection(epsg) ellipsoid = proj.ellipsoid - # Create directory to store SET temp results + # Get LUT spacing + rg_step = lut_par.range_spacing + az_step = lut_par.azimuth_spacing + + # Create directory to temporary results output_path = f'{scratch_path}/corrections' os.makedirs(output_path, exist_ok=True) - # Compute Geometrical Steering Doppler - geometrical_steering_doppler = \ - burst.doppler_induced_range_shift(range_step=rg_step, az_step=az_step) - - # Compute bistatic delay - bistatic_delay = burst.bistatic_delay(range_step=rg_step, az_step=az_step) - - # Run rdr2geo to obtain the required layers - # return contents: lon_path, lat_path, height_path, inc_path, head_path - rdr2geo_raster_paths = compute_rdr2geo_rasters(burst, dem_raster, - output_path, rg_step, - az_step) - - # Open rdr2geo layers - lon, lat, height, inc_angle, head_angle = \ - [open_raster(raster_path) for raster_path in rdr2geo_raster_paths] - - # Compute azimuth FM-rate mismatch - az_fm_mismatch = burst.az_fm_rate_mismatch_from_llh(lat, lon, height, + # If any of the following corrections is enabled + # generate rdr2geo layers + rdr2geo_enabled = lut_par.azimuth_fm_rate or \ + lut_par.solid_earth_tides or \ + lut_par.ionosphere_tec or \ + lut_par.static_troposphere or \ + lut_par.weather_model_troposphere + if rdr2geo_enabled: + # return contents: lon_path, lat_path, height_path, inc_path, head_path + rdr2geo_raster_paths = compute_rdr2geo_rasters(burst, dem_raster, + output_path, rg_step, + az_step) + # Open rdr2geo layers + lon, lat, height, inc_angle, head_angle = \ + [open_raster(raster_path) for raster_path in rdr2geo_raster_paths] + + # Get the shape of the correction LUT and create empty numpy array + lut = burst.bistatic_delay(range_step=rg_step, + az_step=az_step) + lut_shape = lut.data.shape + rg_data = np.zeros(lut_shape, dtype=np.float32) + az_data = np.zeros(lut_shape, dtype=np.float32) + + # Initialize data list and description to save corrections + data_dict = { + 'doppler': [rg_data, 'Slant range geometry and steering doppler'], + 'bistatic_delay': [rg_data, 'Bistatic delay'], + 'azimuth_fm_rate': [rg_data, 'Azimuth FM rate mismatch'], + 'rg_set': [rg_data, 'Slant range Solid Earth Tides'], + 'az_set': [rg_data, 'Azimuth Solid Earth Tides'], + 'static_tropo': [rg_data, 'Static troposphere delay'], + 'tec_iono': [rg_data, 'Slant range ionosphere delay'], + 'dry_tropo': [rg_data, 'Dry troposphere delay from weather model'], + 'wet_tropo': [rg_data, 'Wet troposphere delay from weather model']} + + # Check which corrections are requested and accumulate corresponding data + # Geometrical and steering doppler + if lut_par.geometry_steering_doppler: + doppler = burst.doppler_induced_range_shift(range_step=rg_step, + az_step=az_step) + doppler_meter = doppler.data * isce3.core.speed_of_light * 0.5 + rg_data += doppler_meter + data_dict['doppler'][0] = doppler_meter + + # Bistatic delay + if lut_par.bistatic_delay: + bistatic_delay = burst.bistatic_delay(range_step=rg_step, + az_step=az_step).data + az_data -= bistatic_delay + data_dict['bistatic_delay'][0] = -bistatic_delay + + # Azimuth FM-rate mismatch + if lut_par.azimuth_fm_rate: + az_fm_rate = burst.az_fm_rate_mismatch_from_llh(lat, lon, height, ellipsoid, burst.as_isce3_radargrid( az_step=az_step, - rg_step=rg_step) - ) - - # compute Solid Earth Tides using pySolid. Decimate the rdr2geo layers. - # compute decimation factor assuming a 5 km spacing along slant range - dec_factor = int(np.round(5000.0 / rg_step)) - dec_slice = np.s_[::dec_factor] - rg_set_temp, az_set_temp = solid_earth_tides(burst, lat[dec_slice], - lon[dec_slice], - inc_angle[dec_slice], - head_angle[dec_slice]) - - # Resize SET to the size of the correction grid - out_shape = bistatic_delay.data.shape - kwargs = dict(order=1, mode='edge', anti_aliasing=True, - preserve_range=True) - rg_set = resize(rg_set_temp, out_shape, **kwargs) - az_set = resize(az_set_temp, out_shape, **kwargs) - - # Compute ionosphere delay - los_ionosphere = ionosphere_delay(burst.sensing_mid, - burst.wavelength, - tec_path, lon, lat, inc_angle) - - # Compute wet and dry troposphere delays using RAiDER - wet_los_tropo, dry_los_tropo, los_static_tropo =\ - [np.zeros(out_shape) for _ in range(3)] - - if weather_model_path is None: - # Compute static troposphere correction + rg_step=rg_step)).data + az_data -= az_fm_rate + data_dict['azimuth_fm_rate'][0] = -az_fm_rate + + # Solid Earth tides + if lut_par.solid_earth_tides: + dec_factor = int(np.round(5000.0 / rg_step)) + dec_slice = np.s_[::dec_factor] + rg_set_temp, az_set_temp = solid_earth_tides(burst, lat[dec_slice], + lon[dec_slice], + inc_angle[dec_slice], + head_angle[dec_slice]) + + # Resize SET to the size of the correction grid + kwargs = dict(order=1, mode='edge', anti_aliasing=True, + preserve_range=True) + rg_set = resize(rg_set_temp, lut_shape, **kwargs) + az_set = resize(az_set_temp, lut_shape, **kwargs) + rg_data += rg_set + az_data += az_set + data_dict['rg_set'][0] = rg_set + data_dict['az_set'][0] = az_set + + # Static troposphere + if lut_par.static_troposphere: los_static_tropo = compute_static_troposphere_delay(inc_angle, height) - - else: + rg_data += los_static_tropo + data_dict['static_tropo'][0] = los_static_tropo + + # Ionosphere TEC correction + if lut_par.ionosphere_tec: + los_iono = ionosphere_delay(burst.sensing_mid, + burst.wavelength, + tec_path, lon, lat, inc_angle) + rg_data += los_iono + data_dict['tec_iono'][0] = los_iono + + # Weather model troposphere correction + if lut_par.weather_model_troposphere.enabled: + delay_type = lut_par.weather_model_troposphere.delay_type # Instantiate an "aoi" object to read lat/lon/height files aoi = RasterRDR(rdr2geo_raster_paths[1], rdr2geo_raster_paths[0], rdr2geo_raster_paths[2]) @@ -260,15 +144,29 @@ def compute_geocoding_correction_luts(burst, dem_path, tec_path, 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], los_ionosphere, [wet_los_tropo, dry_los_tropo], los_static_tropo + if 'wet' in delay_type: + rg_data += wet_los_tropo + data_dict['wet_tropo'][0] = wet_los_tropo + if 'dry' in delay_type: + rg_data += dry_los_tropo + data_dict['dry_tropo'][0] = dry_los_tropo + + # Create the range and azimuth LUT2d + rg_lut = isce3.core.LUT2d(lut.x_start, lut.y_start, + lut.x_spacing, lut.y_spacing, + rg_data) + az_lut = isce3.core.LUT2d(lut.x_start, lut.y_start, + lut.x_spacing, lut.y_spacing, + az_data) + # TO DO: save corrections + + return rg_lut, az_lut def solid_earth_tides(burst, lat_radar_grid, lon_radar_grid, inc_angle, From a7f797fa5a0b474f31400a746c83c5aefcfb5c8d Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 22:10:47 -0700 Subject: [PATCH 04/16] Save corrections --- src/compass/utils/lut.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 294390be..cc4baa57 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -12,7 +12,6 @@ from compass.utils.geometry_utils import enu2los, en2az from compass.utils.iono import ionosphere_delay 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 @@ -164,7 +163,17 @@ def correction_luts(burst, lut_par, dem_path, tec_path, lut.x_spacing, lut.y_spacing, az_data) - # TO DO: save corrections + # Save corrections + driver = gdal.GetDriverByName('ENVI') + out_ds = driver.Create(f'{output_path}/corrections', + lut_shape[1], lut_shape[0], len(data_dict), + gdal.GDT_Float32) + band = 0 + for key in data_dict.keys(): + band += 1 + raster_band = out_ds.GetRasterBand(band) + raster_band.SetDescription(data_dict[key][1]) + raster_band.WriteArray(data_dict[key][0]) return rg_lut, az_lut From fcad310745d59f03324d55e4b3dc17dff8b2aa02 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 22:11:10 -0700 Subject: [PATCH 05/16] Remove write raster function --- src/compass/utils/helpers.py | 52 ------------------------------------ 1 file changed, 52 deletions(-) diff --git a/src/compass/utils/helpers.py b/src/compass/utils/helpers.py index 777eb463..3a45bc10 100644 --- a/src/compass/utils/helpers.py +++ b/src/compass/utils/helpers.py @@ -345,58 +345,6 @@ def open_raster(filename, band=1): return raster -def write_raster(filename, data_list, descriptions, - data_type=gdal.GDT_Float32, data_format='ENVI'): - ''' - Write a multiband GDAL-friendly raster to disk. - Each dataset allocated in the output file contains - a description of the dataset allocated for that band - - Parameters - ---------- - filename: str - File path where to store output dataset - data_list: list[np.ndarray] - List of numpy.ndarray to allocate for each - raster band. All datasets within the list - are assumed to have the same shape - descriptions: list[str] - List of strings containing a description - for the bands to allocate - data_type: gdal.dtype - GDAL dataset type - format: gdal.Format - Format for GDAL output file - ''' - - error_channel = journal.error('helpers.write_raster') - - # Check number of datasets match number of descriptions - if len(data_list) != len(descriptions): - err_str = f'Number of datasets to write does not match' \ - f'the number of descriptions ' \ - f'{len(data_list)} != {len(descriptions)}' - error_channel.log(err_str) - raise ValueError(err_str) - - # Get the shape of a dataset within the list. All the datasets - # are assumed to have the same shape - length, width = data_list[0].shape - nbands = len(data_list) - - driver = gdal.GetDriverByName(data_format) - out_ds = driver.Create(filename, width, length, nbands, data_type) - - band = 0 - for data, description in zip(data_list, descriptions): - band += 1 - raster_band = out_ds.GetRasterBand(band) - raster_band.SetDescription(description) - raster_band.WriteArray(data) - - out_ds.FlushCache() - - def bursts_grouping_generator(bursts): ''' Dict to group bursts with the same burst ID but different polarizations From 494f9ff172cf3f25a9e0572685b3d7987f584b5e Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 22:41:05 -0700 Subject: [PATCH 06/16] Save correction LUTs --- src/compass/s1_geocode_slc.py | 43 +++++++++++++------------------ src/compass/utils/h5_helpers.py | 45 ++++++++++----------------------- 2 files changed, 31 insertions(+), 57 deletions(-) diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index e3775ea4..0bd89876 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -22,7 +22,7 @@ init_geocoded_dataset, metadata_to_h5group) from compass.utils.helpers import bursts_grouping_generator, get_module_name -from compass.utils.lut import cumulative_correction_luts +from compass.utils.lut import correction_luts from compass.utils.yaml_argparse import YamlArgparse @@ -70,20 +70,12 @@ def run(cfg: GeoRunConfig): # Create scratch as needed scratch_path = out_paths.scratch_directory - - # If enabled, get range and azimuth LUTs - if cfg.lut_params.enabled: - rg_lut, az_lut = \ - cumulative_correction_luts(burst, dem_path=cfg.dem, - tec_path=cfg.tec_file, - 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, - delay_type=cfg.tropo_params.delay_type) - else: - rg_lut = isce3.core.LUT2d() - az_lut = isce3.core.LUT2d() + # Compute correction LUTs + rg_lut, az_lut = correction_luts(burst, cfg.lut_params, + dem_path=cfg.dem, + tec_path=cfg.tec_file, + scratch_path=scratch_path, + weather_model_path=cfg.weather_model_file) radar_grid = burst.as_isce3_radargrid() native_doppler = burst.doppler.lut2d @@ -181,11 +173,10 @@ def run(cfg: GeoRunConfig): cslc_group = geo_burst_h5.require_group(f'{root_path}/CSLC') metadata_to_h5group(cslc_group, burst, cfg) - if cfg.lut_params.enabled: - 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) + + # Save corrections + corrections_to_h5group(cslc_group, burst, rg_lut, az_lut, + scratch_path) # If needed, make browse image and compute CSLC raster stats browse_params = cfg.browse_image_params @@ -199,12 +190,12 @@ def run(cfg: GeoRunConfig): # If needed, perform QA and write results to JSON if cfg.quality_assurance_params.perform_qa: cslc_qa = QualityAssuranceCSLC() - if cfg.lut_params.enabled: - # apply tropo corrections if weather file provided - apply_tropo_corrections = cfg.weather_model_file is not None - cslc_qa.compute_correction_stats( - geo_burst_h5, apply_tropo_corrections, - cfg.tropo_params.delay_type) + # if cfg.lut_params.enabled: + # # apply tropo corrections if weather file provided + # apply_tropo_corrections = cfg.weather_model_file is not None + # cslc_qa.compute_correction_stats( + # geo_burst_h5, apply_tropo_corrections, + # cfg.tropo_params.delay_type) cslc_qa.compute_CSLC_raster_stats(geo_burst_h5, bursts) cslc_qa.populate_rfi_dict(geo_burst_h5) cslc_qa.valid_pixel_percentages(geo_burst_h5) diff --git a/src/compass/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py index bbe48d47..0dea884f 100644 --- a/src/compass/utils/h5_helpers.py +++ b/src/compass/utils/h5_helpers.py @@ -580,9 +580,8 @@ def poly1d_to_h5(group, poly1d_name, poly1d): poly1d_to_h5(burst_meta_group, 'doppler', burst.doppler.poly1d) -def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut, - scratch_path, weather_model_path=None, - delay_type='dry'): +def corrections_to_h5group(parent_group, burst, rg_lut, + az_lut, scratch_path): ''' Write azimuth, slant range, and EAP (if needed) correction LUT2ds to HDF5 @@ -592,28 +591,14 @@ def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut, HDF5 group where correction data will be written to burst: Sentinel1BurstSlc Burst containing corrections - cfg: types.SimpleNamespace - SimpleNamespace containing run configuration rg_lut: isce3.core.LUT2d() LUT2d along slant direction az_lut: isce3.core.LUT2d() 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 - if not cfg.lut_params.enabled: - return - # Open GDAL dataset to fetch corrections ds = gdal.Open(f'{scratch_path}/corrections/corrections', gdal.GA_ReadOnly) @@ -638,32 +623,30 @@ def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut, {'units': 'seconds'}), Meta('zero_doppler_time_spacing',rg_lut.y_spacing, 'spacing of azimuth time of LUT data', {'units': 'seconds'}), - Meta('bistatic_delay', ds.GetRasterBand(2).ReadAsArray(), - f'bistatic delay (azimuth) {desc}', {'units': 'seconds'}), Meta('geometry_steering_doppler', ds.GetRasterBand(1).ReadAsArray(), f'geometry steering doppler (range) {desc}', {'units': 'meters'}), + Meta('bistatic_delay', ds.GetRasterBand(2).ReadAsArray(), + f'bistatic delay (azimuth) {desc}', {'units': 'seconds'}), Meta('azimuth_fm_rate_mismatch', ds.GetRasterBand(3).ReadAsArray(), f'azimuth FM rate mismatch mitigation (azimuth) {desc}', {'units': 'seconds'}), Meta('los_solid_earth_tides', ds.GetRasterBand(4).ReadAsArray(), f'Solid Earth tides (range) {desc}', {'units': 'meters'}), - Meta('los_ionospheric_delay', ds.GetRasterBand(5).ReadAsArray(), + Meta('azimuth_solid_earth_tides', ds.GetRasterBand(5).ReadAsArray(), + f'Solid Earth tides (range) {desc}', + {'units': 'seconds'}), + Meta('los_ionospheric_delay', ds.GetRasterBand(7).ReadAsArray(), f'Ionospheric delay (range) {desc}', {'units': 'meters'}), + Meta('wet_los_troposphere_delay', ds.GetRasterBand(8).ReadAsArray(), + f'Wet LOS troposphere delay {desc}', + {'units': 'meters'}), + Meta('dry_los_troposphere_delay', ds.GetRasterBand(9).ReadAsArray(), + f'Dry LOS troposphere delay {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(6).ReadAsArray(), - f'Wet LOS troposphere delay {desc}', - {'units': 'meters'})) - if 'dry' in delay_type: - correction_items.append(Meta('dry_los_troposphere_delay', - ds.GetRasterBand(7).ReadAsArray(), - f'Dry LOS troposphere delay {desc}', - {'units': 'meters'})) for meta_item in correction_items: add_dataset_and_attrs(correction_group, meta_item) From b633ed501e9eff1e1f63e3ed9d8489b4630fdaee Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 22:41:27 -0700 Subject: [PATCH 07/16] Doc string correction function --- src/compass/utils/lut.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index cc4baa57..04251550 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -20,6 +20,32 @@ def correction_luts(burst, lut_par, dem_path, tec_path, scratch_path=None, weather_model_path=None): + ''' + Compute correction look-up tables (LUTs) + + Parameters + ---------- + burst: Sentinel1BurstSlc + S1-A/B burst SLC object + lut_par: dict + Dictionary with LUT parameters + dem_path: str + File path to DEM + tec_path: str + File path to ionosphere TEC file + scratch_path: str + File path to scratch directory (default: None) + weather_model_path: str + File path to weather model file (default: None) + + Returns + ------- + rg_lut: isce3.core.LUT2d + Cumulative LUT in slant range direction (meters) + qz_lut: isce3.core.LUT2d + Cumulative LUT in azimuth direction (seconds) + ''' + # Dem info dem_raster = isce3.io.Raster(dem_path) epsg = dem_raster.get_epsg() From 3f1a20157667ef6659af6444a43c6cc70b36ba20 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Fri, 31 Mar 2023 23:13:47 -0700 Subject: [PATCH 08/16] Address codacy issue --- src/compass/utils/lut.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 04251550..1ebdf8d4 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -45,7 +45,6 @@ def correction_luts(burst, lut_par, dem_path, tec_path, qz_lut: isce3.core.LUT2d Cumulative LUT in azimuth direction (seconds) ''' - # Dem info dem_raster = isce3.io.Raster(dem_path) epsg = dem_raster.get_epsg() From 8b01cf4a0a03f204724ae351f780d55884b6e0a5 Mon Sep 17 00:00:00 2001 From: vbrancat Date: Tue, 4 Apr 2023 14:13:12 -0700 Subject: [PATCH 09/16] Fix decimation factor along columns --- src/compass/utils/lut.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 1ebdf8d4..7e6271d1 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -124,10 +124,10 @@ def correction_luts(burst, lut_par, dem_path, tec_path, if lut_par.solid_earth_tides: dec_factor = int(np.round(5000.0 / rg_step)) dec_slice = np.s_[::dec_factor] - rg_set_temp, az_set_temp = solid_earth_tides(burst, lat[dec_slice], - lon[dec_slice], - inc_angle[dec_slice], - head_angle[dec_slice]) + rg_set_temp, az_set_temp = solid_earth_tides(burst, lat[dec_slice, dec_slice], + lon[dec_slice, dec_slice], + inc_angle[dec_slice, dec_slice], + head_angle[dec_slice, dec_slice]) # Resize SET to the size of the correction grid kwargs = dict(order=1, mode='edge', anti_aliasing=True, From 271526b2c5c34c6de4d1f2506bc9c92bd04a0fb6 Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Mon, 15 May 2023 15:43:24 -0700 Subject: [PATCH 10/16] tweak loops for list/dict comprehension typo fixes --- src/compass/schemas/s1_cslc_geo.yaml | 2 +- src/compass/utils/lut.py | 34 ++++++++++++++-------------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/compass/schemas/s1_cslc_geo.yaml b/src/compass/schemas/s1_cslc_geo.yaml index 9f41cf3f..f373dbe4 100644 --- a/src/compass/schemas/s1_cslc_geo.yaml +++ b/src/compass/schemas/s1_cslc_geo.yaml @@ -102,7 +102,7 @@ lut_options: # Enable/disable static troposphere correction static_troposphere: bool(required=False) # Troposphere delay using weather model - weather_model_troposphere: bool(required=False) + weather_model_troposphere: include('troposphere_options', required=False) troposphere_options: enabled: bool(required=False) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 7e6271d1..00ab9acb 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -42,7 +42,7 @@ def correction_luts(burst, lut_par, dem_path, tec_path, ------- rg_lut: isce3.core.LUT2d Cumulative LUT in slant range direction (meters) - qz_lut: isce3.core.LUT2d + az_lut: isce3.core.LUT2d Cumulative LUT in azimuth direction (seconds) ''' # Dem info @@ -83,16 +83,18 @@ def correction_luts(burst, lut_par, dem_path, tec_path, az_data = np.zeros(lut_shape, dtype=np.float32) # Initialize data list and description to save corrections - data_dict = { - 'doppler': [rg_data, 'Slant range geometry and steering doppler'], - 'bistatic_delay': [rg_data, 'Bistatic delay'], - 'azimuth_fm_rate': [rg_data, 'Azimuth FM rate mismatch'], - 'rg_set': [rg_data, 'Slant range Solid Earth Tides'], - 'az_set': [rg_data, 'Azimuth Solid Earth Tides'], - 'static_tropo': [rg_data, 'Static troposphere delay'], - 'tec_iono': [rg_data, 'Slant range ionosphere delay'], - 'dry_tropo': [rg_data, 'Dry troposphere delay from weather model'], - 'wet_tropo': [rg_data, 'Wet troposphere delay from weather model']} + data_dict_key_dscrs = ( + ['doppler', 'Slant range geometry and steering doppler'], + ['bistatic_delay', 'Bistatic delay'], + ['azimuth_fm_rate', 'Azimuth FM rate mismatch'], + ['rg_set', 'Slant range Solid Earth Tides'], + ['az_set', 'Azimuth Solid Earth Tides'], + ['static_tropo', 'Static troposphere delay'], + ['tec_iono', 'Slant range ionosphere delay'], + ['dry_tropo', 'Dry troposphere delay from weather model'], + ['wet_tropo', 'Wet troposphere delay from weather model']) + data_dict = {key: (np.zeros_like(rg_data), dscr) + for (key, dscr) in data_dict_key_dscrs} # Check which corrections are requested and accumulate corresponding data # Geometrical and steering doppler @@ -193,12 +195,10 @@ def correction_luts(burst, lut_par, dem_path, tec_path, out_ds = driver.Create(f'{output_path}/corrections', lut_shape[1], lut_shape[0], len(data_dict), gdal.GDT_Float32) - band = 0 - for key in data_dict.keys(): - band += 1 - raster_band = out_ds.GetRasterBand(band) - raster_band.SetDescription(data_dict[key][1]) - raster_band.WriteArray(data_dict[key][0]) + for band, (key, (data, desc)) in enumerate(data_dict.items()): + raster_band = out_ds.GetRasterBand(band + 1) + raster_band.SetDescription(desc) + raster_band.WriteArray(data) return rg_lut, az_lut From bbd83065195fe509e7bb38c886072c70236b519a Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Fri, 19 May 2023 15:32:27 -0700 Subject: [PATCH 11/16] write corrections to HDF5 when computed and flags for corrections computed --- src/compass/s1_cslc_qa.py | 22 ++---- src/compass/s1_geocode_slc.py | 27 ++++---- src/compass/utils/h5_helpers.py | 27 -------- src/compass/utils/lut.py | 116 ++++++++++++++++++++++++-------- 4 files changed, 105 insertions(+), 87 deletions(-) diff --git a/src/compass/s1_cslc_qa.py b/src/compass/s1_cslc_qa.py index 5cd89a93..7671a565 100644 --- a/src/compass/s1_cslc_qa.py +++ b/src/compass/s1_cslc_qa.py @@ -137,8 +137,7 @@ def compute_static_layer_stats(self, cslc_h5py_root, rdr2geo_params): static_layers) - def compute_correction_stats(self, cslc_h5py_root, apply_tropo_corrections, - tropo_delay_type): + def compute_correction_stats(self, cslc_h5py_root): ''' Compute correction stats. Stats written to HDF5 and saved to class dict for later JSON output @@ -147,26 +146,13 @@ def compute_correction_stats(self, cslc_h5py_root, apply_tropo_corrections, ---------- cslc_h5py_root: h5py.File Root of CSLC HDF5 - apply_tropo_corrections: bool - Whether or not to compute troposhpere correction stats - tropo_delay_type: str - Type of troposphere delay. Any between 'dry', or 'wet', or - 'wet_dry' for the sum of wet and dry troposphere delays. Only used - apply_tropo_corrections is true. ''' # path to source group corrections_src_path = f'{ROOT_PATH}/corrections' - # names of datasets to compute stats for - corrections = ['bistatic_delay', 'geometry_steering_doppler', - 'azimuth_fm_rate_mismatch', 'los_ionospheric_delay', - 'los_solid_earth_tides'] - - # check if tropo corrections need to be computed and saved - if apply_tropo_corrections: - for delay_type in ['wet', 'dry']: - if delay_type in tropo_delay_type: - corrections.append(f'{delay_type}_los_troposphere_delay') + # compute stats for corrections flagged true + corrections= [k for k, v in cslc_h5py_root[ + 'science/SENTINEL1/CSLC/metadata/processing_information/corrections'].items() if v] self.compute_stats_from_float_hdf5_dataset(cslc_h5py_root, corrections_src_path, diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index 49978ffd..02fd4828 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -72,13 +72,6 @@ def run(cfg: GeoRunConfig): # Create scratch as needed scratch_path = out_paths.scratch_directory - # Compute correction LUTs - rg_lut, az_lut = correction_luts(burst, cfg.lut_params, - dem_path=cfg.dem, - tec_path=cfg.tec_file, - scratch_path=scratch_path, - weather_model_path=cfg.weather_model_file) - radar_grid = burst.as_isce3_radargrid() native_doppler = burst.doppler.lut2d orbit = burst.orbit @@ -111,17 +104,26 @@ def run(cfg: GeoRunConfig): # add type to root for GDAL recognition of datasets ctype = h5py.h5t.py_create(np.complex64) - ctype.commit(geo_burst_h5['/'].id, np.string_('complex64')) + #ctype.commit(geo_burst_h5['/'].id, np.string_('complex64')) grid_path = f'{root_path}/CSLC/grids' grid_group = geo_burst_h5.require_group(grid_path) check_eap = is_eap_correction_necessary(burst.ipf_version) + + # Compute correction LUTs + rg_lut, az_lut = correction_luts(burst, cfg.lut_params, + dem_path=cfg.dem, + tec_path=cfg.tec_file, + h5_file_obj=geo_burst_h5, + scratch_path=scratch_path, + weather_model_path=cfg.weather_model_file) + for b in bursts: pol = b.polarization # Load the input burst SLC temp_slc_path = f'{scratch_path}/{out_paths.file_name_pol}_temp.vrt' - burst.slc_to_vrt_file(temp_slc_path) + #burst.slc_to_vrt_file(temp_slc_path) # Apply EAP correction if necessary if check_eap.phase_correction: @@ -189,13 +191,8 @@ def run(cfg: GeoRunConfig): # If needed, perform QA and write results to JSON if cfg.quality_assurance_params.perform_qa: cslc_qa = QualityAssuranceCSLC() - # if cfg.lut_params.enabled: - # # apply tropo corrections if weather file provided - # apply_tropo_corrections = cfg.weather_model_file is not None - # cslc_qa.compute_correction_stats( - # geo_burst_h5, apply_tropo_corrections, - # cfg.tropo_params.delay_type) cslc_qa.compute_CSLC_raster_stats(geo_burst_h5, bursts) + cslc_qa.compute_correction_stats(geo_burst_h5) cslc_qa.populate_rfi_dict(geo_burst_h5) cslc_qa.valid_pixel_percentages(geo_burst_h5) cslc_qa.set_orbit_type(cfg, geo_burst_h5) diff --git a/src/compass/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py index d69ac4cb..28a5f326 100644 --- a/src/compass/utils/h5_helpers.py +++ b/src/compass/utils/h5_helpers.py @@ -404,9 +404,6 @@ def metadata_to_h5group(parent_group, burst, cfg): cfg: types.SimpleNamespace SimpleNamespace containing run configuration ''' - if 'metadata' in parent_group: - del parent_group['metadata'] - # create metadata group to write datasets to meta_group = parent_group.require_group('metadata') @@ -626,31 +623,7 @@ def corrections_to_h5group(parent_group, burst, rg_lut, {'units': 'seconds'}), Meta('zero_doppler_time_spacing',rg_lut.y_spacing, 'spacing of azimuth time of LUT data', {'units': 'seconds'}), - Meta('geometry_steering_doppler', ds.GetRasterBand(1).ReadAsArray(), - f'geometry steering doppler (range) {desc}', - {'units': 'meters'}), - Meta('bistatic_delay', ds.GetRasterBand(2).ReadAsArray(), - f'bistatic delay (azimuth) {desc}', {'units': 'seconds'}), - Meta('azimuth_fm_rate_mismatch', ds.GetRasterBand(3).ReadAsArray(), - f'azimuth FM rate mismatch mitigation (azimuth) {desc}', - {'units': 'seconds'}), - Meta('los_solid_earth_tides', ds.GetRasterBand(4).ReadAsArray(), - f'Solid Earth tides (range) {desc}', - {'units': 'meters'}), - Meta('azimuth_solid_earth_tides', ds.GetRasterBand(5).ReadAsArray(), - f'Solid Earth tides (range) {desc}', - {'units': 'seconds'}), - Meta('los_ionospheric_delay', ds.GetRasterBand(7).ReadAsArray(), - f'Ionospheric delay (range) {desc}', - {'units': 'meters'}), - Meta('wet_los_troposphere_delay', ds.GetRasterBand(8).ReadAsArray(), - f'Wet LOS troposphere delay {desc}', - {'units': 'meters'}), - Meta('dry_los_troposphere_delay', ds.GetRasterBand(9).ReadAsArray(), - f'Dry LOS troposphere delay {desc}', - {'units': 'meters'}) ] - for meta_item in correction_items: add_dataset_and_attrs(correction_group, meta_item) diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 83c6a726..146172ad 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -10,12 +10,12 @@ from skimage.transform import resize from compass.utils.geometry_utils import enu2los, en2az -from compass.utils.iono import ionosphere_delay +from compass.utils.h5_helpers import Meta, add_dataset_and_attrs from compass.utils.helpers import open_raster -from compass.utils.helpers import write_raster +from compass.utils.iono import ionosphere_delay -def correction_luts(burst, lut_par, dem_path, tec_path, +def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, scratch_path=None, weather_model_path=None): ''' @@ -80,7 +80,7 @@ def correction_luts(burst, lut_par, dem_path, tec_path, rg_data = np.zeros(lut_shape, dtype=np.float32) az_data = np.zeros(lut_shape, dtype=np.float32) - # Initialize data list and description to save corrections + # Initialize data list and lut_description to save corrections data_dict_key_dscrs = ( ['doppler', 'Slant range geometry and steering doppler'], ['bistatic_delay', 'Bistatic delay'], @@ -94,23 +94,49 @@ def correction_luts(burst, lut_par, dem_path, tec_path, data_dict = {key: (np.zeros_like(rg_data), dscr) for (key, dscr) in data_dict_key_dscrs} + # Dict of meta correction items to be written to HDF5 + correction_lut_items = [] + + # Common string to all lut_descriptions + lut_desc = 'correction as a function of slant range and azimuth time' + + # Dict indicating if a correction item has been applied + correction_application_items = [] + + # Common string to all lut_descriptions + corr_desc = 'correction has been applied' + # Check which corrections are requested and accumulate corresponding data # Geometrical and steering doppler + correction_application_items.append( + Meta('geometry_steering_doppler', lut_par.geometry_steering_doppler, + f'Boolean if geometry steering doppler {corr_desc}')) if lut_par.geometry_steering_doppler: doppler = burst.doppler_induced_range_shift(range_step=rg_step, az_step=az_step) doppler_meter = doppler.data * isce3.core.speed_of_light * 0.5 rg_data += doppler_meter - data_dict['doppler'][0] = doppler_meter + correction_lut_items.append( + Meta('geometry_steering_doppler', doppler_meter, + f'geometry steering doppler (range) {lut_desc}', + {'units': 'meters'})) # Bistatic delay + correction_application_items.append( + Meta('bistatic_delay', lut_par.bistatic_delay, + f'Boolean if bistatic delay {corr_desc}')) if lut_par.bistatic_delay: bistatic_delay = burst.bistatic_delay(range_step=rg_step, az_step=az_step).data az_data -= bistatic_delay - data_dict['bistatic_delay'][0] = -bistatic_delay + correction_lut_items.append( + Meta('bistatic_delay', bistatic_delay, + f'bistatic delay (azimuth) {lut_desc}', {'units': 'seconds'})) # Azimuth FM-rate mismatch + correction_application_items.append( + Meta('azimuth_fm_rate_mismatch', lut_par.azimuth_fm_rate, + f'Boolean if azimuth FM rate mismatch mitigation {corr_desc}')) if lut_par.azimuth_fm_rate: az_fm_rate = burst.az_fm_rate_mismatch_from_llh(lat, lon, height, ellipsoid, @@ -118,9 +144,18 @@ def correction_luts(burst, lut_par, dem_path, tec_path, az_step=az_step, rg_step=rg_step)).data az_data -= az_fm_rate - data_dict['azimuth_fm_rate'][0] = -az_fm_rate + correction_lut_items.append( + Meta('azimuth_fm_rate_mismatch', az_fm_rate, + f'azimuth FM rate mismatch mitigation (azimuth) {lut_desc}', + {'units': 'seconds'})) # Solid Earth tides + correction_application_items.append( + Meta('los_solid_earth_tides', lut_par.solid_earth_tides, + f'Boolean if LOS solid Earth tides {corr_desc}')) + correction_application_items.append( + Meta('azimuth_solid_earth_tides', lut_par.solid_earth_tides, + f'Boolean if azimuth solid Earth tides {corr_desc}')) if lut_par.solid_earth_tides: dec_factor = int(np.round(5000.0 / rg_step)) dec_slice = np.s_[::dec_factor] @@ -136,30 +171,52 @@ def correction_luts(burst, lut_par, dem_path, tec_path, az_set = resize(az_set_temp, lut_shape, **kwargs) rg_data += rg_set az_data += az_set - data_dict['rg_set'][0] = rg_set - data_dict['az_set'][0] = az_set + correction_lut_items.append( + Meta('los_solid_earth_tides', rg_set, + f'Solid Earth tides (range) {lut_desc}', {'units': 'meters'})) + correction_lut_items.append( + Meta('azimuth_solid_earth_tides', az_set, + f'Solid Earth tides (azimuth) {lut_desc}', {'units': 'seconds'})) # Static troposphere + correction_application_items.append( + Meta('static_los_tropospheric_delay', lut_par.static_troposphere, + f'Boolean if static tropospheric delay {corr_desc}')) if lut_par.static_troposphere: los_static_tropo = compute_static_troposphere_delay(inc_angle, height) rg_data += los_static_tropo - data_dict['static_tropo'][0] = los_static_tropo + correction_lut_items.append( + Meta('static_los_tropospheric_delay', los_static_tropo, + f'Static tropospheric delay (range) {lut_desc}', + {'units': 'meters'})) # Ionosphere TEC correction + correction_application_items.append( + Meta('los_ionospheric_delay', lut_par.ionosphere_tec, + f'Boolean if ionospheric delay {corr_desc}')) if lut_par.ionosphere_tec: los_iono = ionosphere_delay(burst.sensing_mid, burst.wavelength, tec_path, lon, lat, inc_angle) rg_data += los_iono - data_dict['tec_iono'][0] = los_iono + correction_lut_items.append( + Meta('los_ionospheric_delay', los_iono, + f'Ionospheric delay (range) {lut_desc}', {'units': 'meters'})) # Weather model troposphere correction + tropo_enabled = lut_par.weather_model_troposphere.enabled + delay_type = lut_par.weather_model_troposphere.delay_type + correction_application_items.append( + Meta('wet_los_troposphere_delay', tropo_enabled and 'wet' in delay_type, + f'Boolean if wet LOS troposphere delay {corr_desc}')) + correction_application_items.append( + Meta('dry_los_troposphere_delay', tropo_enabled and 'dry' in delay_type, + f'Boolean if dry LOS troposphere delay {corr_desc}')) if lut_par.weather_model_troposphere.enabled: from RAiDER.delay import tropo_delay from RAiDER.llreader import RasterRDR from RAiDER.losreader import Zenith - delay_type = lut_par.weather_model_troposphere.delay_type # Instantiate an "aoi" object to read lat/lon/height files aoi = RasterRDR(rdr2geo_raster_paths[1], rdr2geo_raster_paths[0], rdr2geo_raster_paths[2]) @@ -172,17 +229,32 @@ def correction_luts(burst, lut_par, dem_path, tec_path, 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)) - if 'wet' in delay_type: + wet_los_tropo = 2.0 * zen_wet / np.cos(np.deg2rad(inc_angle)) rg_data += wet_los_tropo - data_dict['wet_tropo'][0] = wet_los_tropo + correction_lut_items.append( + Meta('wet_los_troposphere_delay', wet_los_tropo, + f'Wet LOS troposphere delay {lut_desc}', + {'units': 'meters'})) if 'dry' in delay_type: + dry_los_tropo = 2.0 * zen_dry / np.cos(np.deg2rad(inc_angle)) rg_data += dry_los_tropo - data_dict['dry_tropo'][0] = dry_los_tropo + correction_lut_items.append( + Meta('dry_los_troposphere_delay', dry_los_tropo, + f'Dry LOS troposphere delay {lut_desc}', + {'units': 'meters'})) + + proc_nfo_group = \ + h5_file_obj.require_group('science/SENTINEL1/CSLC/metadata/processing_information/corrections') + for meta_item in correction_application_items: + add_dataset_and_attrs(proc_nfo_group, meta_item) + + correction_group = h5_file_obj.require_group('science/SENTINEL1/CSLC/corrections') + for meta_item in correction_lut_items: + add_dataset_and_attrs(correction_group, meta_item) # Create the range and azimuth LUT2d rg_lut = isce3.core.LUT2d(lut.x_start, lut.y_start, @@ -192,16 +264,6 @@ def correction_luts(burst, lut_par, dem_path, tec_path, lut.x_spacing, lut.y_spacing, az_data) - # Save corrections - driver = gdal.GetDriverByName('ENVI') - out_ds = driver.Create(f'{output_path}/corrections', - lut_shape[1], lut_shape[0], len(data_dict), - gdal.GDT_Float32) - for band, (key, (data, desc)) in enumerate(data_dict.items()): - raster_band = out_ds.GetRasterBand(band + 1) - raster_band.SetDescription(desc) - raster_band.WriteArray(data) - return rg_lut, az_lut From 2ac7fedbd5e9947c849a009c5e59d76a0f0429ce Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Mon, 22 May 2023 12:45:02 -0700 Subject: [PATCH 12/16] restore data type assignment path updates fix dataset access to correction flags update documentation remove unused code --- src/compass/s1_cslc_qa.py | 14 ++++++++++---- src/compass/s1_geocode_slc.py | 2 +- src/compass/utils/h5_helpers.py | 5 +---- src/compass/utils/lut.py | 23 +++++------------------ 4 files changed, 17 insertions(+), 27 deletions(-) diff --git a/src/compass/s1_cslc_qa.py b/src/compass/s1_cslc_qa.py index e0138189..2d13045f 100644 --- a/src/compass/s1_cslc_qa.py +++ b/src/compass/s1_cslc_qa.py @@ -7,7 +7,7 @@ import isce3 import numpy as np -from compass.utils.h5_helpers import (DATA_PATH, METADATA_PATH, +from compass.utils.h5_helpers import (DATA_PATH, PROCESSING_INFO_PATH, QA_PATH, add_dataset_and_attrs, Meta) @@ -148,15 +148,16 @@ def compute_correction_stats(self, cslc_h5py_root): Root of CSLC HDF5 ''' # path to source group - corrections_src_path = f'{METADATA_PATH}/processing_information/timing_corrections' + corrections_src_path = f'{DATA_PATH}/timing_corrections' # compute stats for corrections flagged true corrections= [k for k, v in cslc_h5py_root[ - 'science/SENTINEL1/CSLC/metadata/processing_information/corrections'].items() if v] + f'{PROCESSING_INFO_PATH}/corrections'].items() if v[()]] self.compute_stats_from_float_hdf5_dataset(cslc_h5py_root, corrections_src_path, - 'timing_corrections', corrections) + 'timing_corrections', + corrections) def compute_stats_from_float_hdf5_dataset(self, cslc_h5py_root, @@ -171,8 +172,13 @@ def compute_stats_from_float_hdf5_dataset(self, cslc_h5py_root, cslc_h5py_root: h5py.File Root of CSLC HDF5 src_group_path: str + Path to HDF5 group with datasets whose stats are to be computed qa_group_name: str + Group to be created in QA statistics HDF5 group to contain stats + from datasts in src_group_path qa_item_names: list[str] + Names of datasets in src_group_path path whose stats are to be + computed ''' # init dict to save all QA item stats to self.stats_dict[qa_group_name] = {} diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index cf55757a..0918dddc 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -107,7 +107,7 @@ def run(cfg: GeoRunConfig): # add type to root for GDAL recognition of datasets ctype = h5py.h5t.py_create(np.complex64) - #ctype.commit(geo_burst_h5['/'].id, np.string_('complex64')) + ctype.commit(geo_burst_h5['/'].id, np.string_('complex64')) grid_group = geo_burst_h5.require_group(DATA_PATH) check_eap = is_eap_correction_necessary(burst.ipf_version) diff --git a/src/compass/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py index 5a56251d..5d839648 100644 --- a/src/compass/utils/h5_helpers.py +++ b/src/compass/utils/h5_helpers.py @@ -21,6 +21,7 @@ DATA_PATH = '/data' QA_PATH = '/quality_assurance' METADATA_PATH = '/metadata' +PROCESSING_INFO_PATH = f'{METADATA_PATH}/processing_information' @dataclass @@ -619,10 +620,7 @@ def corrections_to_h5group(parent_group, burst, rg_lut, scratch_path: str Path to the scratch directory ''' - # Open GDAL dataset to fetch corrections - ds = gdal.Open(f'{scratch_path}/corrections/corrections', - gdal.GA_ReadOnly) correction_group = parent_group.require_group('timing_corrections') # create slant range and azimuth vectors shared by the LUTs @@ -634,7 +632,6 @@ def corrections_to_h5group(parent_group, burst, rg_lut, az_lut.length, dtype=np.float64) # correction LUTs axis and doppler correction LUTs - desc = 'correction as a function of slant range and azimuth time' correction_items = [ Meta('slant_range', slant_range, 'slant range of LUT data', {'units': 'meters'}), diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 146172ad..b772c2d7 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -10,7 +10,8 @@ from skimage.transform import resize from compass.utils.geometry_utils import enu2los, en2az -from compass.utils.h5_helpers import Meta, add_dataset_and_attrs +from compass.utils.h5_helpers import (Meta, add_dataset_and_attrs, + DATA_PATH, PROCESSING_INFO_PATH) from compass.utils.helpers import open_raster from compass.utils.iono import ionosphere_delay @@ -80,20 +81,6 @@ def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, rg_data = np.zeros(lut_shape, dtype=np.float32) az_data = np.zeros(lut_shape, dtype=np.float32) - # Initialize data list and lut_description to save corrections - data_dict_key_dscrs = ( - ['doppler', 'Slant range geometry and steering doppler'], - ['bistatic_delay', 'Bistatic delay'], - ['azimuth_fm_rate', 'Azimuth FM rate mismatch'], - ['rg_set', 'Slant range Solid Earth Tides'], - ['az_set', 'Azimuth Solid Earth Tides'], - ['static_tropo', 'Static troposphere delay'], - ['tec_iono', 'Slant range ionosphere delay'], - ['dry_tropo', 'Dry troposphere delay from weather model'], - ['wet_tropo', 'Wet troposphere delay from weather model']) - data_dict = {key: (np.zeros_like(rg_data), dscr) - for (key, dscr) in data_dict_key_dscrs} - # Dict of meta correction items to be written to HDF5 correction_lut_items = [] @@ -212,7 +199,7 @@ def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, correction_application_items.append( Meta('dry_los_troposphere_delay', tropo_enabled and 'dry' in delay_type, f'Boolean if dry LOS troposphere delay {corr_desc}')) - if lut_par.weather_model_troposphere.enabled: + if tropo_enabled: from RAiDER.delay import tropo_delay from RAiDER.llreader import RasterRDR from RAiDER.losreader import Zenith @@ -248,11 +235,11 @@ def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, {'units': 'meters'})) proc_nfo_group = \ - h5_file_obj.require_group('science/SENTINEL1/CSLC/metadata/processing_information/corrections') + h5_file_obj.require_group(f'{PROCESSING_INFO_PATH}/corrections') for meta_item in correction_application_items: add_dataset_and_attrs(proc_nfo_group, meta_item) - correction_group = h5_file_obj.require_group('science/SENTINEL1/CSLC/corrections') + correction_group = h5_file_obj.require_group(f'{DATA_PATH}/timing_corrections') for meta_item in correction_lut_items: add_dataset_and_attrs(correction_group, meta_item) From 31935c4f35bf80340641c67f905d80734591a61f Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Mon, 22 May 2023 14:12:58 -0700 Subject: [PATCH 13/16] replace grid with data --- tests/conftest.py | 4 ++-- tests/test_s1_geocode_slc.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 780861e2..3bd31e02 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -93,8 +93,8 @@ def geocode_slc_params(): test_params.output_hdf5_path = f'{output_path}/{output_file_name}' # path to groups and datasets in output HDF5 - test_params.grid_group_path = DATA_PATH - test_params.raster_path = f'{test_params.grid_group_path}/VV' + test_params.data_group_path = DATA_PATH + test_params.raster_path = f'{test_params.data_group_path}/VV' return test_params diff --git a/tests/test_s1_geocode_slc.py b/tests/test_s1_geocode_slc.py index 352fc8cc..1e043480 100644 --- a/tests/test_s1_geocode_slc.py +++ b/tests/test_s1_geocode_slc.py @@ -56,7 +56,7 @@ def _get_reflectors_bounding_slice(geocode_slc_params): ''' # extract from HDF5 with h5py.File(geocode_slc_params.output_hdf5_path, 'r') as h5_obj: - grid_group = h5_obj[geocode_slc_params.grid_group_path] + grid_group = h5_obj[geocode_slc_params.data_group_path] # create projection to covert from UTM to LLH epsg = int(grid_group['projection'][()]) From a6dc9876fd1ea9290a3ec4a9e86c5032e0eb0b16 Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Tue, 23 May 2023 09:13:01 -0700 Subject: [PATCH 14/16] uncomment vrt writing --- requirements.txt | 2 +- src/compass/s1_geocode_slc.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index 1f23ea1d..bead4130 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ python>=3.9 # sentinel1-reader requirement numpy # sentinel1-reader requirement lxml # sentinel1-reader requirement gdal>=3 -#isce3 # since the conda-installed isce3 is not the most updated version, installing isce3 from stratch is recommended, to stay in sync with isce3 development. +isce3 # since the conda-installed isce3 is not the most updated version, installing isce3 from stratch is recommended, to stay in sync with isce3 development. #journal # as of Mar 2022, journal from conda does not support python3.9; since it is included during isce3 installation above, comment this out temporarily. pandas pyproj diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index 0918dddc..67d683c2 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -125,7 +125,7 @@ def run(cfg: GeoRunConfig): # Load the input burst SLC temp_slc_path = f'{scratch_path}/{out_paths.file_name_pol}_temp.vrt' - #burst.slc_to_vrt_file(temp_slc_path) + burst.slc_to_vrt_file(temp_slc_path) # Apply EAP correction if necessary if check_eap.phase_correction: From 679674c408343c38c7501aa28b6dadbfdc5e0130 Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Tue, 23 May 2023 09:34:13 -0700 Subject: [PATCH 15/16] fix paths update/add comments and variable names --- src/compass/s1_geocode_slc.py | 11 +++++------ src/compass/schemas/s1_cslc_geo.yaml | 2 +- src/compass/utils/h5_helpers.py | 7 ++++--- src/compass/utils/lut.py | 5 ++++- 4 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index 67d683c2..d35e7440 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -23,9 +23,9 @@ corrections_to_h5group, identity_to_h5group, init_geocoded_dataset, - METADATA_PATH, - metadata_to_h5group, - ROOT_PATH) + METADATA_PATH, PROCESSING_INFO_PATH, + ROOT_PATH, metadata_to_h5group, + ) from compass.utils.helpers import bursts_grouping_generator, get_module_name from compass.utils.lut import correction_luts from compass.utils.yaml_argparse import YamlArgparse @@ -112,7 +112,7 @@ def run(cfg: GeoRunConfig): grid_group = geo_burst_h5.require_group(DATA_PATH) check_eap = is_eap_correction_necessary(burst.ipf_version) - # Compute correction LUTs + # Get cumulative correction LUTs and save individual LUTs to HDF5 rg_lut, az_lut = correction_luts(burst, cfg.lut_params, dem_path=cfg.dem, tec_path=cfg.tec_file, @@ -175,8 +175,7 @@ def run(cfg: GeoRunConfig): identity_to_h5group(root_group, burst, cfg) metadata_to_h5group(root_group, burst, cfg) - correction_group = geo_burst_h5.require_group( - f'{METADATA_PATH}/processing_information') + correction_group = geo_burst_h5.require_group(PROCESSING_INFO_PATH) corrections_to_h5group(correction_group, burst, rg_lut, az_lut, scratch_path) diff --git a/src/compass/schemas/s1_cslc_geo.yaml b/src/compass/schemas/s1_cslc_geo.yaml index 663a07ed..e6c272b9 100644 --- a/src/compass/schemas/s1_cslc_geo.yaml +++ b/src/compass/schemas/s1_cslc_geo.yaml @@ -103,7 +103,7 @@ lut_options: ionosphere_tec: bool(required=False) # Enable/disable static troposphere correction static_troposphere: bool(required=False) - # Troposphere delay using weather model + # Options for troposphere delay using weather model weather_model_troposphere: include('troposphere_options', required=False) troposphere_options: diff --git a/src/compass/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py index 5d839648..18f44a6d 100644 --- a/src/compass/utils/h5_helpers.py +++ b/src/compass/utils/h5_helpers.py @@ -631,8 +631,9 @@ def corrections_to_h5group(parent_group, burst, rg_lut, azimuth = np.linspace(az_lut.y_start, y_end, az_lut.length, dtype=np.float64) - # correction LUTs axis and doppler correction LUTs - correction_items = [ + # correction LUTs axis info only. LUTs written as needed in + # compass.utils.lut.correction_luts + correction_axis_items = [ Meta('slant_range', slant_range, 'slant range of LUT data', {'units': 'meters'}), Meta('slant_range_spacing', rg_lut.x_spacing, @@ -642,7 +643,7 @@ def corrections_to_h5group(parent_group, burst, rg_lut, Meta('zero_doppler_time_spacing',rg_lut.y_spacing, 'spacing of azimuth time of LUT data', {'units': 'seconds'}), ] - for meta_item in correction_items: + for meta_item in correction_axis_items: add_dataset_and_attrs(correction_group, meta_item) # Extended FM rate and doppler centroid polynomial coefficients for azimuth diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 187ac634..3b44ee48 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -20,7 +20,8 @@ def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, scratch_path=None, weather_model_path=None): ''' - Compute correction look-up tables (LUTs) + Compute correction look-up tables (LUTs), write as needed to HDF5, then + return cumulative LUTs Parameters ---------- @@ -234,11 +235,13 @@ def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, f'Dry LOS troposphere delay {lut_desc}', {'units': 'meters'})) + # Write correction flags to HDF5 proc_nfo_group = \ h5_file_obj.require_group(f'{PROCESSING_INFO_PATH}/corrections') for meta_item in correction_application_items: add_dataset_and_attrs(proc_nfo_group, meta_item) + # Write LUTs of enabled corrections to HDF5 correction_group = h5_file_obj.require_group(f'{DATA_PATH}/timing_corrections') for meta_item in correction_lut_items: add_dataset_and_attrs(correction_group, meta_item) From 0b128608a1403bb7c5acc1423f1278a81fe98660 Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Wed, 24 May 2023 08:39:34 -0700 Subject: [PATCH 16/16] Revert "Merge branch 'main' into single_lut" This reverts commit 3c4d8d40ccc5f4654643f8562c31ebc4a37637f4, reversing changes made to 31935c4f35bf80340641c67f905d80734591a61f. --- src/compass/utils/geometry_utils.py | 163 +--------------------------- src/compass/utils/lut.py | 15 ++- 2 files changed, 11 insertions(+), 167 deletions(-) diff --git a/src/compass/utils/geometry_utils.py b/src/compass/utils/geometry_utils.py index aff62526..6711679d 100644 --- a/src/compass/utils/geometry_utils.py +++ b/src/compass/utils/geometry_utils.py @@ -7,7 +7,7 @@ import numpy as np -import isce3 + def los2orbit_azimuth_angle(los_az_angle, look_direction='right'): """ @@ -265,164 +265,3 @@ def get_unit_vector4component_of_interest(los_inc_angle, los_az_angle, comp='enu ] return unit_vec - - -def enu2rgaz(radargrid_ref, orbit, ellipsoid, - lon_arr, lat_arr, hgt_arr, - e_arr, n_arr, u_arr, - geo2rdr_params=None): - ''' - Convert ENU displacement into range / azimuth displacement, - based on the idea mentioned in ETAD ATBD, available in the link below: - https://sentinels.copernicus.eu/documents/247904/4629150/ETAD-DLR-DD-0008_Algorithm-Technical-Baseline-Document_2.3.pdf/5cb45b43-76dc-8dec-04ef-ca1252ace434?t=1680181574715 # noqa - - Algorithm description - --------------------- - For all lon / lat / height of the array; - 1. Calculate the ECEF coordinates before applying SET - 2. Calculate the unit vector of east / north / up directions of the point (i.e. ENU vectors) - 3. Scale the ENU vectors in 2 with ENU displacement to - get the displacement in ECEF - 4. Add the vectors calculated in 3 into 1. - This will be the ECEF coordinates after applying SET - 5. Convert 4 into lat / lon / hgt. - This will be LLH coordinates after applying SET - 6. Calculate the radar coordinate before SET applied using `geo2rdr` - 7. Calculate the radar coordinate AFTER SET applied using `geo2rdr` - 8. Calculate the difference between (7) and (6), - which will be the displacement in radargrid by SET - - Parameters - ---------- - radargrid_ref: isce3.product.RadarGridParameters - Radargrid of the burst - orbit: isce3.core.Orbit - Orbit of the burst - ellipsoid: isce3.core.Ellipsoid - Ellipsoid definition - lon_arr, lat_arr, hgt_arr: np.nadrray - Arrays for longitude, latitude, and height. - Units for longitude and latitude are degree; unit for height is meters. - e_arr, n_arr, u_arr: np.ndarray - Displacement in east, north, and up direction in meters - geo2rdr_params: SimpleNameSpace - Parameters for geo2rdr - - Returns - ------- - rg_arr: np.ndarray - Displacement in slant range direction in meters. - az_arr: np.ndarray - Displacement in azimuth direction in seconds. - - Notes - ----- - When `geo2rdr_params` is not provided, then the iteration - threshold and max # iterations are set to - `1.0e-8` and `25` respectively. - - ''' - if geo2rdr_params is None: - # default threshold and # iteration for geo2rdr - threshold = 1.0e-8 - maxiter = 25 - else: - threshold = geo2rdr_params.threshold - maxiter = geo2rdr_params.numiter - - shape_arr = lon_arr.shape - rg_arr = np.zeros(shape_arr) - az_arr = np.zeros(shape_arr) - - # Calculate the ENU vector in ECEF - for i, lon_deg in enumerate(np.nditer(lon_arr)): - index_arr = np.unravel_index(i, lon_arr.shape) - lat_deg = lat_arr[index_arr] - hgt = hgt_arr[index_arr] - - vec_e, vec_n, vec_u = get_enu_vector_ecef(lon_deg, lat_deg) - - llh_ref = np.array([np.deg2rad(lon_deg), - np.deg2rad(lat_deg), - hgt]) - - xyz_before = ellipsoid.lon_lat_to_xyz(llh_ref) - xyz_after_set = (xyz_before - + vec_e * e_arr[index_arr] - + vec_n * n_arr[index_arr] - + vec_u * u_arr[index_arr]) - llh_displaced = ellipsoid.xyz_to_lon_lat(xyz_after_set) - - aztime_ref, slant_range_ref =\ - isce3.geometry.geo2rdr(llh_ref, - ellipsoid, - orbit, - isce3.core.LUT2d(), - radargrid_ref.wavelength, - radargrid_ref.lookside, - threshold=threshold, - maxiter=maxiter) - - aztime_displaced, slant_range_displaced =\ - isce3.geometry.geo2rdr(llh_displaced, - ellipsoid, - orbit, - isce3.core.LUT2d(), - radargrid_ref.wavelength, - radargrid_ref.lookside, - threshold=threshold, - maxiter=maxiter) - - rg_arr[index_arr] = slant_range_displaced - slant_range_ref - az_arr[index_arr] = aztime_displaced - aztime_ref - - return rg_arr, az_arr - - -def get_enu_vector_ecef(lon, lat, units='degrees'): - ''' - Calculate the east, north, and up vectors in ECEF for lon / lat provided - - Parameters - ---------- - lon: np.ndarray - Longitude of the points to calculate ENU vectors - lat: np.ndarray - Latitude of the points to calculate ENU vectors - units: str - Units of the `lon` and `lat`. - Acceptable units are `radians` or `degrees`, (Default: degrees) - - Returns - ------- - vec_e: np.ndarray - unit vector of "east" direction in ECEF - vec_n: np.ndarray - unit vector of "north" direction in ECEF - vec_u: np.ndarray - unit vector of "up" direction in ECEF - ''' - if units == 'degrees': - lon_rad = np.deg2rad(lon) - lat_rad = np.deg2rad(lat) - elif units == 'radians': - lon_rad = lon - lat_rad = lat - else: - raise ValueError(f'"{units}" was provided for `units`, ' - 'which needs to be either `degrees` or `radians`') - - # Calculate up, north, and east vectors - # reference: https://github.com/isce-framework/isce3/blob/944eba17f4a5b1c88c6a035c2d58ddd0d4f0709c/cxx/isce3/core/Ellipsoid.h#L154-L157 # noqa - # https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU # noqa - vec_u = np.array([np.cos(lon_rad) * np.cos(lat_rad), - np.sin(lon_rad) * np.cos(lat_rad), - np.sin(lat_rad)]) - - vec_n = np.array([-np.cos(lon_rad) * np.sin(lat_rad), - -np.sin(lon_rad) * np.sin(lat_rad), - np.cos(lat_rad)]) - - vec_e = np.cross(vec_n, vec_u, axis=0) - - return vec_e, vec_n, vec_u diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 3b44ee48..4f583d28 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -256,8 +256,9 @@ def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, return rg_lut, az_lut -def solid_earth_tides(burst, lat_radar_grid, lon_radar_grid, hgt_radar_grid, - ellipsoid, geo2rdr_params=None): + +def solid_earth_tides(burst, lat_radar_grid, lon_radar_grid, inc_angle, + head_angle): ''' Compute displacement due to Solid Earth Tides (SET) in slant range and azimuth directions @@ -325,9 +326,13 @@ def solid_earth_tides(burst, lat_radar_grid, lon_radar_grid, hgt_radar_grid, for set_enu in [set_e, set_n, set_u]] # Convert SET from ENU to range/azimuth coordinates - set_rg, set_az = enu2rgaz(burst.as_isce3_radargrid(), burst.orbit, ellipsoid, - lon_radar_grid, lat_radar_grid, hgt_radar_grid, - rdr_set_e, rdr_set_n, rdr_set_u, geo2rdr_params) + # Note: rdr2geo heading angle is measured wrt to the East and it is positive + # anti-clockwise. To convert ENU to LOS, we need the azimuth angle which is + # measured from the north and positive anti-clockwise + # azimuth_angle = heading + 90 + set_rg = enu2los(rdr_set_e, rdr_set_n, rdr_set_u, inc_angle, + az_angle=head_angle + 90.0) + set_az = en2az(rdr_set_e, rdr_set_n, head_angle - 90.0) return set_rg, set_az