Erstellen von Höhenangabedaten und Höhenangabediensten

2023-03-27

In dieser Anleitung wird beschrieben, wie Sie die weltweiten USGS-DEM-Daten ihrer SRTM-Mission mit einer Genauigkeit von 30 m verwenden, um einen Höhenangabendienst in der Microsoft Azure-Cloud zu erstellen.

Dieser Artikel beschreibt, wie Sie:

Erstellen von Konturlinienvektorkacheln und RGB-codierten DEM-Kacheln.
Erstellen einer Höhenangaben-API mithilfe von Azure Function- und RGB-codierten DEM-Kacheln aus Azure Blob Storage.
Erstellen eines Konturlinienvektorkachel-Dienstes mithilfe von Azure Function und PostgreSQL.

Voraussetzungen

Diese Anleitung erfordert die Verwendung der folgenden Drittanbietersoftware und -daten:

USGS-Daten. DEM-Daten können als GeoTiff-Datei mit einer Abdeckung von 1 Bogensekunde pro Kachel über den USGS EarthExplorer heruntergeladen werden. Dies erfordert ein EarthExplorer-Konto, aber die Daten können kostenlos heruntergeladen werden.
Die GIS-Desktopanwendung QGIS wird verwendet, um die Rasterkacheln zu verarbeiten und zu glätten. QGIS kann kostenlos heruntergeladen und verwendet werden. In dieser Anleitung wird QGIS Version 3.26.2-Buenos Aires verwendet.
Das von MapBox entwickelte Python-Paket rio-rgbify wird verwendet, um die GeoTIFF-Datei als RGB zu codieren.
PostgreSQL-Datenbank mit der räumlichen Erweiterung PostGIS.

Erstellen von Konturlinienvektorkacheln und RGB-codierten DEM-Kacheln

In dieser Anleitung werden die 36 Kacheln verwendet, die den Bundesstaat Washington abdecken und im USGS EarthExplorer verfügbar sind.

Herunterladen von Rasterkacheln aus USGS EarthExplorer

Suchkriterien

Wählen Sie die Region aus, für die Sie Rasterkacheln haben möchten. Zur Veranschaulichung wird in dieser Anleitung die Methode „Polygon“ verwendet, um die Region auf der Karte auszuwählen.

Navigieren Sie zum USGS EarthExplorer.
Wählen Sie auf der Registerkarte Suchkriterien die Option Polygon aus, und wählen Sie dann einen beliebigen Teil der Karte aus, um die Grenze zu erstellen.

Datasets

Wählen Sie die Registerkarte Datasets aus.
Wählen Sie im Abschnitt Digitale Höhenangaben die Option SRTM 1 Bogensekunde Global aus.

Ergebnisse

Wählen Sie Ergebnisse >> aus, um die Kacheln für die ausgewählte Region und das ausgewählte Dataset anzuzeigen.
Die Liste der herunterladbaren Kacheln wird auf der Ergebnisseite angezeigt. Wenn Sie nur die gewünschten Kacheln herunterladen möchten, wählen Sie auf der Ergebniskarte für jede Kachel die Schaltfläche Downloadoptionen und dann die Option GeoTIFF 1 Bogensekunde aus, und wiederholen Sie diesen Schritt für die verbleibenden Kacheln.
Verwenden Sie alternativ die Massendownloadoption, und wählen Sie GeoTIFF 1 Bogensekunde aus.

Hinzufügen von Rasterkacheln zu QGIS

Sobald Sie über die benötigten Rasterkacheln verfügen, können Sie sie in QGIS importieren.

Fügen Sie Rasterkacheln zu QGIS hinzu, indem Sie die Dateien auf die Registerkarte QGIS-Ebene ziehen oder im Menü Ebene die Option Ebene hinzufügen auswählen.
Wenn die Rasterebenen in QGIS geladen werden, können die Kacheln verschiedene Schattierungen haben. Beheben Sie dies, indem Sie die Rasterebenen zusammenführen. Dies führt zu einem einzelnen glatten Rasterbild im GeoTIFF-Format. Wählen Sie dazu im Menü Raster die Option Verschiedenes und dann Zusammenführen... aus.
Projizieren Sie die zusammengeführte Rasterebene in EPSG:3857 (WGS84/Pseudo-Mercator) erneut, indem Sie Rasterebene speichern unter durch Klicken mit der rechten Maustaste auf die zusammengeführte Rasterebene im Inhaltsverzeichnis verwenden und die Option >Exportieren ->Speichern unter auswählen. EPSG:3857 ist für das Azure Maps-Web-SDK erforderlich.
Wenn Sie nur Konturlinienvektorkacheln erstellen möchten, können Sie die folgenden Schritte überspringen und mit Erstellen eines Konturlinienvektorkachel-Dienstes mithilfe einer Azure Function und von PostgreSQL fortfahren.
Zum Erstellen einer Höhenangaben-API wird die GeoTIFF-Datei im nächsten Schritt RGB-codiert. Dies kann mithilfe von rio-rgbify, entwickelt von MapBox, erfolgen. Da die Ausführung dieses Tools direkt in Windows mit gewissen Schwierigkeiten verbunden ist, ist es einfacher, es über WSL auszuführen. Im Folgenden finden Sie die Schritte in Ubuntu unter 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/
```
Mit der RGB-codierten GeoTIFF-Datei können Sie R-, G- und B-Werte für ein Pixel abrufen und die Höhe anhand dieser Werte berechnen:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Erstellen Sie als Nächstes einen Kachelsatz, der mit dem Kartensteuerelement verwendet werden soll, und/oder verwenden Sie ihn, um die Höhe für alle geografischen Koordinaten innerhalb der Kartenausdehnung des Kachelsatzes abzurufen. Der Kachelsatz kann in QGIS mit dem Tool XYZ-Kacheln (Verzeichnis) generieren erstellt werden.
Speichern Sie den Speicherort des Kachelsatzes. Er wird im nächsten Abschnitt verwendet.

Erstellen einer Höhenangaben-API mithilfe einer Azure Function und von RGB-codierten DEM-Kacheln aus Azure Blob Storage

Die RGB-codierten DEM-Kacheln müssen in einen Datenbankspeicher hochgeladen werden, bevor sie mit den Azure Functions zum Erstellen einer API verwendet werden können.

Laden Sie die Kacheln in Azure Blob Storage hoch. Azure Storage-Explorer ist hierfür ein nützliches Tool.

Das Hochladen von Kacheln in Azure Blob Storage kann einige Minuten dauern.
Sobald der Upload abgeschlossen ist, können Sie die Azure Function erstellen, um eine API zu erstellen, die die Höhe für eine bestimmte geografische Koordinate zurückgibt.

Diese Funktion empfängt ein Koordinatenpaar, bestimmt die Kachel, die es bei Zoomfaktor 14 abdeckt, und bestimmt dann die Pixelkoordinaten innerhalb dieser Kachel, die mit den geografischen Koordinaten übereinstimmen. Anschließend werden die Kachel und die RGB-Werte für dieses Pixel abgerufen und die folgende Formel verwendet, um die Höhe zu bestimmen:

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
        )

Führen Sie es lokal aus, um die Ergebnisse des Codebeispiels anzuzeigen:

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

Erstellen eines Konturlinienvektorkachel-Dienstes mithilfe von Azure Function und PostgreSQL

In diesem Abschnitt werden die Schritte beschrieben, um Konturlinien in QGIS zu erstellen und zu verarbeiten, sie in PostgreSQL hochzuladen und dann eine Azure Function zu erstellen, um PostgreSQL abzufragen, damit Vektorkacheln zurückgegeben werden.

Öffnen Sie in QGIS die zusammengeführten Rasterkacheln in der EPSG:4326-Projektion, die in Schritt 3 von Erstellen von Konturlinienvektorkacheln und RGB-codierten DEM-Kacheln erstellt wurde.
Wählen Sie im Menü Raster> die Option Extraktion – Kontur aus, um das Konturtool zu öffnen.

Wenn Sie Ausführen auswählen, werden Konturlinien erstellt und als Ebene zur Karte hinzugefügt. Einige der Konturlinienränder können etwas rau erscheinen. Dies wird im nächsten Schritt behoben.
Wählen Sie im Menü Verarbeitung die Option Toolbox aus, um die Verarbeitungstoolbox aufzurufen.
Wählen Sie dann im Abschnitt Vektorgeometrie der Verarbeitungstoolbox die Option Glatt aus.

Hinweis

Die Konturlinienglättung kann erheblich verbessert werden, indem eine größere Dateigröße in Kauf genommen wird.
Laden Sie die Konturlinien in die Datenbank. In dieser Anleitung wird die kostenlose Version der PostgreSQL-Datenbank verwendet, die auf localhost ausgeführt wird. Sie können sie auch in die Azure Database for PostgreSQL laden.

Der nächste Schritt erfordert eine PostgreSQL-Datenbank mit der PostGIS-Erweiterung.
Um eine Verbindung zwischen QGIS und PostgreSQL herzustellen, wählen Sie im Menü Ebene die Option >Ebene hinzufügen –PostGIS-Ebenen hinzufügen aus, und wählen Sie dann die Schaltfläche Neu aus.
Laden Sie als Nächstes Daten aus QGIS in PostgreSQL mit dem Datenbank-Manager in QGIS. Wählen Sie hierzu im Menü Datenbank die Option DB-Manager aus.
Stellen Sie eine Verbindung mit der PostGIS-Datenbank her, und wählen Sie Ebene/Datei importieren... aus, um Konturlinien in die Datenbank zu importieren.

Sie können jetzt eine Azure Function verwenden, um PostgreSQL abzufragen und Vektorkacheln für die Konturlinien zurückzugeben. Der Kachelserver kann mit dem Azure Maps-Web-SDK verwendet werden, um eine Web-App zu erstellen, die Konturlinien auf der Karte anzeigt.

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
        )