A tile proxy server for USGS 3DEP Lidar Digital Surface Model dataset from Microsoft Planetary Computer

3dep-tile-proxy-api

from fastapi import FastAPI
from fastapi.middleware.cors import CORSMiddleware
from starlette.responses import StreamingResponse

import rasterio
import io
from PIL import Image
import numpy as np

app = FastAPI()

origins = ["*"]

app.add_middleware(
    CORSMiddleware,
    allow_origins=origins,
    allow_credentials=True,
    allow_methods=["*"],
    allow_headers=["*"],
)


@app.get("/{z}/{x}/{y}")
async def read_item(x, y, z):
    searchid = "dcff24ae650c9ffa0c392cf964efb0c4"
    url = f"https://planetarycomputer.microsoft.com/api/data/v1/mosaic/tiles/{searchid}/{z}/{x}/{y}?scale=1&format=tif&pixel_selection=first&TileMatrixSetId=WebMercatorQuad&assets=data&collection=3dep-lidar-dsm&unscale=False&resampling=nearest&return_mask=True"

    with rasterio.open(url) as src:
        dem = src.read(1)
        dem[dem == -9999] = 0

        r = np.zeros(dem.shape)
        g = np.zeros(dem.shape)
        b = np.zeros(dem.shape)

        v = dem + 32768
        r += np.floor(v / 256.0)
        g += np.floor(v % 256.0)
        b += np.floor((v - np.floor(v)) * 256.0)

        stack = np.dstack(
            [
                r.astype(rasterio.uint8),
                g.astype(rasterio.uint8),
                b.astype(rasterio.uint8),
            ]
        )

        im = Image.fromarray(stack)

        img_bytes = io.BytesIO()
        im.save(img_bytes, "PNG")
        img_bytes.seek(0)

        headers = {
            "Content-Type": "image/png",
            "Access-Control-Allow-Origin": "*",
            "Access-Control-Allow-Methods": "OPTIONS,GET",
            "Access-Control-Allow-Headers": "Content-Type",
        }

        return StreamingResponse(img_bytes, headers=headers)