Merge pull request #86 from sertit/0.19.4

0.19.4

CHANGES.md

@@ -1,5 +1,18 @@
 # Release History
 
+## 0.19.4 (2023-04-12)
+
+### Bug Fixes
+
+- FIX: Removing calibration step from SNAP pre-processing graph for multi-swath `Cosmo-SkyMed 1st GEN` products (to avoid ending up with empty images after pre-process)
+- FIX: Fixing the paths to Sentinel-2 quicklooks: using PVI instead of TCI file if no .jpg preview file is found ([#84](https://github.com/sertit/eoreader/issues/84)
+  , [#85](https://github.com/sertit/eoreader/issues/85), thanks a lot @floriandeboissieu)
+
+### Other
+
+- STAC: Updates in STAC management
+- INTERNAL: Use `geopandas.estimate_utm_crs()` when possible
+
 ## 0.19.3 (2023-03-24)
 
 ### Bug Fixes
@@ -19,7 +32,7 @@
 
 ### Bug Fixes
 
-- FIX: Fixing stack when saved as integer for some special cases 
+- FIX: Fixing stack when saved as integer for some special cases
 - FIX: Clipping negative reflectances to 0 ([#79](https://github.com/sertit/eoreader/issues/79))
 - FIX: Fixing nodata mangement for Theia product
 - FIX: Fixing handling of SCS multi-swath `Cosmo-SkyMed` products ([#78](https://github.com/sertit/eoreader/issues/78))

eoreader/__meta__.py

@@ -17,7 +17,7 @@
 """
 **EOReader** library
 """
-__version__ = "0.19.3"
+__version__ = "0.19.4"
 __title__ = "eoreader"
 __description__ = (
     "Remote-sensing opensource python library reading optical and SAR constellations, "

eoreader/products/custom_product.py

@@ -30,7 +30,7 @@
 from lxml.builder import E
 from rasterio import crs
 from rasterio.enums import Resampling
-from sertit import files, misc, rasters, vectors
+from sertit import files, misc, rasters
 from sertit.misc import ListEnum
 
 from eoreader import DATETIME_FMT, EOREADER_NAME, cache, utils
@@ -307,12 +307,9 @@ def crs(self) -> crs.CRS:
             extent_wgs84 = rasters.get_extent(self.get_default_band_path())
 
             # Get upper-left corner and deduce UTM proj from it
-            crs_str = vectors.corresponding_utm_projection(
-                extent_wgs84.bounds.minx, extent_wgs84.bounds.maxy
-            )
             raise InvalidProductError(
                 "Only stacks with projected CRS can be processed! "
-                f"Please reproject it to the corresponding UTM projection ({crs_str})!"
+                f"Please reproject it to the corresponding UTM projection ({extent_wgs84.estimate_utm_crs()})!"
             )
 
         return def_crs

eoreader/products/optical/s2_product.py

@@ -1603,20 +1603,14 @@ def get_quicklook_path(self) -> str:
                 else:
                     quicklook_path = str(next(self.path.glob("**/preview.jpg")))
             except (StopIteration, FileNotFoundError):
-                # Use the TCI
+                # Use the PVI
                 try:
-                    if self.product_type == S2ProductType.L2A:
-                        tci_regex = "TCI_60m"
-                    else:
-                        tci_regex = "TCI"
                     if self.is_archived:
                         quicklook_path = files.get_archived_rio_path(
-                            self.path, file_regex=rf".*{tci_regex}\.jp2"
+                            self.path, file_regex=r".*PVI\.jp2"
                         )
                     else:
-                        quicklook_path = str(
-                            next(self.path.glob(f"**/{tci_regex}.jp2"))
-                        )
+                        quicklook_path = str(next(self.path.glob("**/*PVI.jp2")))
                 except (StopIteration, FileNotFoundError):
                     LOGGER.warning(f"No quicklook found in {self.condensed_name}")
 

eoreader/products/sar/capella_product.py

@@ -273,11 +273,7 @@ def wgs84_extent(self) -> gpd.GeoDataFrame:
                     geometry=[Point(target_position)],
                     crs={"proj": "geocent", "ellps": "WGS84", "datum": "WGS84"},
                 ).to_crs(WGS84)
-
-                crs = vectors.corresponding_utm_projection(
-                    center_pix.x.iat[0], center_pix.x.iat[1]
-                )
-                center_pix.to_crs(crs)
+                center_pix.to_crs(center_pix.extent_wgs84.estimate_utm_crs())
 
                 # Get pixel spacing in meters
                 pixel_spacing_h = float(root.findtext(".//pixel_spacing_row"))

eoreader/products/sar/cosmo_product.py

@@ -116,17 +116,12 @@ def _pre_init(self, **kwargs) -> None:
         # SNAP cannot process its archive
         self.needs_extraction = True
 
-        # Post init done by the super class
-        super()._pre_init(**kwargs)
-
-    def _post_init(self, **kwargs) -> None:
-        """
-        Function used to post_init the products
-        (setting product-type, band names and so on)
-        """
+        # Get the number of swaths of this product
+        with h5netcdf.File(self._img_path, phony_dims="access") as raw_h5:
+            self.nof_swaths = len(list(raw_h5.groups))
 
         # Post init done by the super class
-        super()._post_init(**kwargs)
+        super()._pre_init(**kwargs)
 
     @cache
     def wgs84_extent(self) -> gpd.GeoDataFrame:
@@ -496,8 +491,7 @@ def _pre_process_sar(self, band, resolution: float = None, **kwargs) -> str:
             str: Band path
         """
         with h5netcdf.File(self._img_path, phony_dims="access") as raw_h5:
-            swaths = list(raw_h5.groups)
-            if self.sar_prod_type == SarProductType.GDRG or len(swaths) == 1:
+            if self.sar_prod_type == SarProductType.GDRG or self.nof_swaths == 1:
                 return super()._pre_process_sar(band, resolution, **kwargs)
             else:
                 LOGGER.warning(

eoreader/products/sar/sar_product.py

@@ -32,7 +32,7 @@
 from cloudpathlib import AnyPath, CloudPath
 from rasterio import crs
 from rasterio.enums import Resampling
-from sertit import files, misc, rasters, snap, strings, vectors
+from sertit import files, misc, rasters, snap, strings
 from sertit.misc import ListEnum
 
 from eoreader import EOREADER_NAME, cache, utils
@@ -179,6 +179,9 @@ def __init__(
         self.snap_filename = None
         """Path used by SNAP to process this product"""
 
+        self.nof_swaths = None
+        """Number of swaths of the current SAR product"""
+
         # Private attributes
         self._band_folder = None
         self._raw_band_regex = None
@@ -348,13 +351,8 @@ def extent(self) -> gpd.GeoDataFrame:
         # Get WGS84 extent
         extent_wgs84 = self.wgs84_extent()
 
-        # Get upper-left corner and deduce UTM proj from it
-        utm = vectors.corresponding_utm_projection(
-            extent_wgs84.bounds.minx, extent_wgs84.bounds.maxy
-        )
-        extent = extent_wgs84.to_crs(utm)
-
-        return extent
+        # Convert to UTM
+        return extent_wgs84.to_crs(self.crs())
 
     @cache
     def crs(self) -> crs.CRS:
@@ -375,12 +373,8 @@ def crs(self) -> crs.CRS:
         # Get WGS84 extent
         extent_wgs84 = self.wgs84_extent()
 
-        # Get upper-left corner and deduce UTM proj from it
-        crs_str = vectors.corresponding_utm_projection(
-            extent_wgs84.bounds.minx, extent_wgs84.bounds.maxy
-        )
-
-        return crs.CRS.from_string(crs_str)
+        # Estimate UTM from extent
+        return extent_wgs84.estimate_utm_crs()
 
     @abstractmethod
     def _set_sensor_mode(self) -> None:
@@ -695,11 +689,10 @@ def _pre_process_sar(self, band: sab, resolution: float = None, **kwargs) -> str
                             "cplx_no_calib_preprocess_default.xml"
                         )
                     elif (
-                        self.constellation == Constellation.CSK
-                        and self.sensor_mode.name == "HR"
+                        self.constellation == Constellation.CSK and self.nof_swaths > 1
                     ):
                         LOGGER.debug(
-                            "SNAP Error: Calibration currently fails for CSK HR data. Removing this step."
+                            "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"

eoreader/stac/__init__.py

@@ -74,6 +74,6 @@
     "ProjExt",
 ]
 
-from .stac_item import OPTICAL_STAC_EXTENSIONS, SAR_STAC_EXTENSIONS, StacItem
+from .stac_item import StacItem
 
-__all__ += ["StacItem", "OPTICAL_STAC_EXTENSIONS", "SAR_STAC_EXTENSIONS"]
+__all__ += ["StacItem"]

eoreader/stac/stac_extensions.py

@@ -64,6 +64,14 @@ class EoExt:
     """
 
     def __init__(self, prod, **kwargs):
+        try:
+            from pystac.extensions.eo import EOExtension
+
+            self.extension = EOExtension.get_schema_uri()
+        except ImportError:
+
+            self.extension = None
+
         self.cloud_cover = None
         self._prod = prod
 
@@ -123,9 +131,14 @@ def add_to_item(self, item) -> None:
             )
         # Add the EO extension
         eo_ext = EOExtension.ext(item, add_if_missing=True)
+
         if self.cloud_cover is not None:
             eo_ext.cloud_cover = self.cloud_cover
 
+        # TODO
+        # if self.snow_cover is not None:
+        #     eo_ext.snow_cover = self.snow_cover
+
         # Add band asset
         band_paths = self._prod.get_raw_band_paths()
         for band_name, band_path in band_paths.items():
@@ -181,6 +194,13 @@ class ProjExt:
     """
 
     def __init__(self, prod, **kwargs):
+        try:
+            from pystac.extensions.projection import ProjectionExtension
+
+            self.extension = ProjectionExtension.get_schema_uri()
+        except ImportError:
+            self.extension = None
+
         self._prod = prod
         self.epsg = self.crs().to_epsg()
         self.wkt2 = self.crs().to_wkt()
@@ -291,6 +311,12 @@ class ViewExt:
     """
 
     def __init__(self, prod, **kwargs):
+        try:
+            from pystac.extensions.view import ViewExtension
+
+            self.extension = ViewExtension.get_schema_uri()
+        except ImportError:
+            self.extension = None
 
         self.extension_fields = {
             "sun_az": VIEW_SUN_AZIMUTH,

eoreader/stac/stac_item.py

@@ -54,16 +54,6 @@
     repr_multiline_str,
 )
 
-SAR_STAC_EXTENSIONS = [
-    "https://stac-extensions.github.io/eo/v1.0.0/schema.json",
-    "https://stac-extensions.github.io/projection/v1.1.0/schema.json",
-]
-OPTICAL_STAC_EXTENSIONS = [
-    "https://stac-extensions.github.io/eo/v1.0.0/schema.json",
-    "https://stac-extensions.github.io/projection/v1.1.0/schema.json",
-    "https://stac-extensions.github.io/view/v1.0.0/schema.json",
-]
-
 
 class StacItem:
     """
@@ -92,11 +82,13 @@ def __init__(self, prod, **kwargs):
         self.datetime = self._prod.datetime
         self.geometry = gdf_to_geometry(self.geometry_fct())
         self.bbox = gdf_to_bbox(self.bbox_fct())
-        self.extensions = (
-            SAR_STAC_EXTENSIONS
-            if self._prod.sensor_type.value == "SAR"
-            else OPTICAL_STAC_EXTENSIONS
-        )
+
+        # Will be added if missing with xxxExtension.ext(item, add_if_missing=True)
+        self.extensions = [
+            ext.extension
+            for ext in [self.eo, self.proj, self.view]
+            if ext.extension is not None
+        ]
 
         # Common mtd
         self.gsd = self._prod.resolution
@@ -138,7 +130,6 @@ def create_item(self):
             geometry=self.geometry,
             bbox=self.bbox,
             properties=self.properties.copy(),
-            stac_extensions=self.extensions,
         )
 
         # Add assets

eoreader/stac/stac_utils.py

@@ -35,7 +35,7 @@ def fill_common_mtd(asset: Any, prod, **kwargs) -> None:
 
     Args:
         asset (pystac.Asset, pystac.Item): Asset or item
-        prod (product): EOReader product
+        prod (Product): EOReader product
         **kwargs: Additional arguments
     """
     # Basics
@@ -44,7 +44,7 @@ def fill_common_mtd(asset: Any, prod, **kwargs) -> None:
 
     # Date and Time
     asset.common_metadata.created = datetime.utcnow()
-    asset.common_metadata.updated = None  # TODO
+    asset.common_metadata.updated = datetime.utcnow()
 
     # Licensing
     # asset.common_metadata.license = None  # Collection level if possible
@@ -53,12 +53,14 @@ def fill_common_mtd(asset: Any, prod, **kwargs) -> None:
     # asset.common_metadata.providers = None  # Collection level if possible
 
     # Date and Time Range
-    asset.common_metadata.start_datetime = None  # TODO
-    asset.common_metadata.end_datetime = None  # TODO
+    asset.common_metadata.start_datetime = prod.datetime
+    asset.common_metadata.end_datetime = prod.datetime  # TODO
 
     # Instrument
     asset.common_metadata.platform = None  # TODO
-    asset.common_metadata.instruments = None  # TODO
+    asset.common_metadata.instruments = [
+        prod.instrument if isinstance(prod.instrument, str) else prod.instrument.value
+    ]
     asset.common_metadata.constellation = prod.constellation.value.lower()
     asset.common_metadata.mission = None
     asset.common_metadata.gsd = kwargs.get(GSD)

requirements.txt

@@ -34,7 +34,7 @@ rtree
 geopandas>=0.11.0
 
 # Spectral indices
-git+https://github.com/awesome-spectral-indices/spyndex.git
+spyndex>=0.3.0
 
 # STAC
 pystac[validation]