Skapa utökade data och tjänster

Den här guiden beskriver hur du använder USGS globala DEM-data från deras SRTM-uppdrag med 30 m noggrannhet för att skapa en Elevation-tjänst i Microsoft Azure Cloud.

I den här artikeln beskrivs hur du:

  • Skapa konturlinjevektorpaneler och RGB-kodade DEM-paneler.
  • Skapa Elevation API med Azure Function och RGB-kodade DEM-paneler från Azure Blob Storage.
  • Skapa vektortiletjänst för konturlinjer med Azure Function och PostgreSQL.

Förutsättningar

Den här guiden kräver användning av följande programvara och data från tredje part:

  • USGS-data. DEM-data kan laddas ned som GeoTiff med 1 bågsekunders täckning per panel via USGS EarthExplorer. Detta kräver ett EarthExplorer-konto, men data kan laddas ned kostnadsfritt.
  • GIS-programmet för QGIS-skrivbordet används för att bearbeta och jämna ut Raster-panelerna. QGIS är kostnadsfritt att ladda ned och använda. Den här guiden använder QGIS version 3.26.2-Buenos Aires.
  • Python-paketet rio-rgbify, som utvecklats av MapBox, används för att koda GeoTIFF som RGB.
  • PostgreSQL-databas med PostGIS spatial extension.

Skapa konturlinjevektorpaneler och RGB-kodade DEM-paneler

Den här guiden använder de 36 panelerna som täcker delstaten Washington, som är tillgängliga från USGS EarthExplorer.

Ladda ned rasterpaneler från USGS EarthExplorer

Sökvillkor

Välj den region som du vill använda rasterpaneler för. I demonstrationssyfte använder den här guiden metoden "Polygon" för att välja regionen på kartan.

  1. Gå till USGS EarthExplorer.

  2. På fliken Sökvillkor väljer du Polygon och sedan någonstans på kartan för att skapa gränsen.

    En skärmbild som visar fliken sökvillkor på webbplatsen för USGS Earth Explorer.

Datauppsättningar

  1. Välj fliken Datauppsättningar .

  2. Välj SRTM 1 Arc-Second Global i avsnittet Digitala höjder.

    En skärmbild som visar fliken datauppsättningar på webbplatsen usgs earth explorer.

Resultat

  1. Välj Resultat >> för att visa panelerna för den valda regionen och datauppsättningen.

  2. Listan över nedladdningsbara paneler visas på resultatsidan. Om du bara vill ladda ned paneler väljer du knappen Ladda ned alternativ på resultatkortet för varje panel, väljer alternativet GeoTIFF 1 Arc-Second och upprepar det här steget för de återstående panelerna.

    En skärmbild som visar resultatfliken på webbplatsen usgs earth explorer.

  3. Du kan också använda alternativet massnedladdning och välja GeoTIFF 1 Arc-second.

Lägga till rasterpaneler i QGIS

När du har de rasterpaneler du behöver kan du importera dem i QGIS.

  1. Lägg till rasterpaneler i QGIS genom att dra filerna till fliken QGIS-lager eller välja Lägg till lager på menyn Lager .

    En skärmbild som visar rasterpaneler i QGIS.

  2. När rasterskikten läses in i QGIS kan det finnas olika nyanser av paneler. Åtgärda detta genom att sammanfoga rasterskikten, vilket resulterar i en enda jämn rasterbild i GeoTIFF-format. Det gör du genom att välja Diverse på Raster-menyn och sedan Sammanfoga...

    En skärmbild som visar menyn sammanfoga raster i QGIS.

  3. Återskapa det sammanfogade rasterskiktet till EPSG:3857 (WGS84/Pseudo-Mercator) med Spara Raster Layer som högerklicka för att komma åt det sammanfogade rasterskiktet i innehållstabellen -> Exportera ->Spara som-alternativet. EPSG:3857 krävs för att använda den med Azure Maps Web SDK.

    En skärmbild som visar hur sammanfoga rasterskikt-menyn i QGIS.

  4. Om du bara vill skapa konturlinjevektorpaneler kan du hoppa över följande steg och gå till Skapa konturlinjevektorpaneltjänst med Hjälp av Azure Function och PostgreSQL.

  5. För att skapa ett Elevation API är nästa steg att RGB-Koda GeoTIFF. Detta kan göras med hjälp av rio-rgbify, som utvecklats av MapBox. Det finns vissa utmaningar med att köra verktyget direkt i Windows, så det är enklare att köra från WSL. Nedan följer stegen i Ubuntu på 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/

    En skärmbild som visar RGB-kodad GeoTIFF i QGIS.

    Med RGB-kodade GeoTIFF kan du hämta R-, G- och B-värden för en pixel och beräkna höjningen från dessa värden:

    elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)

  6. Skapa sedan en uppsättning av plattor som ska användas med kartkontroll och/eller använd den för att hämta höjd för valfria geografiska koordinater inom kartomfattningen av uppsättningen av plattor. Paneluppsättningen kan skapas i QGIS med verktyget Generera XYZ-paneler (Katalog).

    En skärmbild som visar verktyget Generera XYZ-paneler (Katalog) i QGIS.

  7. Spara platsen för paneluppsättningen. Du använder den i nästa avsnitt.

Skapa Elevation API med Azure Function och RGB-kodade DEM-paneler från Azure Blob Storage

RGB-kodade DEM-paneler måste laddas upp till en databaslagring innan de kan användas med Azure Functions för att skapa ett API.

  1. Ladda upp panelerna till Azure Blob Storage. Azure Storage Explorer är ett användbart verktyg för detta ändamål.

    En skärmbild som visar Microsoft Azure Storage Explorer.

    Det kan ta flera minuter att ladda upp paneler till Azure Blob Storage.

  2. När uppladdningen är klar kan du skapa en Azure-funktion för att bygga ett API som returnerar höjd för en viss geografisk koordinat.

    Den här funktionen tar emot ett koordinatpar, fastställer panelen som täcker den på zoomnivå 14 och fastställer sedan pixelkoordinaterna i panelen som matchar de geografiska koordinaterna. Den hämtar sedan panelen, hämtar RGB-värdena för den pixeln och använder sedan följande formel för att fastställa höjningen:

    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
        )

Om du vill se resultatet av kodexemplet kör du det lokalt:

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

Skapa konturlinjevektorpaneltjänst med hjälp av Azure Function och PostgreSQL

I det här avsnittet beskrivs stegen för att skapa och bearbeta konturlinjer i QGIS, ladda upp dem till PostgreSQL och sedan skapa en Azure-funktion för att fråga PostgreSQL för att returnera vektorpaneler.

  1. I QGIS öppnar du de sammanfogade rasterpanelerna i projektionen EPSG:4326 som skapades i steg 3 i Skapa konturlinjevektorpaneler och RGB-kodade DEM-paneler.

  2. Välj Extrahering –> Kontur på Raster-menyn för att öppna konturverktyget.

    En skärmbild som visar konturdialogrutan i QGIS.

    Om du väljer Kör skapas konturlinjer och de läggs till som ett lager på kartan. några av konturlinjens kanter kan verka lite grova. Detta kommer att åtgärdas i nästa steg.

    En skärmbild som visar en karta med konturer i QGIS.

  3. Välj Verktygslåda på menyn Bearbetning för att visa verktygslådan För bearbetning.

  4. Välj Smooth (Utjämna) i avsnittet Vektorgeometri i verktygslådan Bearbetning.

    En skärmbild som visar den smidiga dialogrutan i QGIS.

    Anteckning

    Konturlinjeutjämning kan förbättras avsevärt men på bekostnad av ökad filstorlek.

  5. Läs in konturlinjerna till databasen. Den här guiden använder den kostnadsfria versionen av PostgreSQL-databasen som körs på localhost. Du kan också läsa in dem till Azure Database for PostgreSQL.

    Nästa steg kräver en PostgreSQL-databas med PostGIS-tillägget .

  6. Om du vill skapa en anslutning från QGIS till PostgreSQL väljer du Lägg till lager -Lägg till PostGIS-lager> på menyn Lager och väljer sedan knappen Nytt.

    En skärmbild som visar dialogrutan Skapa ny PostGIG-anslutning i QGIS.

  7. Läs sedan in data från QGIS till PostgreSQL med hjälp av Database Manager i QGIS. Det gör du genom att välja DB Manager på menyn Databas .

    En skärmbild som visar DB Manager i QGIS.

  8. Anslut till PostGIS-databasen och välj Importera lager/fil... för att importera konturlinjer till databasen.

    En skärmbild som visar dialogrutan importvektor i QGIS.

  9. Nu kan du använda en Azure-funktion för att fråga PostgreSQL och returnera vektorpaneler för konturlinjerna. Panelservern kan användas med Azure Maps Web SDK för att skapa en webbapp som visar konturlinjer på kartan.

    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
        )

Om du vill se resultatet av kodexemplet kör du det lokalt:

http://localhost:7071/api/tileserver?zoom={z}&x={x}&y={y}
  • Uppdaterades senast