Criar dados e serviços de elevação

Este guia descreve como usar dados DEM do USGS a nível mundial da sua missão SRTM com precisão de 30m para construir um serviço Elevation na cloud Microsoft Azure.

Este artigo descreve como:

Crie mosaicos vetoriais de linhas de contorno e mosaicos DEM codificados em RGB.
Crie a API de elevação usando Azure Function e tiles DEM codificados em RGB a partir do Azure Blob Storage.
Criar um serviço de azulejos vetoriais de linhas de contorno usando Azure Functions e PostgreSQL.

Prerequisites

Este guia requer a utilização do seguinte software e dados de terceiros:

Dados do USGS. Os dados DEM podem ser descarregados em formato GeoTiff, com uma cobertura de 1 segundo de arco por mosaico, através do USGS EarthExplorer. Isto requer uma conta EarthExplorer, mas os dados podem ser descarregados gratuitamente.
A aplicação SIG de ambiente de trabalho QGIS é usada para processar e suavizar os tiles Raster. O QGIS é gratuito para descarregar e usar. Este guia utiliza a versão QGIS 3.26.2-Buenos Aires.
O pacote Python rio-rgbify , desenvolvido pela MapBox, é usado para codificar o GeoTIFF como RGB.
Base de dados PostgreSQL com a extensão espacial PostGIS.

Criar tiles vetoriais de linhas de contorno e tiles de DEM codificados em RGB

Este guia utiliza os 36 mosaicos que cobrem o estado de Washington, disponíveis no USGS EarthExplorer.

Descarregue blocos raster do USGS EarthExplorer

Critérios de pesquisa

Seleciona a região para a qual quer mosaicos raster. Para fins de demonstração, este guia utiliza o método "Polygon" para selecionar a região no mapa.

Navegue até ao USGS EarthExplorer.
No separador de Critérios de Pesquisa , selecione Polígono e depois selecione qualquer lugar no mapa para criar o limite.

Conjuntos de dados

Selecione o separador Conjuntos de Dados .
Selecione SRTM 1 Arc-Second Global na secção de Elevações Digitais .

Results

Selecione Resultados >> para visualizar os azulejos da região e do conjunto de dados selecionados.
A lista de azulejos descarregáveis aparece na página de resultados. Para descarregar apenas os tiles que quiser, selecione o botão Opções de Download no cartão de resultados para cada tile, selecione a opção GeoTIFF 1 Arc-Second e repita este passo para os tiles restantes.
Alternativamente, use a opção de download em massa e selecione GeoTIFF 1 Arc-second.

Adicionar mosaicos raster ao QGIS

Depois de teres os tiles raster de que precisas, podes importá-los no QGIS.

Adicione blocos raster ao QGIS arrastando os ficheiros para o separador de camadas do QGIS ou selecionando Adicionar Camada no menu de Camadas .
Quando as camadas raster são carregadas no QGIS, podem existir diferentes tons de tiles. Resolva isto fundindo as camadas raster, o que resulta numa única imagem raster suave em formato GeoTIFF. Para isso, selecione Diversos no menu Raster, depois Fundir…
Reprojete a camada raster fundida para EPSG:3857 (WGS84 / Pseudo-Mercator) clicando com o botão direito em Save Raster Layer as para aceder à camada raster fundida na tabela de conteúdo, depois selecione Exportar ->Guardar Como. EPSG:3857 é necessário para usar com o Azure Maps Web SDK.
Se só quiser criar mosaicos vetoriais de linhas de contorno, pode saltar os passos seguintes e ir ao serviço Criar mosaicos vetoriais de linha de contorno usando Azure Function e PostgreSQL.
Para criar uma API de Elevação, o passo seguinte é codificar o GeoTIFF em RGB. Isto pode ser feito usando rio-rgbify, desenvolvido pela MapBox. Existem alguns desafios em correr esta ferramenta diretamente no Windows, por isso é mais fácil de executar a partir do WSL. Abaixo estão os passos no Ubuntu no 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/
```
O GeoTIFF codificado em RGB permite obter valores R, G e B para um pixel e calcular a elevação a partir destes valores:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
De seguida, criar um conjunto de tiles para usar com o controlo do mapa e/ou usá-lo para obter a elevação de quaisquer coordenadas geográficas dentro da extensão do mapa do conjunto de tiles. O conjunto de tiles pode ser criado no QGIS usando a ferramenta Generate XYZ tiles (Diretório ).
Guarda a localização do conjunto de tiles, vais usá-lo na próxima Secção.

Criar API de elevação usando Azure Function e tiles DEM codificados em RGB a partir do Azure Blob Storage

Os Tiles DEM codificados RGB precisam de ser carregados para um armazenamento de base de dados antes de poderem ser usados com as Azure Functions para criar uma API.

Carrega os tiles para o Azure Blob Storage. O Azure Storage Explorer é uma ferramenta útil para este propósito.

O carregamento de tiles para o Azure Blob Storage pode demorar vários minutos a ser concluído.
Uma vez concluído o upload, podes criar uma Azure Function para construir uma API que devolve a elevação para uma dada coordenada geográfica.

Esta função recebe um par de coordenadas, determina o mosaico que o cobre ao nível de zoom 14, depois determina as coordenadas dos píxeis dentro desse mosaico que correspondem às coordenadas geográficas. Depois recupera o mosaico, obtém os valores RGB desse píxel e usa a seguinte fórmula para determinar a elevação:

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 os resultados do exemplo de código, execute-o localmente:

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

Criar um serviço de linhas de contorno em tiles vetoriais usando Azure Functions e PostgreSQL

Esta secção descreve os passos para criar e processar linhas de contorno no QGIS, carregá-las para o PostgreSQL e depois criar uma Função Azure para Consultar o PostgreSQL e devolver tiles vetoriais.

No QGIS, abra os tiles raster fundidos na projeção EPSG:4326 criada no passo 3 de Criar tiles vetoriais de linha de contorno e tiles DEM codificados RGB.
Selecione Extração -> Contorno no menu Raster para abrir a ferramenta Contorno.

Selecionar Executar cria linhas de contorno e adiciona-as como uma camada ao mapa. Algumas das bordas das linhas de contorno podem parecer um pouco ásperas. Isto será abordado na próxima etapa.
Selecione Caixa de Ferramentas no menu de Processamento para abrir a Caixa de Ferramentas de Processamento.
Depois seleciona Suavizar na secção Geometria Vetorial da Caixa de Ferramentas de Processamento.

Note

O suavizamento das linhas de contorno pode ser substancialmente melhorado, mas à custa do aumento do tamanho do ficheiro.
Carregue as linhas de contorno na base de dados. Este guia utiliza a versão gratuita da base de dados PostgreSQL que corre no localhost. Também podes carregá-las na base de dados Azure para PostgreSQL.

O passo seguinte requer uma base de dados PostgreSQL com extensão PostGIS .
Para criar uma ligação do QGIS para o PostgreSQL, selecione Adicionar Camada ->Adicionar Camadas PostGIS no menu de Camadas, depois selecione o botão Novo.
De seguida, carregue Dados do QGIS para PostgreSQL usando o Gestor de Bases de Dados no QGIS. Para isso, selecione Gestor de Bases de Dados no menu Base de Dados .
Liga-te à base de dados PostGIS e seleciona Importar Camada/Ficheiro... para importar linhas de contorno para a base de dados.

Agora pode usar uma Função Azure para consultar o PostgreSQL e devolver tiles vetoriais para as linhas de contorno. O servidor de tiles pode ser usado com o Azure Maps Web SDK para criar uma aplicação web que exibe as linhas de contorno no 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
        )