ENH: Better handling of calibration step in SNAP for SAR data

sertit · Jan 2, 2024 · 35b5c4f · 35b5c4f
1 parent 3d985d4
commit 35b5c4f
Show file tree

Hide file tree

Showing 8 changed files with 194 additions and 41 deletions.
diff --git 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
@@ -0,0 +1,80 @@
+<graph id="Graph">
+  <version>1.0</version>
+  <node id="Read">
+    <operator>Read</operator>
+    <sources/>
+    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
+      <file>${file}</file>
+    </parameters>
+  </node>
+  <node id="Terrain-Correction">
+    <operator>Terrain-Correction</operator>
+    <sources>
+      <sourceProduct refid="LinearToFromdB"/>
+    </sources>
+    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
+      <sourceBands/>
+      <demName>${dem_name}</demName>
+      <externalDEMFile>${dem_path}</externalDEMFile>
+      <externalDEMNoDataValue>0.0</externalDEMNoDataValue>
+      <externalDEMApplyEGM>true</externalDEMApplyEGM>
+      <demResamplingMethod>BILINEAR_INTERPOLATION</demResamplingMethod>
+      <imgResamplingMethod>BILINEAR_INTERPOLATION</imgResamplingMethod>
+      <pixelSpacingInMeter>${res_m}</pixelSpacingInMeter>
+      <pixelSpacingInDegree>${res_deg}</pixelSpacingInDegree>
+      <mapProjection>${crs}</mapProjection>
+      <alignToStandardGrid>false</alignToStandardGrid>
+      <standardGridOriginX>0.0</standardGridOriginX>
+      <standardGridOriginY>0.0</standardGridOriginY>
+      <nodataValueAtSea>false</nodataValueAtSea>
+      <saveDEM>false</saveDEM>
+      <saveLatLon>false</saveLatLon>
+      <saveIncidenceAngleFromEllipsoid>false</saveIncidenceAngleFromEllipsoid>
+      <saveLocalIncidenceAngle>false</saveLocalIncidenceAngle>
+      <saveProjectedLocalIncidenceAngle>false</saveProjectedLocalIncidenceAngle>
+      <saveSelectedSourceBand>true</saveSelectedSourceBand>
+      <applyRadiometricNormalization>false</applyRadiometricNormalization>
+      <saveSigmaNought>false</saveSigmaNought>
+      <saveGammaNought>false</saveGammaNought>
+      <saveBetaNought>false</saveBetaNought>
+      <incidenceAngleForSigma0>Use projected local incidence angle from DEM</incidenceAngleForSigma0>
+      <incidenceAngleForGamma0>Use projected local incidence angle from DEM</incidenceAngleForGamma0>
+      <auxFile>Latest Auxiliary File</auxFile>
+      <externalAuxFile/>
+    </parameters>
+  </node>
+  <node id="LinearToFromdB">
+    <operator>LinearToFromdB</operator>
+    <sources>
+      <sourceProduct refid="Read"/>
+    </sources>
+    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
+      <sourceBands/>
+    </parameters>
+  </node>
+  <node id="Write">
+    <operator>Write</operator>
+    <sources>
+      <sourceProduct refid="Terrain-Correction"/>
+    </sources>
+    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
+      <file>${out}</file>
+      <formatName>BEAM-DIMAP</formatName>
+    </parameters>
+  </node>
+  <applicationData id="Presentation">
+    <Description/>
+    <node id="Read">
+            <displayPosition x="37.0" y="35.0"/>
+    </node>
+    <node id="Terrain-Correction">
+      <displayPosition x="284.0" y="259.0"/>
+    </node>
+    <node id="LinearToFromdB">
+      <displayPosition x="138.0" y="263.0"/>
+    </node>
+    <node id="Write">
+            <displayPosition x="522.0" y="264.0"/>
+    </node>
+  </applicationData>
+</graph>
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
+    )