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/defaults/s1_cslc_geo.yaml b/src/compass/defaults/s1_cslc_geo.yaml index c297a267..7ae27450 100644 --- a/src/compass/defaults/s1_cslc_geo.yaml +++ b/src/compass/defaults/s1_cslc_geo.yaml @@ -57,18 +57,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_tides: True + # Enable/disable ionosphere TEC correction + ionosphere_tec: 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/s1_cslc_qa.py b/src/compass/s1_cslc_qa.py index a73f47e1..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) @@ -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,30 +146,18 @@ 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'{METADATA_PATH}/processing_information/timing_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'] + corrections_src_path = f'{DATA_PATH}/timing_corrections' - # 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[ + 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, @@ -185,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 431d36ca..d35e7440 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -23,11 +23,11 @@ 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 cumulative_correction_luts +from compass.utils.lut import correction_luts from compass.utils.yaml_argparse import YamlArgparse @@ -75,22 +75,6 @@ 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, - geo2rdr_params=cfg.geo2rdr_params) - else: - rg_lut = isce3.core.LUT2d() - az_lut = isce3.core.LUT2d() - radar_grid = burst.as_isce3_radargrid() native_doppler = burst.doppler.lut2d orbit = burst.orbit @@ -127,6 +111,15 @@ def run(cfg: GeoRunConfig): grid_group = geo_burst_h5.require_group(DATA_PATH) check_eap = is_eap_correction_necessary(burst.ipf_version) + + # 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, + h5_file_obj=geo_burst_h5, + scratch_path=scratch_path, + weather_model_path=cfg.weather_model_file) + for b in bursts: pol = b.polarization @@ -182,13 +175,9 @@ def run(cfg: GeoRunConfig): identity_to_h5group(root_group, burst, cfg) metadata_to_h5group(root_group, burst, cfg) - if cfg.lut_params.enabled: - correction_group = geo_burst_h5.require_group( - f'{METADATA_PATH}/processing_information') - corrections_to_h5group(correction_group, burst, cfg, rg_lut, az_lut, - scratch_path, - weather_model_path=cfg.weather_model_file, - delay_type=cfg.tropo_params.delay_type) + correction_group = geo_burst_h5.require_group(PROCESSING_INFO_PATH) + corrections_to_h5group(correction_group, burst, rg_lut, az_lut, + scratch_path) # If needed, make browse image and compute CSLC raster stats browse_params = cfg.browse_image_params @@ -202,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/s1_geocode_stack.py b/src/compass/s1_geocode_stack.py index 620d18cc..6b057ffe 100755 --- a/src/compass/s1_geocode_stack.py +++ b/src/compass/s1_geocode_stack.py @@ -70,8 +70,21 @@ def create_parser(): optional.add_argument('-m', '--metadata', action='store_true', help='If flag is set, generates radar metadata layers for each' ' burst stack (see rdr2geo processing options)') - optional.add_argument('-nc', '--no-corrections', action='store_true', - help='If flag is set, skip the geocoding LUT corrections.') + optional.add_argument('-gsd', '--geometry-steering-doppler-correction', action='store_false', + help='If flag is set, skip application of geometry steering doppler.') + optional.add_argument('-bd', '--bistatic-delay-correction', action='store_false', + help='If flag is set, skip application of bistatic delay.') + optional.add_argument('-afr', '--azimuth-fm-rate-correction', action='store_false', + help='If flag is set, skip application of azimuth FM rate.') + optional.add_argument('-set', '--solid-earth-tides-correction', action='store_false', + help='If flag is set, skip application of solid earth tides.') + optional.add_argument('-it', '--ionosphere-tec-correction', action='store_false', + help='If flag is set, skip application of ionosphere TEC.') + optional.add_argument('-st', '--static-troposphere-correction', action='store_false', + help='If flag is set, skip application of static troposphere.') + optional.add_argument('-wmt', '--weather-model-troposphere-correction', type=str, default=None, + choices=['dry', 'wet', 'wet_dry'], + help='If no value is given, weather model troposphere not applied. Otherwise apply based on given type.') optional.add_argument('--unzipped', action='store_true', help='If flag is set, assumes that the SLCs are unzipped, ' 'and only the SAFE directory is provided.') @@ -226,7 +239,7 @@ def get_common_burst_ids(data): def create_runconfig(burst_map_row, dem_file, work_dir, flatten, pol, x_spac, - y_spac, enable_metadata, enable_corrections, burst_db_file): + y_spac, enable_metadata, corrections, burst_db_file): """ Create runconfig to process geocoded bursts @@ -248,10 +261,10 @@ def create_runconfig(burst_map_row, dem_file, work_dir, flatten, pol, x_spac, Spacing of geocoded burst along Y-direction enable_metadata: bool Flag to enable/disable metadata generation for each burst stack. + corrections: dict + Dict containing flags to enable/disable applying corrections to burst stacks. burst_db_file: str Path to burst database file to use for burst bounding boxes. - enable_corrections: bool - Flag to enable/disable applying corrections to burst stacks. Returns ------- @@ -284,7 +297,7 @@ def create_runconfig(burst_map_row, dem_file, work_dir, flatten, pol, x_spac, # Geocoding process['polarization'] = pol - process['correction_luts']['enabled'] = enable_corrections + process['correction_luts'] = corrections geocode['flatten'] = flatten geocode['x_posting'] = x_spac geocode['y_posting'] = y_spac @@ -480,11 +493,26 @@ def run(slc_dir, dem_file, burst_id=None, common_bursts_only=False, start_date=N print('Elapsed time (min):', (end_time - start_time) / 60.0) +def _get_correction_args(args): + # gather correction related args into dict for easy iteration + corrections = {k:v for k, v in args.__dict__.items() if 'correction' in k} + wmt = corrections['weather_model_troposphere_correction'] + if wmt is not None: + corrections['weather_model_troposphere_correction'] = \ + {'enabled':True, 'delay_type':wmt} + else: + corrections['weather_model_troposphere_correction'] = \ + {'enabled':False, 'delay_type':wmt} + return corrections + def main(): """Create the command line interface and run the script.""" # Run main script args = create_parser() + # gather correction related args for easy grouped iteration later + corrections = _get_correction_args(args) + run( slc_dir=args.slc_dir, dem_file=args.dem_file, @@ -504,7 +532,7 @@ def main(): burst_db_file=args.burst_db_file, flatten=not args.no_flatten, enable_metadata=args.metadata, - enable_corrections=not args.no_corrections, + corrections=corrections, using_zipped=not args.unzipped, ) diff --git a/src/compass/schemas/s1_cslc_geo.yaml b/src/compass/schemas/s1_cslc_geo.yaml index 35110317..e6c272b9 100644 --- a/src/compass/schemas/s1_cslc_geo.yaml +++ b/src/compass/schemas/s1_cslc_geo.yaml @@ -87,17 +87,27 @@ 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) - # Troposphere delay using weather model - troposphere: include('troposphere_options', required=False) + # Enable/disable geometry steering doppler correction + geometry_steering_doppler: bool(required=False) + # Enable/disable bistatic delay correction + bistatic_delay: bool(required=False) + # Enable/disable azimuth FM rate mismatch correction + azimuth_fm_rate: bool(required=False) + # Enable/disable Solid Earth tides correction + solid_earth_tides: bool(required=False) + # Enable/disable ionosphere TEC correction + ionosphere_tec: bool(required=False) + # Enable/disable static troposphere correction + static_troposphere: bool(required=False) + # Options for troposphere delay using weather model + weather_model_troposphere: include('troposphere_options', 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) 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/h5_helpers.py b/src/compass/utils/h5_helpers.py index c8b1d7c1..18f44a6d 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 @@ -414,9 +415,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') @@ -604,9 +602,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 @@ -616,31 +613,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) correction_group = parent_group.require_group('timing_corrections') # create slant range and azimuth vectors shared by the LUTs @@ -651,9 +631,9 @@ def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut, azimuth = np.linspace(az_lut.y_start, y_end, 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 = [ + # 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, @@ -662,37 +642,8 @@ 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('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 (azimuth) {desc}', - {'units': 'seconds'}), - Meta('los_ionospheric_delay', ds.GetRasterBand(6).ReadAsArray(), - f'Ionospheric delay (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(7).ReadAsArray(), - f'Wet LOS troposphere delay {desc}', - {'units': 'meters'})) - if 'dry' in delay_type: - correction_items.append(Meta('dry_los_troposphere_delay', - ds.GetRasterBand(8).ReadAsArray(), - f'Dry LOS troposphere delay {desc}', - {'units': 'meters'})) - - 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/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 diff --git a/src/compass/utils/lut.py b/src/compass/utils/lut.py index 461bb034..4f583d28 100644 --- a/src/compass/utils/lut.py +++ b/src/compass/utils/lut.py @@ -9,251 +9,202 @@ from scipy.interpolate import RegularGridInterpolator as RGI from skimage.transform import resize -from compass.utils.geometry_utils import enu2rgaz -from compass.utils.iono import ionosphere_delay +from compass.utils.geometry_utils import enu2los, en2az +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.helpers import write_raster +from compass.utils.iono import ionosphere_delay -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', - geo2rdr_params=None): +def correction_luts(burst, lut_par, dem_path, tec_path, h5_file_obj, + scratch_path=None, + weather_model_path=None): ''' - Sum correction LUTs and returns cumulative correction LUT in slant range - and azimuth directions + Compute correction look-up tables (LUTs), write as needed to HDF5, then + return cumulative LUTs Parameters ---------- burst: Sentinel1BurstSlc - Sentinel-1 A/B burst SLC object + S1-A/B burst SLC object + lut_par: dict + Dictionary with LUT parameters dem_path: str - Path to the DEM file + File path to DEM tec_path: str - Path to the TEC file in IONEX format + File path to ionosphere TEC file scratch_path: str - Path to the scratch directory + File path to scratch directory (default: None) 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. + File path to weather model file (default: None) Returns ------- rg_lut: isce3.core.LUT2d - Sum of slant range correction LUTs in meters as a function of azimuth - time and slant range + Cumulative LUT in slant range direction (meters) 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, tide_az], \ - 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, - geo2rdr_params=geo2rdr_params) - - # 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, tide_az, los_ionosphere] - descr = ['slant range geometrical doppler', 'azimuth bistatic delay', - 'azimuth FM rate mismatch', 'slant range Solid Earth tides', - 'azimuth time 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, - geo2rdr_params=None): + Cumulative LUT in azimuth direction (seconds) ''' - 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 + # 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) + + # 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 + 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 + 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, 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, ::dec_factor] - rg_set_temp, az_set_temp = solid_earth_tides(burst, - lat[dec_slice], - lon[dec_slice], - height[dec_slice], - ellipsoid, - geo2rdr_params) - - # 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 + 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] + 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, + 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 + 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) - - else: + rg_data += 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 + 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 tropo_enabled: from RAiDER.delay import tropo_delay from RAiDER.llreader import RasterRDR from RAiDER.losreader import Zenith + # Instantiate an "aoi" object to read lat/lon/height files aoi = RasterRDR(rdr2geo_raster_paths[1], rdr2geo_raster_paths[0], rdr2geo_raster_paths[2]) @@ -269,18 +220,45 @@ def compute_geocoding_correction_luts(burst, dem_path, tec_path, # 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 + 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 + correction_lut_items.append( + Meta('dry_los_troposphere_delay', dry_los_tropo, + 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) + + # 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) - return ( - geometrical_steering_doppler, bistatic_delay, az_fm_mismatch, - [rg_set, az_set], los_ionosphere, - [wet_los_tropo, dry_los_tropo], los_static_tropo - ) + 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 @@ -348,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 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'][()])