update rio-tiler version and add Copernicus DEM

.pre-commit-config.yaml

@@ -31,3 +31,4 @@ repos:
         exclude: tests/.*
         additional_dependencies:
         - types-attrs
+        - types-cachetools

CHANGES.md

@@ -1,5 +1,11 @@
 # Release Notes
 
+## 0.9.0 (2023-07-13)
+
+* update rio-tiler requirement to `>=5.0,<6.0`
+* add `rio_tiler_pds.copernicus.aws.Dem30Reader` and `rio_tiler_pds.copernicus.aws.Dem90Reader` **mosaic** readers
+* add `boto3` in dependencies
+
 ## 0.8.0 (2023-04-11)
 
 * remove Landsat 8 Collection 1

README.md

@@ -62,6 +62,7 @@ $ pip install git+https://github.com/cogeotiff/rio-tiler-pds.git
 | [CBERS 4/4A][cbers_cog]                   | L2/L4                                       | COG                        | AMS Kepler / AWS           | us-east-1    | **Requester-pays** |
 | [MODIS (modis-pds)][modis_pds]            | MCD43A4, MOD09GQ, MYD09GQ, MOD09GA, MYD09GA | GTiff (External Overviews) | -                          | us-west-2    | Public             |
 | [MODIS (astraea-opendata)][modis_astraea] | MCD43A4, MOD11A1, MOD13A1, MYD11A1 MYD13A1  | COG                        | Astraea / AWS              | us-west-2    | **Requester-pays** |
+| [Copernicus Digital Elevation Model][copernicus_dem] | GLO-30, GLO-90                   | COG                        | Sinergise / AWS            | eu-central-1    | Public |
 
 [s2_l1c_jp2]: https://registry.opendata.aws/sentinel-2/
 [s2_l2a_jp2]: https://registry.opendata.aws/sentinel-2/
@@ -71,6 +72,7 @@ $ pip install git+https://github.com/cogeotiff/rio-tiler-pds.git
 [cbers_cog]: https://registry.opendata.aws/cbers/
 [modis_pds]: https://docs.opendata.aws/modis-pds/readme.html
 [modis_astraea]: https://registry.opendata.aws/modis-astraea/
+[copernicus_dem]: https://registry.opendata.aws/copernicus-dem/
 
 **Adding more dataset**:
 
@@ -119,6 +121,7 @@ from rio_tiler_pds.sentinel.aws import (
 
 from rio_tiler_pds.cbers.aws import CBERSReader
 from rio_tiler_pds.modis.aws import MODISPDSReader, MODISASTRAEAReader
+from rio_tiler_pds.copernicus.aws import Dem30Reader, Dem90Reader
 ```
 
 All Readers are subclass of [`rio_tiler.io.BaseReader`](https://github.com/cogeotiff/rio-tiler/blob/f917d0eaf27f8644f3bb18856a63fe45eeb4a2ef/rio_tiler/io/base.py#L17) and inherit its properties/methods.
@@ -324,6 +327,19 @@ with S2COGReader("S2A_L2A_20170729_19UDP_0") as sentinel:
       >> 2.0
 ```
 
+### Mosaic Reader: Copernicus DEM
+
+The Copernicus DEM GLO-30 and GLO-90 readers are not **per scene** but **mosaic** readers. This is possible because the dataset is a global dataset with file names having the `geo-location` of the COG, meaning we can easily contruct a filepath from a coordinate.
+
+```python
+from rio_tiler_pds.copernicus.aws import Dem30Reader
+
+with Dem30Reader() as dem:
+    print(dem.assets_for_point(-57.2, -11.2))
+
+>> ['s3://copernicus-dem-30m/Copernicus_DSM_COG_10_S12_00_W058_00_DEM/Copernicus_DSM_COG_10_S12_00_W058_00_DEM.tif']
+```
+
 ## Changes
 
 See [CHANGES.md](https://github.com/cogeotiff/rio-tiler-pds/blob/main/CHANGES.md).

docs/mkdocs.yml

@@ -28,7 +28,6 @@ nav:
     - Overview: usage/overview.md
     - Sentinel: usage/sentinel.md
     - Landsat Collection 2: usage/landsat-c2.md
-    - Landsat: usage/landsat.md
     - CBERS: usage/cbers.md
     - MODIS: usage/modis.md
   - Development - Contributing: 'contributing.md'

docs/src/usage/overview.md

@@ -14,5 +14,5 @@ from rio_tiler_pds.sentinel.aws import S2JP2Reader
 
 with rasterio.Env(AWS_REQUEST_PAYER="requester"):
     with S2JP2Reader("S2A_L1C_20170729_19UDP_0") as s2:
-        print(s2.preview(bands="B01", width=64, height=64))
+        print(s2.preview(bands="B01", width=64, height=64, max_size=None))
 ```

docs/src/usage/sentinel.md

@@ -81,19 +81,34 @@ with rasterio.Env(AWS_REQUEST_PAYER="requester"):
         assert img.data.shape == (1, 256, 256)
 
         print(sentinel.point(-69.41, 48.25, bands=("B01", "B02")))
-        # Result is in form of
-        # [
-        #   value for band 1 in band B01,
-        #   value for band 1 in band B02
-        # ]
-        > [1230, 875]
+        >> PointData(
+            array=masked_array(data=[1201, 843], mask=[False, False], fill_value=999999, dtype=uint16),
+            band_names=['B01', 'B02'],
+            coordinates=(-69.41, 48.25),
+            crs=CRS.from_epsg(4326),
+            assets=[
+                's3://sentinel-s2-l1c/tiles/19/U/DP/2017/7/29/0/B01.jp2',
+                's3://sentinel-s2-l1c/tiles/19/U/DP/2017/7/29/0/B02.jp2'
+            ],
+            metadata={}
+        )
 
         # Working with Expression
         img = sentinel.tile(77, 89, 8, expression="B01/B02")
         assert igm.data.shape == (1, 256, 256)
 
         print(sentinel.point(-69.41, 48.25, expression="B01/B02"))
-        > [1.424673784104389]
+        >> PointData(
+            array=masked_array(data=[1.424673784104389], mask=[False], fill_value=999999, dtype=float32),
+            band_names=['B01/B02'],
+            coordinates=(-69.41, 48.25),
+            crs=CRS.from_epsg(4326),
+            assets=[
+                's3://sentinel-s2-l1c/tiles/19/U/DP/2017/7/29/0/B01.jp2',
+                's3://sentinel-s2-l1c/tiles/19/U/DP/2017/7/29/0/B02.jp2'
+            ],
+            metadata={}
+        )
 ```
 
 ### L2A - JPEG2000

pyproject.toml

@@ -20,7 +20,7 @@ classifiers = [
     "Typing :: Typed",
 ]
 dynamic = ["version"]
-dependencies = ["rio-tiler>=4.0,<5.0"]
+dependencies = ["rio-tiler>=5.0,<6.0", "boto3"]
 
 [project.optional-dependencies]
 test = ["pytest", "pytest-cov"]

rio_tiler_pds/copernicus/__init__.py

@@ -0,0 +1,3 @@
+"""rio-tiler-pds.copernicus"""
+
+from rio_tiler_pds.copernicus import aws  # noqa

rio_tiler_pds/copernicus/aws/__init__.py

@@ -0,0 +1,3 @@
+"""AWS copernicus dem readers"""
+
+from .dem import Dem30Reader, Dem90Reader  # noqa

rio_tiler_pds/copernicus/aws/dem.py

@@ -0,0 +1,244 @@
+"""Copernicus DEM Mosaic Reader.
+
+https://registry.opendata.aws/copernicus-dem/
+"""
+
+import json
+import math
+from typing import Any, Dict, List, Optional, Tuple, Type
+
+import attr
+from morecantile import TileMatrixSet
+from rasterio.crs import CRS
+from rasterio.errors import RasterioIOError
+from rasterio.warp import transform as transform_coords
+from rasterio.warp import transform_bounds
+
+from rio_tiler.constants import MAX_THREADS, WEB_MERCATOR_TMS, WGS84_CRS
+from rio_tiler.errors import PointOutsideBounds, TileOutsideBounds
+from rio_tiler.io import BaseReader, Reader
+from rio_tiler.models import BandStatistics, ImageData, Info, PointData
+from rio_tiler.mosaic import mosaic_point_reader, mosaic_reader
+from rio_tiler.types import BBox
+
+
+@attr.s
+class Dem30Reader(BaseReader):
+    """Simple Mosaic reader for copernicus DEM"""
+
+    input: str = attr.ib(default="copernicus-dem-30m")
+
+    tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
+
+    minzoom: int = attr.ib(default=7)
+    maxzoom: int = attr.ib(default=8)
+
+    bounds: BBox = attr.ib(init=False, default=(-180, -90, 180, 90))
+    crs: CRS = attr.ib(init=False, default=WGS84_CRS)
+
+    reader: Type[Reader] = attr.ib(default=Reader)
+    reader_options: Dict = attr.ib(factory=dict)
+
+    colormap: Dict = attr.ib(init=False, factory=dict)
+
+    _scheme: str = "s3"
+    bucket: str = attr.ib(default="copernicus-dem-30m")
+    prefix_pattern: str = attr.ib(
+        default="Copernicus_DSM_COG_10_{northsouth}{lat:02d}_00_{eastwest}{lon:03d}_00_DEM/Copernicus_DSM_COG_10_{northsouth}{lat:02d}_00_{eastwest}{lon:03d}_00_DEM"
+    )
+
+    def tile(
+        self,
+        tile_x: int,
+        tile_y: int,
+        tile_z: int,
+        reverse: bool = False,
+        threads: int = MAX_THREADS,
+        **kwargs: Any,
+    ) -> ImageData:
+        """Get Tile."""
+        assets = self.assets_for_tile(tile_x, tile_y, tile_z)
+        mosaic_assets = list(reversed(assets)) if reverse else assets
+
+        def _reader(
+            asset: str, tile_x: int, tile_y: int, tile_z: int, **kwargs: Any
+        ) -> ImageData:
+            with Reader(asset, tms=self.tms) as src:
+                return src.tile(tile_x, tile_y, tile_z, **kwargs)
+
+        return mosaic_reader(
+            mosaic_assets,
+            _reader,
+            tile_x,
+            tile_y,
+            tile_z,
+            threads=threads,
+            allowed_exceptions=(TileOutsideBounds, RasterioIOError),
+            **kwargs,
+        )[0]
+
+    def point(
+        self,
+        lon: float,
+        lat: float,
+        coord_crs: CRS = WGS84_CRS,
+        reverse: bool = False,
+        threads: int = MAX_THREADS,
+        **kwargs: Any,
+    ) -> PointData:
+        """Get Point value."""
+        assets = self.assets_for_point(lon, lat, coord_crs=coord_crs)
+        mosaic_assets = list(reversed(assets)) if reverse else assets
+
+        def _reader(asset: str, lon: float, lat: float, **kwargs) -> PointData:
+            with Reader(asset, tms=self.tms) as src:
+                return src.point(lon, lat, **kwargs)
+
+        return mosaic_point_reader(
+            mosaic_assets,
+            _reader,
+            lon,
+            lat,
+            threads=threads,
+            allowed_exceptions=(PointOutsideBounds, RasterioIOError),
+            **kwargs,
+        )[0]
+
+    def info(self) -> Info:
+        """info."""
+        meta = {
+            "bounds": self.geographic_bounds,
+            "minzoom": self.minzoom,
+            "maxzoom": self.maxzoom,
+            "band_metadata": [("b1", {})],
+            "band_descriptions": [("b1", "")],
+            "dtype": "float32",
+            "colorinterp": ["grey"],
+            "nodata_type": "None",
+            "driver": "GTiff",
+            "count": 1,
+        }
+        return Info(**meta)
+
+    def statistics(self, **kwargs: Any) -> Dict[str, BandStatistics]:
+        """Return Dataset's statistics."""
+        # FOR NOW WE ONLY RETURN VALUE FROM THE FIRST FILE
+        dataset = self.prefix_pattern.format(northsouth="N", eastwest="E", lat=0, lon=6)
+        with Reader(f"{self._scheme}://{self.bucket}/{dataset}.tif") as src:
+            return src.statistics(**kwargs)
+
+    ############################################################################
+    # Not Implemented methods
+    # BaseReader required those method to be implemented
+    def preview(self, *args, **kwargs):
+        """Placeholder for BaseReader.preview."""
+        raise NotImplementedError
+
+    def part(self, *args, **kwargs):
+        """Placeholder for BaseReader.part."""
+        raise NotImplementedError
+
+    def feature(self, *args, **kwargs):
+        """Placeholder for BaseReader.feature."""
+        raise NotImplementedError
+
+    def assets_for_tile(self, x: int, y: int, z: int) -> List[str]:
+        """Retrieve assets for tile."""
+        xmin, ymin, xmax, ymax = self.tms.bounds(x, y, z)
+        return self.assets_for_bbox(
+            xmin, ymin, xmax, ymax, coord_crs=self.tms.rasterio_geographic_crs
+        )
+
+    def assets_for_bbox(
+        self,
+        xmin: float,
+        ymin: float,
+        xmax: float,
+        ymax: float,
+        coord_crs: Optional[CRS] = WGS84_CRS,
+    ) -> List[str]:
+        """Retrieve assets for bbox."""
+        if coord_crs != WGS84_CRS:
+            xmin, ymin, xmax, ymax = transform_bounds(
+                coord_crs,
+                WGS84_CRS,
+                xmin,
+                ymin,
+                xmax,
+                ymax,
+            )
+
+        xmin = int(xmin + 360)
+        xmax = int(xmax + 360)
+        ymax = int(ymax + 180)
+        ymin = int(ymin + 180)
+
+        files = []
+        for x in range(xmin, xmax + 1, 1):
+            for y in range(ymin, ymax + 1, 1):
+                lon = x - 360
+                lat = y - 180
+                files.append(self._get_dataset_url(lon, lat))
+
+        return files
+
+    def assets_for_point(
+        self,
+        lon: float,
+        lat: float,
+        coord_crs: CRS = WGS84_CRS,
+    ) -> List[str]:
+        """Retrieve assets for point."""
+        if coord_crs != WGS84_CRS:
+            xs, ys = transform_coords(coord_crs, WGS84_CRS, [lon], [lat])
+            lon, lat = xs[0], ys[0]
+
+        return [self._get_dataset_url(lon, lat)]
+
+    def _get_dataset_url(self, lon: float, lat: float) -> str:
+        """Return dataset url."""
+        northsouth = "N" if lat >= 0 else "S"
+        eastwest = "W" if lon < 0 else "E"
+
+        lat = abs(math.floor(lat))
+        lon = abs(math.floor(lon))
+
+        dataset = self.prefix_pattern.format(
+            northsouth=northsouth, eastwest=eastwest, lat=lat, lon=lon
+        )
+        return f"{self._scheme}://{self.bucket}/{dataset}.tif"
+
+
+@attr.s
+class Dem90Reader(Dem30Reader):
+    """Simple Mosaic reader for copernicus DEM"""
+
+    input: str = attr.ib(default="copernicus-dem-90m")
+
+    tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
+
+    minzoom: int = attr.ib(default=6)
+    maxzoom: int = attr.ib(default=7)
+
+    bounds: BBox = attr.ib(init=False, default=(-180, -90, 180, 90))
+    crs: CRS = attr.ib(init=False, default=WGS84_CRS)
+
+    reader: Type[Reader] = attr.ib(default=Reader)
+    reader_options: Dict = attr.ib(factory=dict)
+
+    colormap: Dict = attr.ib(init=False, factory=dict)
+
+    _scheme: str = "s3"
+    bucket: str = attr.ib(default="copernicus-dem-90m")
+    prefix_pattern: str = attr.ib(
+        default="Copernicus_DSM_COG_30_{northsouth}{lat:02d}_00_{eastwest}{lon:03d}_00_DEM/Copernicus_DSM_COG_30_{northsouth}{lat:02d}_00_{eastwest}{lon:03d}_00_DEM"
+    )
+
+    def statistics(self, **kwargs: Any) -> Dict[str, BandStatistics]:
+        """Return Dataset's statistics."""
+        # FOR NOW WE ONLY RETURN VALUE FROM THE FIRST FILE
+        dataset = self.prefix_pattern.format(
+            northsouth="S", eastwest="W", lat=90, lon=164
+        )
+        with Reader(f"{self._scheme}://{self.bucket}/{dataset}.tif") as src:
+            return src.statistics(**kwargs)

rio_tiler_pds/sentinel/aws/sentinel1.py

@@ -41,7 +41,7 @@ class S1L1CReader(MultiBandReader):
     maxzoom: int = attr.ib(default=14)
 
     reader: Type[Reader] = attr.ib(default=Reader)
-    reader_options: Dict = attr.ib(default={"options": {"nodata": 0}})
+    reader_options: Dict = attr.ib(factory=dict)
 
     productInfo: Dict = attr.ib(init=False)
     datageom: Dict = attr.ib(init=False)