ENH: Adding the support of Pleides Neo SEN and PRJ products

bastiencyr · Sep 13, 2022 · 22464f4 · 22464f4
1 parent 2dacdc3
commit 22464f4
Show file tree

Hide file tree

Showing 5 changed files with 58 additions and 55 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -6,6 +6,7 @@
 
 - **ENH: Adding the support of RapidEye constellation**
 - **ENH: Adding the support of Landsat Level-2 products** ([#49](https://github.com/sertit/eoreader/issues/49))
+- **ENH: Adding the support of Pleides Neo SEN and PRJ products** *(needs GDAL 3.5+ or rasterio 1.3.0+)*
 - **ENH: Adding the function `bands.is_thermal_band`**
 
 ### Bug Fixes

diff --git a/docs/optical.md b/docs/optical.md
@@ -24,7 +24,7 @@ The Data Access Portfolio Document presents the offer of the datasets and data a
 | PlanetScope        | {meth}`~eoreader.products.PlaProduct`                   | L3A & L3B                 | 3m                                   | ✅                 |
 | SkySat             | {meth}`~eoreader.products.SkyProduct`                   | (Collect Product) ortho_* | 0.5m                                 | ✅                 |
 | RapidEye           | {meth}`~eoreader.products.ReProduct`                    | L3A                       | 5m                                   | ✅                 |
-| Pleiades-Neo       | {meth}`~eoreader.products.PneoProduct`                  | ORT                       | 0.3 (PAN), 1.2m (MS)                 | ✅                 |
+| Pleiades-Neo       | {meth}`~eoreader.products.PneoProduct`                  | SEN, PRJ, ORT & MOS                       | 0.3 (PAN), 1.2m (MS)                 | ✅                 |
 | Pleiades           | {meth}`~eoreader.products.PldProduct`                   | SEN, PRJ, ORT & MOS       | 0.5 (PAN), 2m (MS)                   | ✅                 |
 | Vision-1           | {meth}`~eoreader.products.Vis1Product`                  | PRJ & ORTP                | 0.9 (PAN), 3.5m (MS)                 | ✅                 |
 | SPOT 4             | {meth}`~eoreader.products.Spot45Product`                | L1A, M1B, L2A             | 10 (PAN), 20m (MS)                   | ✅                 |

diff --git a/eoreader/products/optical/dimap_product.py b/eoreader/products/optical/dimap_product.py
@@ -38,6 +38,7 @@
 from rasterio import features, transform
 from sertit import files, rasters_rio, vectors
 from sertit.misc import ListEnum
+from sertit.vectors import WGS84
 from shapely.geometry import Polygon, box
 
 from eoreader import cache, utils
@@ -995,6 +996,11 @@ def open_mask(self, mask_str: str, **kwargs) -> gpd.GeoDataFrame:
                 DimapProductType.SEN,
                 DimapProductType.PRJ,
             ]:
+                # Sometimes the GML mask lacks crs (why ?)
+                if not mask.crs:
+                    mask.crs = self._get_raw_crs()
+
+                mask.crs = WGS84
                 LOGGER.info(f"Orthorectifying {mask_str}")
                 with rasterio.open(str(self._get_tile_path())) as dim_dst:
                     # Rasterize mask (no transform as we have the vector in image geometry)
@@ -1006,24 +1012,31 @@ def open_mask(self, mask_str: str, **kwargs) -> gpd.GeoDataFrame:
                         default_value=self._mask_true,  # Inside vector
                         dtype=np.uint8,
                     )
+                    # Check mask validity (to avoid reprojecting)
+                    # All null
+                    if mask_raster.max() == 0:
+                        mask = gpd.GeoDataFrame(geometry=[], crs=crs)
+                    # All valid
+                    elif mask_raster.min() == 1:
+                        pass
+                    else:
+                        # Reproject mask raster
+                        LOGGER.debug(f"\tReprojecting {mask_str}")
+                        dem_path = self._get_dem_path(**kwargs)
+                        reproj_data = self._reproject(
+                            mask_raster, dim_dst.meta, dim_dst.rpcs, dem_path, **kwargs
+                        )
 
-                    # Reproject mask raster
-                    LOGGER.debug(f"\tReprojecting {mask_str}")
-                    dem_path = self._get_dem_path(**kwargs)
-                    reproj_data = self._reproject(
-                        mask_raster, dim_dst.meta, dim_dst.rpcs, dem_path, **kwargs
-                    )
-
-                    # Vectorize mask raster
-                    LOGGER.debug(f"\tRevectorizing {mask_str}")
-                    mask = rasters_rio.vectorize(
-                        reproj_data,
-                        values=self._mask_true,
-                        default_nodata=self._mask_false,
-                    )
+                        # Vectorize mask raster
+                        LOGGER.debug(f"\tRevectorizing {mask_str}")
+                        mask = rasters_rio.vectorize(
+                            reproj_data,
+                            values=self._mask_true,
+                            default_nodata=self._mask_false,
+                        )
 
-                    # Do not keep pixelized mask
-                    mask = utils.simplify_footprint(mask, self.resolution)
+                        # Do not keep pixelized mask
+                        mask = utils.simplify_footprint(mask, self.resolution)
 
             # Sometimes the GML mask lacks crs (why ?)
             elif (
@@ -1036,8 +1049,7 @@ def open_mask(self, mask_str: str, **kwargs) -> gpd.GeoDataFrame:
                 ]
             ):
                 # Convert to target CRS
-                mask.crs = self._get_raw_crs()
-                mask = mask.to_crs(self.crs())
+                mask.crs = self.crs()
 
             # Save to file
             if mask.empty:
@@ -1066,15 +1078,19 @@ def _load_nodata(
         """
         nodata_det = self.open_mask("ROI", **kwargs)
 
-        # Rasterize nodata
-        return features.rasterize(
-            nodata_det.geometry,
-            out_shape=(height, width),
-            fill=self._mask_true,  # Outside ROI = nodata (inverted compared to the usual)
-            default_value=self._mask_false,  # Inside ROI = not nodata
-            transform=transform,
-            dtype=np.uint8,
-        )
+        if all(nodata_det.is_empty):
+            return np.zeros((1, height, width), dtype=np.uint8)
+        else:
+
+            # Rasterize nodata
+            return features.rasterize(
+                nodata_det.geometry,
+                out_shape=(height, width),
+                fill=self._mask_true,  # Outside ROI = nodata (inverted compared to the usual)
+                default_value=self._mask_false,  # Inside ROI = not nodata
+                transform=transform,
+                dtype=np.uint8,
+            )
 
     def _get_tile_path(self) -> Union[CloudPath, Path]:
         """

diff --git a/eoreader/products/optical/pneo_product.py b/eoreader/products/optical/pneo_product.py
@@ -20,10 +20,6 @@
 for more information.
 """
 import logging
-from pathlib import Path
-from typing import Union
-
-from cloudpathlib import CloudPath
 
 from eoreader.bands import SpectralBand
 from eoreader.bands import spectral_bands as spb
@@ -53,26 +49,6 @@ def _pre_init(self, **kwargs) -> None:
         # Post init done by the super class
         super()._pre_init(**kwargs)
 
-    def _get_ortho_path(self, **kwargs) -> Union[CloudPath, Path]:
-        """
-        Get the orthorectified path of the bands.
-
-        Returns:
-            Union[CloudPath, Path]: Orthorectified path
-        """
-        if self.product_type in self._proj_prod_type:
-            # ERROR BECAUSE RASTERIO/GDAL DOES NOT HANDLE PLEIADES-NEO RPCs
-            raise NotImplementedError(
-                "Pleiades-Neo RPCs file nomenclature is not yet handled by rasterio. "
-                "See https://github.com/rasterio/rasterio/issues/2388. "
-                "GDAL PR is here: https://github.com/OSGeo/gdal/pull/5090"
-            )
-
-        else:
-            ortho_path = self._get_tile_path()
-
-        return ortho_path
-
     def _map_bands(self) -> None:
         """Set products type"""
         # Create spectral bands

diff --git a/eoreader/products/optical/vhr_product.py b/eoreader/products/optical/vhr_product.py
@@ -31,7 +31,7 @@
 import rasterio
 import xarray as xr
 from cloudpathlib import AnyPath, CloudPath
-from rasterio import rpc, warp
+from rasterio import MemoryFile, rpc, warp
 from rasterio.crs import CRS
 from rasterio.enums import Resampling
 from sertit import files, rasters, rasters_rio
@@ -279,6 +279,7 @@ def _reproject(
             "RPC_DEM": dem_path,
             "RPC_DEM_MISSING_VALUE": 0,
             "OSR_USE_ETMERC": "YES",
+            "BIGTIFF": "IF_NEEDED",
         }
 
         # Reproject
@@ -312,9 +313,18 @@ def _reproject(
         # Just in case, read the array with the most appropriate resolution
         # as the warping sometimes gives not the closest resolution possible to the wanted one
         if not math.isclose(dst_transform.a, self.resolution, rel_tol=1e-4):
-            out_arr, meta = rasters_rio.read(
-                (out_arr, meta), resolution=self.resolution
+            LOGGER.debug(
+                f"Reproject resolution ({dst_transform.a}) is too far from the wanted one ({self.resolution}). "
+                "Resampling the stack."
+                ""
             )
+            with MemoryFile() as memfile:
+                with memfile.open(
+                    **meta, BIGTIFF=rasters_rio.bigtiff_value(out_arr)
+                ) as dst:
+                    dst.write(out_arr)
+                    out_arr = None  # free memory
+                    out_arr, meta = rasters_rio.read(dst, resolution=self.resolution)
 
         return out_arr, meta