add copernicus-dem reader ?

[vincentsarago]

https://registry.opendata.aws/copernicus-dem/

[vincentsarago]

"""custom readers."""

import json
import pathlib
from typing import Any, Dict, List, Tuple, Type

import attr
from morecantile import TileMatrixSet

from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
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.errors import TileOutsideBounds, PointOutsideBounds
from rasterio.errors import RasterioIOError

template = "s3://copernicus-dem-30m/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.tif"


@attr.s
class CopernicusDemReader(BaseReader):
    """Simple Mosaic reader for copernicus DEM"""
    input: str = attr.ib(default="copernicus-dem-30m")
    tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)

    reader: Type[Reader] = attr.ib(default=Reader)

    colormap: Dict = attr.ib(init=False)

    def __attrs_post_init__(self):
        """Fetch Reference band to get the bounds."""
        self.minzoom = 7
        self.maxzoom = 8
        self.crs = WGS84_CRS
        self.bounds = (-180, -90, 180, 90)
        self.colormap = {}

    def tile(
        self,
        tile_x: int,
        tile_y: int,
        tile_z: int,
        reverse: bool = False,
        **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) 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=0,
            allowed_exceptions=(TileOutsideBounds, RasterioIOError),
            **kwargs,
        )[0]

    def point(
        self,
        lon: float,
        lat: float,
        reverse: bool = False,
        **kwargs: Any,
    ) -> PointData:
        """Get Point value."""
        assets = self.assets_for_point(lon, lat)
        mosaic_assets = list(reversed(assets)) if reverse else assets

        def _reader(asset: str, lon: float, lat: float, **kwargs) -> PointData:
            with Reader(asset) as src:
                return src.point(lon, lat, **kwargs)

        return mosaic_point_reader(
            mosaic_assets,
            _reader,
            lon,
            lat,
            threads=0,
            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
        item = "s3://copernicus-dem-30m/Copernicus_DSM_COG_10_N00_00_E006_00_DEM/Copernicus_DSM_COG_10_N00_00_E006_00_DEM.tif"
        with Reader(item) 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."""
        bbox = self.tms.bounds(x, y, z)

        left = int(bbox.left + 360)
        right = int(bbox.right + 360)

        top = int(bbox.top + 180)
        bottom = int(bbox.bottom + 180)

        files = []
        for x in range(left, right + 1, 1):
            for y in range(bottom, top + 1, 1):
                lng = x - 360
                lat = y - 180

                northsouth = "N" if lat >=0 else "S"
                eastwest = "W" if lng < 0 else "E"

                files.append(
                    template.format(
                        northsouth=northsouth,
                        eastwest=eastwest,
                        lat=round(abs(lat)),
                        lon=round(abs(lng)),
                    )
                )

        return files

    def assets_for_point(self, lng: float, lat: float) -> List[str]:
        """Retrieve assets for point."""
        northsouth = "N" if lat >=0 else "S"
        eastwest = "W" if lng < 0 else "E"
        return [
            template.format(
                northsouth=northsouth,
                eastwest=eastwest,
                lat=round(abs(lat)),
                lon=round(abs(lng)),
            )
        ]

works well ;-)