Creare dati e servizi di elevazione

2025-04-30

Questa guida descrive come usare i dati DEM globali USGS dalla loro missione SRTM con un'accuratezza di 30m per creare un servizio di elevazione in Microsoft Azure Cloud.

In questo articolo viene descritto come:

Creazione di riquadri vettoriali di linee di contorno e riquadri DEM codificati in RGB.
Creare l'API di elevazione dei privilegi usando i riquadri DEM con codifica RGB e funzione di Azure da Archiviazione BLOB di Azure.
Creare un servizio di tile vettoriali di linee di contorno utilizzando le Funzioni di Azure e PostgreSQL.

Prerequisiti

Questa guida richiede l'uso dei seguenti dati e software di terze parti:

Dati USGS. I dati DEM possono essere scaricati come GeoTiff con una copertura di 1 secondo di arco per riquadro tramite USGS EarthExplorer. Ciò richiede un account EarthExplorer, ma i dati possono essere scaricati gratuitamente.
L'applicazione GIS QGIS desktop viene usata per elaborare e smussare i riquadri Raster. QGIS è gratuito da scaricare e usare. Questa guida usa QGIS versione 3.26.2-Buenos Aires.
Il pacchetto Python rio-rgbify, sviluppato da MapBox, viene usato per codificare GeoTIFF come RGB.
Database PostgreSQL con l'estensione spaziale PostGIS.

Creare riquadri vettoriali linea contorno e riquadri DEM con codifica RGB

Questa guida usa i 36 riquadri che coprono lo stato di Washington, disponibili da USGS EarthExplorer.

Scaricare riquadri raster da USGS EarthExplorer

Criteri di ricerca

Selezionare l'area per cui si desiderano riquadri raster. A scopo dimostrativo, questa guida usa il metodo "Poligono" per selezionare l'area sulla mappa.

Navigare su USGS EarthExplorer.
Nella scheda Criteri di ricerca selezionare Poligono, quindi selezionare un punto qualsiasi della mappa per creare il limite.

Set di dati

Selezionare la scheda Set di dati.
Selezionare SRTM 1 Arc-Second Global nella sezione Elevazioni digitali.

Risultati

Selezionare Risultati >> per visualizzare i riquadri per l'area e il set di dati selezionati.
L'elenco dei riquadri scaricabili viene visualizzato nella pagina dei risultati. Per scaricare solo i riquadri desiderati, selezionare il pulsante Opzioni di download nella scheda dei risultati per ogni riquadro, selezionare l'opzione GeoTIFF 1 Arc-Second e ripetere questo passaggio per i riquadri rimanenti.
In alternativa, usare l'opzione di download bulk e selezionare GeoTIFF 1 Arc-second.

Aggiungere riquadri raster a QGIS

Dopo aver creato i riquadri raster necessari, è possibile importarli in QGIS.

Aggiungere riquadri raster a QGIS trascinando i file nella scheda Livello QGIS o selezionando Aggiungi livello nel menu Livello.
Quando i livelli raster vengono caricati in QGIS, possono esserci sfumature diverse di riquadri. Risolvere questo problema unendo i livelli raster, cosa che dà vita a un'unica immagine raster uniforme in formato GeoTIFF. A tale scopo, selezionare Varie dal menu Raster, quindi Merge...
Riprogettare il livello raster unito a EPSG:3857 (WGS84 / Pseudo-Mercator) usando Salva livello raster come, fare clic con il pulsante destro del mouse per accedere al livello raster unito nell'opzione tabella di contenuto ->Esporta ->Salva con nome opzione. EPSG:3857 è necessario per l'uso con l'SDK Web di Mappe di Azure.
Se si vogliono creare solo tessere vettoriali delle linee di contorno, è possibile ignorare i passaggi seguenti e passare a Creare un servizio di tessere vettoriali delle linee di contorno usando Funzioni di Azure e PostgreSQL.
Per creare un'API Elevazione, il passaggio successivo consiste nel codificare in RGB il GeoTIFF. Questa operazione può essere eseguita usando rio-rgbify, sviluppato da MapBox. L'esecuzione di questo strumento direttamente in Windows presenta alcuni problemi, quindi è più facile eseguirlo da WSL. Di seguito sono riportati i passaggi su Ubuntu su 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/
```
Il GeoTIFF con codifica RGB consente di recuperare i valori R, G e B per un pixel e di calcolare l'elevazione da questi valori:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Creare quindi un set di tile da utilizzare con il controllo della mappa e/o utilizzarlo per ottenere l'altitudine di qualsiasi coordinate geografica all'interno dell'estensione della mappa del set di tile. Il set di riquadri può essere creato in QGIS usando lo strumento Genera riquadri XYZ (Directory).
Salva la posizione del set di riquadri, che userai nella sezione successiva.

Creare l'API di elevazione dei privilegi usando i riquadri DEM con codifica RGB e funzione di Azure da Archiviazione BLOB di Azure

I riquadri DEM con codifica RGB devono essere caricati in un archivio di database prima di poterli usare con Funzioni di Azure per creare un'API.

Caricare i riquadri in Archiviazione BLOB di Azure. Azure Storage Explorer è uno strumento utile per questo scopo.

Il completamento del caricamento dei riquadri in Archiviazione BLOB di Azure può richiedere alcuni minuti.
Al termine del caricamento, è possibile creare una funzione di Azure per creare un'API che restituisca l'elevazione per una determinata coordinata geografica.

Questa funzione riceve una coppia di coordinate, determina il riquadro che la copre al livello di zoom 14, quindi determina le coordinate pixel all'interno di tale riquadro che corrispondono alle coordinate geografiche. Recupera quindi il riquadro, ottiene i valori RGB per tale pixel e usa la formula seguente per determinare l'elevazione:

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
        )

Per visualizzare i risultati dell'esempio di codice, eseguirlo in locale:

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

Creare un servizio riquadro vettoriale linea contorno con Funzioni di Azure e PostgreSQL

Questa sezione descrive i passaggi per creare ed elaborare linee di contorno in QGIS, caricarle in PostgreSQL e quindi creare una funzione di Azure in Query PostgreSQL per restituire riquadri vettoriali.

In QGIS aprire i riquadri raster uniti nella proiezione EPSG:4326 creata nel passaggio 3 di Creare riquadri vettoriali linea di contorno e riquadri DEM con codifica RGB.
Selezionare Estrazione -> Contorno dal menu Raster per aprire lo strumento Contorno.

Selezionando Esegui si creano linee di contorno e le si aggiunge come livello alla mappa. Alcuni dei bordi della linea di contorno potrebbero apparire un po' grezzi. Questo fatto verrà risolto nel passaggio successivo.
Selezionare Casella degli strumenti dal menu Elaborazione per visualizzare la Casella degli strumenti di elaborazione.
Selezionare quindi Fluida nella sezione Geometria vettore della Casella degli strumenti di elaborazione.

Nota

Lo smussamento della linea di contorno può essere notevolmente migliorato, ma a costo di aumentare le dimensioni del file.
Caricare le linee di contorno nel database. Questa guida usa la versione gratuita del database di PostgreSQL eseguito in localhost. È anche possibile caricarli nel Database di Azure per PostgreSQL.

Il passaggio successivo richiede un database PostgreSQL con estensione PostGIS.
Per creare una connessione da QGIS a PostgreSQL, selezionare Aggiungi livello ->Aggiungi livelli PostGIS dal menu Livello, quindi selezionare il pulsante Nuovo.
Carica quindi i dati da QGIS a PostgreSQL usando il Gestione Database in QGIS. A tale scopo, selezionare DB Manager dal menu Database.
Connettersi al database PostGIS e selezionare Importa livello/file... per importare righe di contorno nel database.

È ora possibile usare una funzione di Azure per eseguire query su PostgreSQL e restituire riquadri vettoriali per le linee di contorno. Il server riquadri può essere usato con l'SDK Web di Mappe di Azure per creare un'app Web che visualizza linee di contorno sulla mappa.

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
        )