From 35b5c4f1366dde8a87e57395bf90fd887ace78dc Mon Sep 17 00:00:00 2001 From: BRAUN REMI Date: Tue, 2 Jan 2024 16:36:53 +0100 Subject: [PATCH] ENH: Better handling of calibration step in SNAP for SAR data --- CHANGES.md | 1 + .../data/grd_no_calib_preprocess_default.xml | 80 +++++++++++++++++++ eoreader/products/sar/cosmo_product.py | 2 +- eoreader/products/sar/csg_product.py | 14 ++++ eoreader/products/sar/csk_product.py | 15 ++++ eoreader/products/sar/sar_product.py | 71 ++++++++-------- eoreader/products/sar/tsx_product.py | 18 +++++ eoreader/utils.py | 34 +++++++- 8 files changed, 194 insertions(+), 41 deletions(-) create mode 100644 eoreader/data/grd_no_calib_preprocess_default.xml diff --git a/CHANGES.md b/CHANGES.md index f1bd9ead..5917ab47 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -11,6 +11,7 @@ - **ENH: Calibration step for `Capella` products now exists in ESA SNAP. Add it in pre-processing.** - **ENH: Handling of Sentinel-1 [ASF](https://hyp3-docs.asf.alaska.edu/guides/rtc_product_guide/#readme-file) and [MPC](https://planetarycomputer.microsoft.com/dataset/sentinel-1-rtc) RTC products. ([#112](https://github.com/sertit/eoreader/issues/112), [#118](https://github.com/sertit/eoreader/issues/118))** - **ENH: Handling of Sentinel-1 SM products.** +- **ENH: Better handling of calibration step in SNAP for SAR data.** - FIX: Fix jpg, png... quicklooks management when plotting - FIX: Fix an `xarray` issue when trying to compute percentiles when stacking bands - DEPS: Remove as many mention as possible to `cloudpathlib` diff --git a/eoreader/data/grd_no_calib_preprocess_default.xml b/eoreader/data/grd_no_calib_preprocess_default.xml new file mode 100644 index 00000000..3a7e3be0 --- /dev/null +++ b/eoreader/data/grd_no_calib_preprocess_default.xml @@ -0,0 +1,80 @@ + + 1.0 + + Read + + + ${file} + + + + Terrain-Correction + + + + + + ${dem_name} + ${dem_path} + 0.0 + true + BILINEAR_INTERPOLATION + BILINEAR_INTERPOLATION + ${res_m} + ${res_deg} + ${crs} + false + 0.0 + 0.0 + false + false + false + false + false + false + true + false + false + false + false + Use projected local incidence angle from DEM + Use projected local incidence angle from DEM + Latest Auxiliary File + + + + + LinearToFromdB + + + + + + + + + Write + + + + + ${out} + BEAM-DIMAP + + + + + + + + + + + + + + + + + + diff --git a/eoreader/products/sar/cosmo_product.py b/eoreader/products/sar/cosmo_product.py index 38951a21..b6f51820 100644 --- a/eoreader/products/sar/cosmo_product.py +++ b/eoreader/products/sar/cosmo_product.py @@ -486,7 +486,7 @@ def _pre_process_sar(self, band, pixel_size: float = None, **kwargs) -> str: str: Band path """ with h5netcdf.File(self._img_path, phony_dims="access") as raw_h5: - if self.sar_prod_type == SarProductType.GDRG or self.nof_swaths == 1: + if self.nof_swaths == 1: return super()._pre_process_sar(band, pixel_size, **kwargs) else: LOGGER.warning( diff --git a/eoreader/products/sar/csg_product.py b/eoreader/products/sar/csg_product.py index 4fe31ebc..68d3763c 100644 --- a/eoreader/products/sar/csg_product.py +++ b/eoreader/products/sar/csg_product.py @@ -211,6 +211,20 @@ def _set_pixel_size(self) -> None: self.pixel_size = def_pixel_size self.resolution = def_res + def _post_init(self, **kwargs) -> None: + """ + Function used to post_init the products + (setting product-type, band names and so on) + """ + # Calibration fails with CSG data + LOGGER.debug( + "SNAP Error: Calibration currently fails for CSG data. Removing this step." + ) + self._calibrate = False + + # Post init done by the super class + super()._post_init(**kwargs) + def _set_instrument(self) -> None: """ Set instrument diff --git a/eoreader/products/sar/csk_product.py b/eoreader/products/sar/csk_product.py index 6725ce1f..b60adc0b 100644 --- a/eoreader/products/sar/csk_product.py +++ b/eoreader/products/sar/csk_product.py @@ -97,6 +97,21 @@ def _set_pixel_size(self) -> None: self.pixel_size = def_pixel_size self.resolution = def_res + def _post_init(self, **kwargs) -> None: + """ + Function used to post_init the products + (setting product-type, band names and so on) + """ + if self.nof_swaths > 1: + # Calibration fails with CSG data + LOGGER.debug( + "SNAP Error: Calibration currently fails for CSK data with multiple swaths. Removing this step." + ) + self._calibrate = False + + # Post init done by the super class + super()._post_init(**kwargs) + def _set_sensor_mode(self) -> None: """Get sensor mode""" # Get MTD XML file diff --git a/eoreader/products/sar/sar_product.py b/eoreader/products/sar/sar_product.py index 5e05a1ca..3af0631f 100644 --- a/eoreader/products/sar/sar_product.py +++ b/eoreader/products/sar/sar_product.py @@ -193,6 +193,9 @@ def __init__( self._snap_no_data = 0 self._raw_no_data = 0 + # Calibrate or not + self._calibrate = True + # Initialization from the super class super().__init__(product_path, archive_path, output_path, remove_tmp, **kwargs) @@ -703,41 +706,18 @@ def _pre_process_sar(self, band: sab, pixel_size: float = None, **kwargs) -> str # Pre-process graph if PP_GRAPH not in os.environ: - if ( - self.constellation == Constellation.CSG - and self.sar_prod_type == SarProductType.CPLX - ): - LOGGER.debug( - "SNAP Error: Calibration currently fails for CSG complex data. Removing this step." - ) - pp_graph = utils.get_data_dir().joinpath( - "cplx_no_calib_preprocess_default.xml" - ) - elif ( - self.constellation == Constellation.CSK and self.nof_swaths > 1 - ): - LOGGER.debug( - "SNAP Error: Calibration currently fails for CSK data with multiple swaths. Removing this step." - ) - pp_graph = utils.get_data_dir().joinpath( - "grd_sar_preprocess_fallback.xml" - ) - + if self.constellation_id == Constellation.S1.name: + sat = "s1" + if self.sensor_mode.value == "SM": + sat += "_sm" + elif not self._calibrate: + sat = "no_calib" else: - if self.constellation_id == Constellation.S1.name: - sat = "s1" - if self.sensor_mode.value == "SM": - sat += "_sm" - else: - sat = "sar" - spt = ( - "grd" - if self.sar_prod_type == SarProductType.GDRG - else "cplx" - ) - pp_graph = utils.get_data_dir().joinpath( - f"{spt}_{sat}_preprocess_default.xml" - ) + sat = "sar" + spt = "grd" if self.sar_prod_type == SarProductType.GDRG else "cplx" + pp_graph = utils.get_data_dir().joinpath( + f"{spt}_{sat}_preprocess_default.xml" + ) else: pp_graph = AnyPath(os.environ[PP_GRAPH]) if not pp_graph.is_file() or not pp_graph.suffix == ".xml": @@ -900,14 +880,27 @@ def interp_na(array, dim): return array - # Get .img file path (readable by rasterio) + # Get the .img path(s) try: - img = rasters.get_dim_img_path(dim_path, f"*{pol}*") + imgs = utils.get_dim_img_path(dim_path, f"*{pol}*") except FileNotFoundError: - img = rasters.get_dim_img_path(dim_path) # Maybe not the good name + imgs = utils.get_dim_img_path(dim_path) # Maybe not the good name + + # Manage cases where multiple swaths are ortho independently + if len(imgs) > 1: + mos_path, exists = self._get_out_path( + path.get_filename(dim_path) + "_mos.vrt" + ) + if not exists: + # Get .img file path (readable by rasterio) + + # Useful for PAZ SC data (multiswath) + rasters.merge_vrt(imgs, mos_path) + else: + mos_path = imgs[0] - # Open SAR image - with rioxarray.open_rasterio(str(img)) as arr: + # Open SAR image and convert it to a clean geotiff + with rioxarray.open_rasterio(mos_path) as arr: arr = arr.where(arr != self._snap_no_data, np.nan) # Interpolate if needed (interpolate na works only 1D-like, sadly) diff --git a/eoreader/products/sar/tsx_product.py b/eoreader/products/sar/tsx_product.py index 13f24f5c..7157b900 100644 --- a/eoreader/products/sar/tsx_product.py +++ b/eoreader/products/sar/tsx_product.py @@ -167,6 +167,24 @@ def __init__( if self.product_type != TsxProductType.SSC: self._geometric_res = getattr(TsxGeometricResolution, self.split_name[3]) + # Calibration (after init) + if ( + self.product_type == TsxProductType.SSC + and self.sensor_mode == TsxSensorMode.SC + ): + root, _ = self.read_mtd() + noise_corr_flag = root.findtext(".//noiseCorrectedFlag") + + if noise_corr_flag == "false": + # Don't calibrate SSC ScanSAR: + # PazCalibrator: org.esa.snap.core.gpf.OperatorException: Noise correction for ScanSAR is currently not supported. + # https://github.com/senbox-org/s1tbx/blob/1daae60d572e3ad0ee98c0ce3538c61ad71d7fa1/s1tbx-op-calibration/src/main/java/org/esa/s1tbx/calibration/gpf/calibrators/TerraSARXCalibrator.java#L215 + # https://github.com/senbox-org/s1tbx/blob/1daae60d572e3ad0ee98c0ce3538c61ad71d7fa1/s1tbx-op-calibration/src/main/java/org/esa/s1tbx/calibration/gpf/calibrators/PazCalibrator.java#L214 + LOGGER.warning( + "Noise correction for ScanSAR is currently not supported for PAZ data. Deactivating calibration." + ) + self._calibrate = False + def _set_pixel_size(self) -> None: """ Set product default pixel size (in meters) diff --git a/eoreader/utils.py b/eoreader/utils.py index 605b8113..14ae769f 100644 --- a/eoreader/utils.py +++ b/eoreader/utils.py @@ -30,7 +30,7 @@ from rasterio.enums import Resampling from rasterio.errors import NotGeoreferencedWarning from rasterio.rpc import RPC -from sertit import AnyPath, geometry, logs, rasters +from sertit import AnyPath, geometry, logs, path, rasters from sertit.types import AnyPathStrType, AnyPathType from eoreader import EOREADER_NAME @@ -456,3 +456,35 @@ def stack( stack = stack.rio.write_nodata(nodata, encoded=True, inplace=True) return stack, dtype + + +def get_dim_img_path(dim_path: AnyPathStrType, img_name: str = "*") -> list: + """ + Get the image path from a :code:`BEAM-DIMAP` data. + + A :code:`BEAM-DIMAP` file cannot be opened by rasterio, although its :code:`.img` file can. + + .. code-block:: python + + >>> dim_path = "path/to/dimap.dim" # BEAM-DIMAP image + >>> img_path = get_dim_img_path(dim_path) + + >>> # Read raster + >>> raster, meta = read(img_path) + + Args: + dim_path (AnyPathStrType): DIM path (.dim or .data) + img_name (str): .img file name (or regex), in case there are multiple .img files (ie. for S3 data) + + Returns: + list: .img files as a list + """ + dim_path = AnyPath(dim_path) + if dim_path.suffix == ".dim": + dim_path = dim_path.with_suffix(".data") + + assert dim_path.suffix == ".data" and dim_path.is_dir() + + return path.get_file_in_dir( + dim_path, img_name, extension="img", exact_name=True, get_list=True + )