diff --git a/CHANGES.md b/CHANGES.md index d1978076..03acd824 100644 --- 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 index 4c6a8102..6ef50d6c 100644 --- 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 index 9e742b48..c5fe7258 100644 --- 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 index 15acee9f..f5a7421e 100644 --- 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 index 6a42333f..7bdd60e4 100644 --- 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