ENH: Handling COSMO product with only the h5 file in it (if missi…

…ng XML metadata file)** sertit#36
bastiencyr · May 27, 2022 · a48c2b5 · a48c2b5
1 parent 0a927de
commit a48c2b5
Show file tree

Hide file tree

Showing 5 changed files with 215 additions and 5 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -32,6 +32,7 @@
 - **ENH: Extending `get_raw_band_paths` to every product ([#31](https://github.com/sertit/eoreader/issues/31))**
 - **ENH: Adding a `is_ortho` attribute corresponding to when the product is already orthorectified/geocoded, in order to avoid computing heavy processes without wanting it (i.e. footprint...)**
 - **ENH: Adding the instrument name of every constellation, under `prod.instrument`**
+- **ENH: Handling `COSMO` product with only the `h5` file in it (if missing XML metadata file)** ([#36](https://github.com/sertit/eoreader/issues/36))
 
 ### Optimizations
 

diff --git a/docs/faq.md b/docs/faq.md
@@ -69,9 +69,22 @@ os.environ["PATH"] += r";C:\Program Files\snap\bin"
 
 ### SNAP known bugs
 
+#### SNAP `secure-processing` not recognized
 Sometimes SNAP process returns `Feature 'http://javax.xml.XMLConstants/feature/secure-processing' is not recognized.`
 
 This is a known SNAP bug.  
 
 Just add the line `-Djavax.xml.parsers.SAXParserFactory=com.sun.org.apache.xerces.internal.jaxp.SAXParserFactoryImpl` to your `gpt.vmoptions` file.  
-Please look at [this issue](https://forum.step.esa.int/t/xmlfactory-error-using-snap-8/26566) for more information.
+Please look at [this issue](https://forum.step.esa.int/t/xmlfactory-error-using-snap-8/26566) for more information.
+
+#### COSMO-SkyMed orthorectified files are empty
+
+For an unknown reason, the SNAP calibration step doesn't work and set nodata everywhere.
+A workaround is to set the file `cplx_no_calib_preprocess_default.xml` stored in `eoreader/data` in the `EOREADER_PP_GRAPH` environment variable.
+Even if the product won't be calibrated, you will be able to work with some orthorectified data.
+
+```python
+import os
+os.env["EOREADER_PP_GRAPH"] = "/home/eoreader/data/cplx_no_calib_preprocess_default.xml"
+prod.load(VV)
+```
diff --git a/eoreader/data/cplx_no_calib_preprocess_default.xml b/eoreader/data/cplx_no_calib_preprocess_default.xml
@@ -0,0 +1,106 @@
+<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="Multilook">
+    <operator>Multilook</operator>
+    <sources>
+      <sourceProduct refid="Read"/>
+    </sources>
+    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
+      <sourceBands/>
+      <nRgLooks>2</nRgLooks>
+      <nAzLooks>2</nAzLooks>
+      <outputIntensity>true</outputIntensity>
+      <grSquarePixel>true</grSquarePixel>
+    </parameters>
+  </node>
+
+
+  <node id="Terrain-Correction">
+    <operator>Terrain-Correction</operator>
+    <sources>
+      <sourceProduct refid="Multilook"/>
+    </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="Terrain-Correction"/>
+    </sources>
+    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
+      <sourceBands/>
+    </parameters>
+  </node>
+
+  <node id="Write">
+    <operator>Write</operator>
+    <sources>
+      <sourceProduct refid="LinearToFromdB"/>
+    </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="7.0" y="131.0"/>
+    </node>
+    <node id="Calibration">
+      <displayPosition x="43.0" y="188.0"/>
+    </node>
+    <node id="Terrain-Correction">
+      <displayPosition x="143.0" y="265.0"/>
+    </node>
+    <node id="LinearToFromdB">
+      <displayPosition x="256.0" y="299.0"/>
+    </node>
+    <node id="Multilook">
+      <displayPosition x="108.0" y="224.0"/>
+    </node>
+    <node id="Write">
+      <displayPosition x="461.0" y="307.0"/>
+    </node>
+  </applicationData>
+</graph>
diff --git a/eoreader/products/sar/cosmo_product.py b/eoreader/products/sar/cosmo_product.py
@@ -26,11 +26,13 @@
 from typing import Union
 
 import geopandas as gpd
+import h5netcdf
 import numpy as np
 import rasterio
 import xarray as xr
 from cloudpathlib import AnyPath, CloudPath
 from lxml import etree
+from lxml.builder import E
 from sertit import files, rasters, strings, vectors
 from sertit.misc import ListEnum
 from shapely.geometry import Polygon, box
@@ -342,9 +344,97 @@ def _read_mtd(self) -> (etree._Element, dict):
         Returns:
             (etree._Element, dict): Metadata XML root and its namespaces
         """
-        mtd_from_path = "DFDN_*.h5.xml"
+        try:
+            mtd_from_path = "DFDN_*.h5.xml"
 
-        return self._read_mtd_xml(mtd_from_path)
+            return self._read_mtd_xml(mtd_from_path)
+        except InvalidProductError:
+            try:
+                field_map = {
+                    # ProductInfo
+                    "ProductName": "Product Filename",
+                    # "ProductId": ,
+                    "MissionId": "Mission ID",
+                    # "UniqueIdentifier": ,
+                    "ProductGenerationDate": "Product Generation UTC",
+                    # "UserRequestId": ,
+                    # "ServiceRequestName": ,
+                    # ProductDefinitionData
+                    "ProductType": "Product Type",
+                    "SceneSensingStartUTC": "Scene Sensing Start UTC",
+                    "SceneSensingStopUTC": "Scene Sensing Stop UTC",
+                    # "GeoCoordTopRightEN": ,
+                    "GeoCoordSceneCentre": "Scene Centre Geodetic Coordinates",
+                    "SatelliteId": "Satellite ID",
+                    "AcquisitionMode": "Acquisition Mode",
+                    "LookSide": "Look Side",
+                    "ProjectionId": "Projection ID",
+                    "DeliveryMode": "Delivery Mode",
+                    "AcquisitionStationId": "Acquisition Station ID",
+                    # ProcessingInfo
+                    # "ProcessingLevel":,
+                    # ProductCharacteristics
+                    "AzimuthGeometricResolution": "Azimuth Geometric Resolution",
+                    "GroundRangeGeometricResolution": "Ground Range Geometric Resolution",
+                }
+
+                sbi_field_map = {
+                    "GeoCoordBottomLeft": "Bottom Left Geodetic Coordinates",
+                    "GeoCoordBottomRight": "Bottom Right Geodetic Coordinates",
+                    "GeoCoordTopLeft": "Top Left Geodetic Coordinates",
+                    "GeoCoordTopRight": "Top Right Geodetic Coordinates",
+                    "GeoCoordTopRightEN": "Top Right East-North",
+                    "NearLookAngle": "Near Look Angle",
+                    "FarLookAngle": "Far Look Angle",
+                }
+
+                def h5_to_str(h5_val):
+                    str_val = str(h5_val)
+                    str_val = str_val.replace("[", "")
+                    str_val = str_val.replace("]", "")
+                    return str_val
+
+                with h5netcdf.File(self._img_path) as netcdf_ds:
+
+                    # Create XML attributes
+                    global_attr = []
+                    for xml_attr, h5_attr in field_map.items():
+                        try:
+                            global_attr.append(
+                                E(xml_attr, h5_to_str(netcdf_ds.attrs[h5_attr]))
+                            )
+                        except KeyError:
+                            # CSG products don't have their ProductName in the h5 file...
+                            if xml_attr == "ProductName":
+                                global_attr.append(
+                                    E(xml_attr, files.get_filename(self._img_path))
+                                )
+
+                    try:
+                        # CSK products
+                        sbi = netcdf_ds.groups["S01"].variables["SBI"]
+                    except KeyError:
+                        # CSG products
+                        sbi = netcdf_ds.groups["S01"].variables["IMG"]
+
+                    for xml_attr, h5_attr in sbi_field_map.items():
+                        global_attr.append(E(xml_attr, h5_to_str(sbi.attrs[h5_attr])))
+
+                    mtd = E.s3_global_attributes(*global_attr)
+                    mtd_el = etree.fromstring(
+                        etree.tostring(
+                            mtd,
+                            pretty_print=True,
+                            xml_declaration=True,
+                            encoding="UTF-8",
+                        )
+                    )
+
+                return mtd_el, {}
+            except KeyError:
+                raise InvalidProductError(
+                    "Missing the XML metadata file. Cannot read the product."
+                )
 
     def get_quicklook_path(self) -> str:
         """

diff --git a/eoreader/reader.py b/eoreader/reader.py
@@ -248,8 +248,8 @@ class Constellation(ListEnum):
         "nested": -1,  # File that can be found at any level (product/**/file)
         "regex": r"\d{8}_\d{6}_ssc\w{1,4}_\w{4,5}_.*metadata.*\.json",
     },
-    Constellation.CSK: rf"{CONSTELLATION_REGEX[Constellation.CSK][1]}\.xml",
-    Constellation.CSG: rf"{CONSTELLATION_REGEX[Constellation.CSG][1]}\.xml",
+    Constellation.CSK: rf"{CONSTELLATION_REGEX[Constellation.CSK][1]}",
+    Constellation.CSG: rf"{CONSTELLATION_REGEX[Constellation.CSG][1]}",
     Constellation.TSX: rf"{CONSTELLATION_REGEX[Constellation.TSX]}\.xml",
     Constellation.TDX: rf"{CONSTELLATION_REGEX[Constellation.TDX]}\.xml",
     Constellation.PAZ: rf"{CONSTELLATION_REGEX[Constellation.PAZ]}\.xml",