Criar dados de elevação e serviços

2023-03-24

Este guia descreve como usar dados DEM mundiais do USGS a partir da missão SRTM com precisão de 30m para criar um serviço Elevação na Nuvem do Microsoft Azure.

Este artigo descreve como:

Crie Peças vetoriais de curvas de nível e Peças DEM codificadas em RGB.
Crie a API de Elevação usando a Função do Azure e as Peças DEM codificadas em RGB do Armazenamento de Blobs do Azure.
Crie um serviço de Peça vetorial de curvas de nível usando a Função do Azure e o PostgreSQL.

Pré-requisitos

Este guia requer o uso do seguinte software e dados de terceiros:

Dados USGS. Os dados DEM podem ser baixados como GeoTiff com cobertura de 1 segundo de arco por peça de mapa por meio do USGS EarthExplorer. Isso requer uma conta no EarthExplorer, mas os dados podem ser baixados gratuitamente.
O aplicativo GIS QGIS Desktop é usado para processar e suavizar as peças de mapa raster. O QGIS é gratuito para download e uso. Este guia usa o QGIS versão 3.26.2-Buenos Aires.
O pacote rio-rgbify do Python, desenvolvido pela MapBox, é usado para codificar o GeoTIFF como RGB.
Banco de dados do PostgreSQL com a extensão espacial PostGIS.

Criar Peças vetoriais de curvas de nível e Peças DEM codificados em RGB

Este guia usa os 36 peças de mapa que abrangem o estado de Washington, disponíveis no USGS EarthExplorer.

Baixar peças raster no USGS EarthExplorer

Critérios de pesquisa

Selecione a região na qual você deseja as peças de mapa raster. Para fins de demonstração, este guia usa o método "Polígono" para selecionar a região no mapa.

Navegue até o site USGS EarthExplorer.
Na guia Search Criteria, selecione Polygon e selecione qualquer ponto no mapa para criar o limite.

Conjuntos de dados

Selecione a guia Data Sets.
Selecione SRTM 1 Arc-Second Global na seção Digital Elevation.

Resultados

Selecione Results>> para exibir as peças de mapa e os conjunto de dados da região selecionada.
A lista das peças de mapa baixáveis aparece na página de resultados. Para baixar apenas as peça de mapa desejadas, selecione o botão Download Options no quadro de resultados de cada peça, selecionando a opção GeoTIFF 1 Arc-Second e repita esta etapa para as peças restantes.
Alternativamente, use a opção download em massa e selecione GeoTIFF 1 Arc-second.

Adicione as peças raster ao QGIS

Depois de ter as peças raster necessárias, você poderá importá-las no QGIS.

Adicione peças de mapa raster ao QGIS arrastando os arquivos para a guia Camada QGIS ou selecionando Adicionar Camada no menu Camada.
Quando as camadas raster são carregadas no QGIS, pode haver diferentes tons de peças de mapa. Corrija isso mesclando as camadas raster, o que resulta em uma única imagem rasterizada uniforme no formato GeoTIFF. Para fazer isso, selecione Miscelânea no menu Raster e, em seguida, Mesclar...
Reprojete a camada raster mesclada para EPSG:3857 (WGS84 /Pseudo-Mercator) usando Salvar Camada Raster como clique com o botão direito do mouse para acessar a camada raster mesclada na tabela de conteúdos as opções –>Exportar –>Salvar Como. O EPSG:3857 é necessário para usá-lo com o SDK Web do Azure Mapas.
Se você deseja apenas criar peças vetoriais de curvas de nível, pule as etapas a seguir e vá para Criar serviço de Peça vetorial de curvas de nível usando a Função do Azure e o PostgreSQL.
Para criar uma API de Elevação, a próxima etapa é codificar o GeoTIFF em RGB. Isso pode ser feito usando o recurso rio-rgbify, desenvolvido pela MapBox. Existem alguns desafios ao executar essa ferramenta diretamente no Windows, portanto, é mais fácil executá-la no WSL. Abaixo estão as etapas no WSL no Ubuntu:
```
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 recuperar os valores R, G e B de um pixel e calcular a elevação a partir desses valores:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Em seguida, crie um conjunto de peças de mapa para usar com o controle de mapa e/ou usá-lo para obter a Elevação de qualquer coordenada geográfica dentro da extensão do mapa do conjunto de peças de mapa. O conjunto de peças de mapa pode ser criado no QGIS usando a ferramenta Gerar peças de mapa XYZ (Diretório).
Salve o local do conjunto de peças de mapa, você o usará na próxima Seção.

Criar API de Elevação usando a Função do Azure e as Peças DEM codificadas em RGB do Armazenamento de Blobs do Azure

Os Peças DEM codificadas em RGB precisam ser carregados em um armazenamento de banco de dados antes de serem usados com o Azure Functions para criar uma API.

Carregue as peças de mapa no Armazenamento de Blobs do Azure. O Gerenciador de Armazenamento do Azure é uma ferramenta útil para essa finalidade.

O carregamento das peça de mapa no Armazenamento de Blobs do Azure pode levar alguns minutos para concluir.
Depois que o upload for concluído, você poderá criar a Função do Azure para desenvolver uma API que retorna a elevação de uma determinada coordenada geográfica.

Essa função recebe um par de coordenadas, determina a peça de mapa que a cobre no nível de zoom 14 e, em seguida, determina as coordenadas de pixel dentro dessa peça que correspondem às coordenadas geográficas. Em seguida, recupera a peça, obtém os valores RGB desse pixel e, em seguida, 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 código de exemplo, execute-o localmente:

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

Crie um serviço de Peça vetorial de curvas de nível usando a Função do Azure e o PostgreSQL

Esta seção descreve as etapas para criar e processar curvas de nível no QGIS, carregá-las no PostgreSQL e, em seguida, criar uma Função do Azure para consultar o PostgreSQL e retornar peças vetoriais.

No QGIS, abra as peças raster mesclados na projeção EPSG:4326 criada na etapa 3 do tutorial Criar Peças vetoriais de curvas de nível e Peças DEM codificadas em RGB.
Selecione Extração –> Contorno no menu Raster para abrir a ferramenta Contorno.

Ao selecionar Executar, as curvas de nível serão criadas e adicionadas como uma camada no mapa. alguns dos contornos das curvas de nível podem parecer um pouco irregulares. Isso será abordado na próxima etapa.
Selecione a Caixa de ferramentas no menu Processamento para abrir a Caixa de Ferramentas de Processamento.
Em seguida, selecione Suavização na seção Geometria do vetor da Caixa de Ferramentas de Processamento.

Observação

A suavização das curvas de nível pode ser substancialmente aprimorada, mas isso aumentará o tamanho do arquivo.
Carregue as curvas de nível no banco de dados. Este guia usa a versão gratuita do banco de dados do PostgreSQL, que é executada no localhost. Você também pode carregá-las no Banco de Dados do Azure para PostgreSQL.

A próxima etapa requer um banco de dados do PostgreSQL com a extensão PostGIS.
Para criar uma conexão do QGIS com o PostgreSQL, selecione Adicionar Camada –>Adicionar Camada PostGIS no menu Camada e, em seguida, selecione o botão Novo.
Em seguida, carregue os dados do QGIS no PostgreSQL usando o Gerenciador de Banco de Dados no QGIS. Para fazer isso, selecione o Gerenciador de BD no menu Banco de dados.
Conecte-se ao banco de dados do PostGIS e selecione Importar Camada/Arquivo... para Importar as curvas de nível para o banco de dados.

Agora você pode usar uma Função do Azure para consultar o PostgreSQL e retornar peças vetoriais para as curvas de nível. O servidor de peças de mapa pode ser usado com o SDK web do Azure Mapas para criar um aplicativo web que exibe as curvas de nível 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
        )