Hoogtegegevens en services maken

In deze handleiding wordt beschreven hoe u USGS worldwide DEM-gegevens van hun SRTM-missie gebruikt met een nauwkeurigheid van 30 meter om een Elevation-service te bouwen in de Microsoft Azure Cloud.

In dit artikel komt het volgende aan de orde:

Maak contourlijnvectortegels en RGB-gecodeerde DEM-tegels.
Maak elevation-API met behulp van Azure Function en RGB-gecodeerde DEM-tegels van Azure Blob Storage.
Een contourlijnvectortegelservice maken met behulp van Azure Function en PostgreSQL.

Vereisten

Voor deze handleiding is het gebruik van de volgende software en gegevens van derden vereist:

USGS-gegevens. DEM-gegevens kunnen worden gedownload als GeoTiff met een dekking van 1 boog seconde per tegel via de USGS EarthExplorer. Hiervoor is een EarthExplorer-account vereist, maar de gegevens kunnen gratis worden gedownload.
De desktop GIS-toepassing QGIS wordt gebruikt om de rastertegels te verwerken en te versoepelen. QGIS is gratis te downloaden en te gebruiken. Deze gids gebruikt QGIS versie 3.26.2-Buenos Aires.
Het Python-pakket rio-rgbify , ontwikkeld door MapBox, wordt gebruikt om geoTIFF als RGB te coderen.
PostgreSQL-database met de ruimtelijke PostGIS-extensie .

Maak contourlijn-, vector- en RGB-gecodeerde DEM-tegels

In deze handleiding worden de 36 tegels gebruikt die betrekking hebben op de staat Washington, beschikbaar via USGS EarthExplorer.

Rastertegels downloaden van USGS EarthExplorer

Zoekcriteria

Selecteer de regio waarvoor u rastertegels wilt gebruiken. Voor demonstratiedoeleinden gebruikt deze handleiding de methode 'Veelhoek' om de regio op de kaart te selecteren.

Navigeer naar de USGS EarthExplorer.
Selecteer op het tabblad Zoekcriteria de optie Veelhoek en selecteer vervolgens een willekeurige plaats op de kaart om de grens te maken.

Gegevenssets

Selecteer het tabblad Gegevenssets .
Selecteer SRTM 1 Arc-Second Global in de sectie Digitale hoogtes.

Resultaten

Selecteer Resultaten >> om de tegels voor de geselecteerde regio en gegevensset weer te geven.
De lijst met downloadbare tegels wordt weergegeven op de resultatenpagina. Als u alleen de gewenste tegels wilt downloaden, selecteert u de knop Downloadopties op de resultatenkaart voor elke tegel, waarbij u de optie GeoTIFF 1 Arc-Second selecteert en deze stap herhaalt voor de resterende tegels.
U kunt ook de optie voor bulksgewijs downloaden gebruiken en GeoTIFF 1 Arc-second selecteren.

Rastertegels toevoegen aan QGIS

Zodra u de rastertegels hebt die u nodig hebt, kunt u ze importeren in QGIS.

Voeg rastertegels toe aan QGIS door de bestanden naar het tabblad QGIS-laag te slepen of Laag toevoegen te selecteren in het menu Laag.
Wanneer de rasterlagen in QGIS worden geladen, kunnen er verschillende tinten tegels zijn. Los dit op door de rasterlagen samen te voegen, wat resulteert in één vloeiende rasterafbeelding in geoTIFF-indeling. Hiervoor selecteert u Diversen in het menu Raster en vervolgens Samenvoegen...
Herprojecteer de samengevoegde rasterlaag naar EPSG:3857 (WGS84 / Pseudo-Mercator) door met de rechtermuisknop te klikken op Rasterlaag opslaan als om toegang te krijgen tot de samengevoegde rasterlaag in de tabel met inhoud -> > -> Opslaan als. EPSG:3857 is vereist voor gebruik met azure Maps Web SDK.
Als u alleen contourlijnvectortegels wilt maken, kunt u de volgende stappen overslaan en gaan naar Contourlijnvector-tegelservice maken met behulp van Azure Function en PostgreSQL.
Voor het maken van een Elevation API is de volgende stap om de GeoTIFF te RGB-coderen. Dit kan worden gedaan met behulp van rio-rgbify, ontwikkeld door MapBox. Er zijn enkele uitdagingen bij het uitvoeren van dit hulpprogramma rechtstreeks in Windows, dus is het eenvoudiger dit vanuit WSL te doen. Hieronder ziet u de stappen in Ubuntu in 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/
```
Met de RGB-gecodeerde GeoTIFF kunt u R-, G- en B-waarden voor een pixel ophalen en de hoogte van deze waarden berekenen:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Maak vervolgens een tegelset die moet worden gebruikt met de kaartcontrole en/of gebruik deze om de hoogte op te halen voor geografische coördinaten binnen het kaartbereik van de tegelset. De tegelset kan worden gemaakt in QGIS met behulp van het hulpprogramma XYZ-tegels genereren (Directory).
Sla de locatie van de tegelset op. U gebruikt deze in de volgende sectie.

Hoogte-API maken met behulp van Azure Function en RGB-gecodeerde DEM-tegels van de Azure Blob Storage

De RGB-gecodeerde DEM-tegels moeten worden geüpload naar een databaseopslag voordat deze kan worden gebruikt met de Azure Functions om een API te maken.

Upload de tegels naar Azure Blob Storage. Azure Storage Explorer is een nuttig hulpprogramma voor dit doel.

Het uploaden van tegels naar Azure Blob Storage kan enkele minuten duren.
Zodra het uploaden is voltooid, kunt u een Azure-functie maken om een API te bouwen die de hoogte voor een bepaalde geografische coördinaat retourneert.

Deze functie ontvangt een coördinaatpaar, bepaalt de tegel die deze bedekt op zoomniveau 14 en bepaalt vervolgens de pixelcoördinaten binnen die tegel die overeenkomt met de geografische coördinaten. Vervolgens wordt de tegel opgehaald, worden de RGB-waarden voor die pixel opgehaald en wordt vervolgens de volgende formule gebruikt om de verhoging te bepalen:

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
        )

Als u de resultaten van het codevoorbeeld wilt zien, voert u het lokaal uit:

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

Een contourlijnvectortegelservice maken met behulp van Azure Function en PostgreSQL

In deze sectie worden de stappen beschreven voor het maken en verwerken van contourlijnen in QGIS, upload ze naar PostgreSQL en maak vervolgens een Azure-functie om Query PostgreSQL om vectortegels te retourneren.

Open in QGIS de samengevoegde rastertegels in de EPSG:4326-projectie die is gemaakt in stap 3 van Contour lijnvectortegels en RGB-gecodeerde DEM-tegels maken.
Selecteer Extractie -> Contour in het menu Raster om het hulpmiddel Contour te openen.

Als u Uitvoeren selecteert, worden contourlijnen gemaakt en als laag aan de kaart toegevoegd. sommige van de contourlijnranden kunnen een beetje ruw lijken. Dit wordt in de volgende stap opgelost.
Selecteer Toolbox in het Verwerkingsmenu om de Verwerkingstoolbox weer te geven.
Vervolgens selecteer Smooth in de sectie Vectorgeometrie van de Verwerkingstoolbox.

Notitie

Het gladmaken van contourlijnen kan aanzienlijk worden verbeterd, maar ten koste van een grotere bestandsgrootte.
Laad de contourlijnen naar de database. In deze handleiding wordt de gratis versie van postgreSQL-database gebruikt die wordt uitgevoerd op localhost. U kunt ze ook laden in Azure Database for PostgreSQL.

Voor de volgende stap is een PostgreSQL-database met postGIS-extensie vereist.
Als u een verbinding wilt maken van QGIS naar PostgreSQL, selecteert u Laag toevoegen -> toevoegen in het menu Laag en selecteert u vervolgens de knop Nieuw.
Laad vervolgens gegevens van QGIS naar PostgreSQL met behulp van de Database Manager in QGIS. Hiervoor selecteert u DB Manager in het menu Database .
Maak verbinding met de PostGIS-database en selecteer Import Layer/File... om contourlijnen naar de database te importeren.

U kunt nu een Azure-functie gebruiken om Query's uit te voeren op PostgreSQL en vectortegels te retourneren voor de contourlijnen. De tegelserver kan worden gebruikt met de Azure Maps Web SDK om een web-app te maken die contourlijnen op de kaart weergeeft.

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
        )