Creación de servicios y datos de elevación

2023-03-23

En esta guía se describe cómo usar los datos de USGS en todo el mundo de DEM de su misión SRTM con una precisión de 30 m para crear un servicio Elevation en la nube de Microsoft Azure.

En este artículo se describe cómo:

Cree mosaicos vectoriales de línea de contorno y mosaicos DEM codificados en RGB.
Cree la API de elevación mediante mosaicos DEM codificados con Azure Functions y RGB a partir de Azure Blob Storage.
Cree un servicio de mosaico de vectores de línea de contorno mediante Azure Function y PostgreSQL.

Prerrequisitos

Esta guía requiere el uso del siguiente software y datos de terceros:

Datos de USGS. Los datos DEM se pueden descargar como GeoTiff con una cobertura de 1 segundo de arco por mosaico a través del USGS EarthExplorer. Esto requiere una cuenta de EarthExplorer, pero los datos se pueden descargar de forma gratuita.
La aplicación GIS de escritorio QGIS se usa para procesar y suavizar los mosaicos de trama. Tanto la descarga como el uso de QGIS son gratuitos. En esta guía se usa QGIS versión 3.26.2-Buenos Aires.
El paquete de Python rio-rgbify, desarrollado por MapBox, se usa para codificar GeoTIFF como RGB.
Base de datos PostgreSQL con la extensión espacial PostGIS.

Cree mosaicos vectoriales de línea de contorno y mosaicos DEM codificados en RGB

En esta guía se usan los 36 mosaicos que cubren el estado de Washington, disponible en USGS EarthExplorer.

Descarga de mosaicos de trama de USGS EarthExplorer

Criterios de búsqueda

Seleccione la región para la que desea los mosaicos de trama. Con fines de demostración, esta guía usa el método "Polygon" para seleccionar la región en el mapa.

Vaya al USGS EarthExplorer.
En la pestaña Criterios de búsqueda, seleccione Polígono y después seleccione cualquier punto del mapa para crear el límite.

Conjuntos de datos

Seleccione la pestaña Conjunto de datos.
Seleccione SRTM 1 Arc-Second Global en la sección Elevaciones digitales.

Results

Seleccione Resultados >> para ver los mosaicos de la región y el conjunto de datos seleccionados.
La lista de mosaicos descargables aparece en la página de resultados. Para descargar solo mosaicos que quiera, seleccione el botón Opciones de descarga de la tarjeta de resultados de cada mosaico, seleccionando la opción GeoTIFF 1 Arc-Second y repita este paso para los mosaicos restantes.
Como alternativa, use la opción de descarga masiva y seleccione GeoTIFF 1 Arc-second.

Adición de mosaicos de trama a QGIS

Una vez que tenga los mosaicos de trama que necesita, puede importarlos en QGIS.

Agregue mosaicos de trama a QGIS arrastrando los archivos a la pestaña capa de QGIS o seleccionando Agregar capa en el menú Capa.
Cuando las capas de trama se cargan en QGIS, puede haber diferentes tonos de mosaicos. Corrija esto combinando las capas de trama, lo que da como resultado una sola imagen de trama suave en formato GeoTIFF. Para ello, seleccione Varios en el menú Trama y, a continuación, Combinar...
Para cambiar el proyecto de la capa de trama combinada a EPSG:3857 (WGS84 / Pseudo-Mercator), haga clic con el botón derecho en Guardar capa de trama como para acceder a la capa de trama combinada de la opción tabla de contenido ->Exportar ->Guardar como. EPSG:3857 es necesario para usarlo con Azure Maps SDK web.
Si solo desea crear mosaicos vectoriales de línea de contorno, puede omitir los pasos siguientes y ir a Creación del servicio de mosaico de vectores de línea de contorno mediante Azure Function y PostgreSQL.
Para crear una API de elevación, el siguiente paso es RGB-Encode GeoTIFF. Esto se puede hacer mediante rio-rgbify, desarrollado por MapBox. Hay algunos desafíos al ejecutar esta herramienta directamente en Windows, por lo que es más fácil ejecutar desde WSL. A continuación se muestran los pasos descritos en Ubuntu en WSL:
```
sudo apt get update
sudo apt get upgrade
sudo apt install python3-pip
pip install rio-rgbify
PATH="$PATH:/home/<user  /.local/bin"
# The following two steps are only necessary when mounting an external hard drive or USB flash drive:
sudo mkdir /mnt/f
sudo mount -t drvfs D: /mnt/f

rio rgbify -b -10000 -i 0.1 wa_1arc_v3_merged_3857.tif wa_1arc_v3_merged_3857_rgb.tif

# The following steps are only necessary when unmounting an external hard drive or USB flash drive:
cd \~
sudo umount /mnt/f/
```
Los valores GeoTIFF codificados en RGB permiten recuperar valores R, G y B para un píxel y calcular la elevación a partir de estos valores:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
A continuación, cree un conjunto de mosaicos para usarlo con el control de mapa o úselo para obtener elevación de las coordenadas geográficas dentro de la extensión del mapa del conjunto de mosaicos. El conjunto de mosaicos se puede crear en QGIS mediante la herramienta Generar iconos XYZ (directorio).
Guarde la ubicación del conjunto de mosaicos y lo usará en la sección siguiente.

Crear la API de elevación mediante mosaicos DEM codificados con Azure Functions y RGB a partir de Azure Blob Storage

Los mosaicos DEM codificados RGB deben cargarse en un almacenamiento de base de datos para poder usarse con el Azure Functions para crear una API.

Cargue los mosaicos en Azure Blob Storage. Explorador de Azure Storage es una herramienta útil para este propósito.

La carga de mosaicos en Azure Blob Storage puede tardar varios minutos en completarse.
Una vez completada la carga, puede crear una función de Azure para crear una API que devuelva la elevación de una coordenada geográfica determinada.

Esta función recibe un par de coordenadas, determina el mosaico que lo cubre en el nivel de zoom 14 y, a continuación, determina las coordenadas de píxel dentro de ese mosaico que coincide con las coordenadas geográficas. A continuación, recupera el mosaico, obtiene los valores RGB de ese píxel y, a continuación, usa la fórmula siguiente para determinar la elevación:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)

import logging
import json
import azure.functions as func
from PIL import Image
import requests
from io import BytesIO
import math

def main(req: func.HttpRequest) -> func.HttpResponse:
    logging.info('Python HTTP trigger function processed a request.')

    # http://localhost:7071/api/GetElevationPoint?lng=-122.01911&lat=47.67091
    zoom = 14
    lng = float(req.params.get('lng'))
    lat = float(req.params.get('lat'))
    logging.info('Lng: ' + str(lng) + ' / lat: ' + str(lat))

    # Calculate th global pixel x and y coordinates for a lng / lat
    gx = (lng + 180) / 360
    sinLat = math.sin(lat * math.pi / 180)
    gy = 0.5 - math.log((1 + sinLat) / (1 - sinLat)) / (4 * math.pi)
    mapSize = math.ceil(256 * math.pow(2, zoom));
    gxc = min(max(gx * mapSize + 0.5, 0), mapSize - 1);            
    gyc = min(max(gy * mapSize + 0.5, 0), mapSize - 1);            

    # Calculate the tile x and y covering the lng / lat
    tileX = int(gxc / 256)
    tileY = int(gyc / 256)

    # Calculate the pixel coordinates for the tile covering the lng / lat
    tilePixelX = math.floor(gxc - (tileX * 256))
    tilePixelY = math.floor(gyc - (tileY * 256))

    response = requests.get("{BlobStorageURL}" + str(zoom) + "/" + str(tileX) + "/" + str(tileY) + ".png")
    im = Image.open(BytesIO(response.content))

    pix = im.load()
    r = pix[tilePixelX,tilePixelY][0]
    g = pix[tilePixelX,tilePixelY][1]
    b = pix[tilePixelX,tilePixelY][2]

    # elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
    ele = -10000 + ((r * 256 * 256 + g * 256 + b) * 0.1)

    jsonRes = {"elevation": + ele}
    logging.info('Response: ' + json.dumps(jsonRes))

    if lng and lat:
        return func.HttpResponse(
            json.dumps(jsonRes),
            mimetype="application/json",
        )
    else: 
        return func.HttpResponse(
            "ERROR: Missing parameter!",
            status_code=400
        )

Para ver los resultados del ejemplo de código, ejecútelo localmente:

localhost:7071/api/GetElevationPoint?lng=-122.01911&lat=47.67091`

Crear un servicio de mosaico de vectores de línea de contorno mediante Azure Function y PostgreSQL

En esta sección se describen los pasos para crear y procesar líneas de contorno en QGIS, cargarlas en PostgreSQL y, a continuación, crear una función de Azure para consultar PostgreSQL para devolver mosaicos vectoriales.

En QGIS, abra los mosaicos de trama combinados en la proyección EPSG:4326 creada en el paso 3 de Crear mosaicos vectoriales de línea de contorno y mosaicos DEM codificados en RGB.
Seleccione Extracción -> Contorno en el menú Trama para abrir la herramienta Contorno.

Al seleccionar Ejecutar se crean líneas de contorno y se agregan como una capa al mapa. algunos de los bordes de la línea de contorno pueden aparecer un poco rugoso. Esto se tratará en el siguiente paso.
Seleccione Cuadro de herramientas en el menú Procesamiento para abrir el Cuadro de herramientas de procesamiento.
A continuación, seleccione Smooth en la sección Geometría vectorial del Cuadro de herramientas de procesamiento.

Nota

El suavizado de líneas de contorno puede mejorarse considerablemente, pero a costa del aumento del tamaño de archivo.
Cargue las líneas de contorno en la base de datos. En esta guía se usa la versión gratuita de la base de datos PostgreSQL que se ejecuta en localhost. También puede cargarlos en el Azure Database for PostgreSQL.

El siguiente paso requiere una base de datos PostgreSQL con la extensión PostGIS.
Para crear una conexión de QGIS a PostgreSQL, seleccione Agregar capa ->Agregar capas PostGIS en el menú Capa y, a continuación, seleccione el botón Nuevo.
A continuación, cargue datos de QGIS a PostgreSQL mediante el Administrador de bases de datos en QGIS. Para ello, seleccione Administrador de bases de datos en el menú Base de datos.
Conéctese a la base de datos de PostGIS y seleccione Importar capa/archivo... para importar líneas de contorno a la base de datos.

Ahora puede usar una función de Azure para consultar PostgreSQL y devolver mosaicos vectoriales para las líneas de contorno. El servidor de mosaicos se puede usar con el SDK web de Azure Maps para crear una aplicación web que muestre líneas de contorno en el mapa.

import logging
from wsgiref import headers
import azure.functions as func
import psycopg2
# Database to connect to
DATABASE = {
    'user':     'postgres',
    'password': '{password}',
    'host':     'localhost',
    'port':     '5432',
    'database': '{database}'
    }
def main(req: func.HttpRequest) -> func.HttpResponse:
    logging.info('Python HTTP trigger function processed a request.')
    DATABASE_CONNECTION = None
    # get url parameters http://localhost:7071/api/tileserver?zoom={z}&x={x}&y={y}
    # http://localhost:7071/api/tileserver?zoom=16&x=10556&y=22870
    zoom = int(req.params.get('zoom'))
    x = int(req.params.get('x'))
    y = int(req.params.get('y'))
    table = req.params.get('table')
    # calculate the envelope of the tile
    # Width of world in EPSG:3857
    worldMercMax = 20037508.3427892
    worldMercMin = -1 * worldMercMax
    worldMercSize = worldMercMax - worldMercMin
    # Width in tiles
    worldTileSize = 2 ** zoom

    # Tile width in EPSG:3857
    tileMercSize = worldMercSize / worldTileSize

    # Calculate geographic bounds from tile coordinates
    # XYZ tile coordinates are in "image space" so origin is
    # top-left, not bottom right
    xmin = worldMercMin + tileMercSize * x
    xmax = worldMercMin + tileMercSize * (x + 1)
    ymin = worldMercMax - tileMercSize * (y + 1)
    ymax = worldMercMax - tileMercSize * y
    # Generate SQL to materialize a query envelope in EPSG:3857.
    # Densify the edges a little so the envelope can be
    # safely converted to other coordinate systems.
    DENSIFY_FACTOR = 4
    segSize = (xmax - xmin)/DENSIFY_FACTOR
    sql01 = 'ST_Segmentize(ST_MakeEnvelope(' + str(xmin) + ', ' + str(ymin) + ', ' + str(xmax) + ', ' + str(ymax) + ', 3857), ' + str(segSize) +')'

    # Generate a SQL query to pull a tile worth of MVT data
    # from the table of interest. 
    # Materialize the bounds
    # Select the relevant geometry and clip to MVT bounds
    # Convert to MVT format
    sql02 = 'WITH bounds AS (SELECT ' + sql01 + ' AS geom, ' + sql01 + '::box2d AS b2d), mvtgeom AS (SELECT ST_AsMVTGeom(ST_Transform(t.geom, 3857), bounds.b2d) AS geom, elev FROM contourlines_smooth t, bounds WHERE ST_Intersects(t.geom, ST_Transform(bounds.geom, 4326))) SELECT ST_AsMVT(mvtgeom.*) FROM mvtgeom' 

    # Run tile query SQL and return error on failure conditions
    # Make and hold connection to database
    if not DATABASE_CONNECTION:
        try:
            DATABASE_CONNECTION = psycopg2.connect(**DATABASE)
            logging.info('Connected to database.')
        except (Exception, psycopg2.Error) as error:
            logging.error('ERROR: Cannot connect to database.')
    # Query for MVT
    with DATABASE_CONNECTION.cursor() as cur:
        cur.execute(sql02)
        if not cur:
            logging.error('ERROR: SQL Query failed.')
        pbf = cur.fetchone()[0]
        logging.info('Queried database')

    if zoom and x and y:
        return func.HttpResponse(
            # f"This HTTP triggered function executed successfully.\n\nzoom={zoom}\nx={x}\ny={y}\n\nxmin={xmin}\nxmax={xmax}\nymin={ymin}\nymax={ymax}\n\nsql01={sql01}\n\nsql02={sql02}",
            bytes(pbf),
            status_code=200,
            headers={"Content-type": "application/vnd.mapbox-vector-tile","Access-Control-Allow-Origin": "*"}
        )
    else: 
        return func.HttpResponse(
            "ERROR: Missing parameter!",
            status_code=400
        )