Tworzenie danych i usług podniesienia uprawnień

2025-04-21

W tym przewodniku opisano, jak używać danych USGS na całym świecie DEM z ich misji SRTM z dokładnością 30m w celu utworzenia usługi podniesienia uprawnień w chmurze Microsoft Azure.

W tym artykule opisano sposób wykonywania następujących czynności:

Tworzenie kafelków wektorów liniowych konturu i kafelków DEM zakodowanych w formacie RGB.
Tworzenie interfejsu API elewacji za pomocą funkcji Azure i kafelków DEM zakodowanych w formacie RGB z usługi Azure Blob Storage.
Utwórz usługę kafli wektorowych linii konturowych, używając Funkcji Azure i bazy danych PostgreSQL.

Wymagania wstępne

Ten przewodnik wymaga użycia następującego oprogramowania i danych innych firm:

Dane USGS. Dane DEM można pobrać jako GeoTiff z pokryciem z dokładnością do 1 sekundy kątowej na kafelkę za pośrednictwem USGS EarthExplorer. Wymaga to konta EarthExplorer, ale dane można pobrać bezpłatnie.
Aplikacja QGIS desktop GIS służy do przetwarzania i wygładzania kafelków Raster. QGIS jest bezpłatny do pobrania i użycia. W tym przewodniku jest używana wersja QGIS 3.26.2-Buenos Aires.
Pakiet rio-rgbify języka Python, opracowany przez usługę MapBox, służy do kodowania geoTIFF jako RGB.
Baza danych PostgreSQL z rozszerzeniem przestrzennym PostGIS .

Tworzenie kafelków wektorów liniowych konturu i kafelków DEM zakodowanych w formacie RGB

W tym przewodniku użyto 36 kafelków obejmujących stan Waszyngton, dostępny z USGS EarthExplorer.

Pobieranie kafelków rasterowych z USGS EarthExplorer

Kryteria wyszukiwania

Wybierz region, dla którego chcesz umieścić kafelki rastrowe. W celach demonstracyjnych w tym przewodniku użyto metody "Polygon" do wybrania regionu na mapie.

Przejdź do USGS EarthExplorer.
Na karcie Kryteria wyszukiwania wybierz pozycję Wielokąt , a następnie wybierz dowolne miejsce na mapie, aby utworzyć granicę.

Zestawy danych

Wybierz kartę Zestawy danych.
Wybierz pozycję SRTM 1 Arc-Second Global w sekcji Cyfrowe wzniesienia.

Wyniki

Wybierz Wyniki >>, aby wyświetlić kafelki dla wybranego regionu i zestawu danych.
Lista kafelków możliwych do pobrania jest wyświetlana na stronie wyników. Aby pobrać tylko żądane kafelki, wybierz przycisk Opcje pobierania na karcie wyników dla każdego kafelka, wybierając opcję GeoTIFF 1 Arc-Second i powtórz ten krok dla pozostałych kafelków.
Alternatywnie użyj opcji pobierania zbiorczego i wybierz pozycję GeoTIFF 1 Arc-second.

Dodawanie kafelków rastrowych do zestawu QGIS

Gdy masz potrzebne kafelki rastrowe, możesz je zaimportować w systemie QGIS.

Dodaj kafelki rastrowe do zestawu QGIS, przeciągając pliki na kartę warstwy QGIS lub wybierając polecenie Dodaj warstwę w menu Warstwa .
Gdy warstwy rastrowe są ładowane do QGIS, mogą istnieć różne odcienie kafelków. Napraw to, scalając warstwy rastrowe, co powoduje utworzenie pojedynczego gładkiego obrazu rasterowego w formacie GeoTIFF. W tym celu wybierz Różne z menu Raster, a następnie Scal...
Przeprojektuj scaloną warstwę rastrową do EPSG:3857 (WGS84 / Pseudo-Mercator) przy użyciu opcji Zapisz warstwę rastrową jako dostępnej po kliknięciu prawym przyciskiem myszy na scaloną warstwę w tabeli zawartości - Eksportuj - Zapisz jako. EpsG:3857 jest wymagany do używania go z zestawem Sdk sieci Web usługi Azure Maps.
Jeśli chcesz tylko utworzyć kafelki wektorowe linii konturowej, możesz pominąć poniższe kroki i przejść do Tworzenie Usługi Kafelków Wektorowych Linii Konturowej z użyciem funkcji Azure i PostgreSQL.
Aby utworzyć API wysokościowe, następnym krokiem jest zakodowanie GeoTIFF przy użyciu RGB. Można to zrobić przy użyciu rio-rgbify, opracowanego przez rozwiązanie MapBox. Istnieją pewne wyzwania związane z uruchamianiem tego narzędzia bezpośrednio w systemie Windows, więc łatwiej jest uruchomić z poziomu programu WSL. Poniżej znajdują się kroki w Ubuntu na 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/
```
Zakodowany w formacie RGB geoTIFF umożliwia pobieranie wartości R, G i B dla piksela i obliczanie podniesienia z następujących wartości:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Następnie utwórz zestaw kafelków do użycia z kontrolą mapy i/lub użyj go, aby uzyskać wysokość terenu dla wszelkich współrzędnych geograficznych w obrębie obszaru mapy zestawu kafelków. Zestaw kafelków można utworzyć w QGIS przy użyciu narzędzia Generowanie kafelków XYZ (katalog).
Zapisz lokalizację zestawu kafelków, użyjesz jej w następnej sekcji.

Tworzenie interfejsu API elewacji przy użyciu funkcji Azure i kafelków DEM zakodowanych w RGB z usługi Azure Blob Storage

Aby można było utworzyć interfejs API za pomocą funkcji Azure, należy przekazać kafelki DEM zakodowane w formacie RGB do magazynu bazy danych.

Przekaż kafelki do usługi Azure Blob Storage. Eksplorator usługi Azure Storage jest przydatnym narzędziem do tego celu.

Przekazywanie kafelków do usługi Azure Blob Storage może potrwać kilka minut.
Po zakończeniu przekazywania możesz utworzyć Funkcję Azure, aby zbudować interfejs API zwracający wysokość dla danej współrzędnej geograficznej.

Ta funkcja otrzymuje parę współrzędnych, określ kafelek, który obejmuje go na poziomie powiększenia 14, a następnie określ współrzędne pikseli w tym kafelku, które pasują do współrzędnych geograficznych. Następnie pobiera płytkę, uzyskuje wartości RGB dla piksela, a potem używa następującej formuły do określenia wysokości.

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
        )

Aby wyświetlić wyniki przykładowego kodu, uruchom go lokalnie:

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

Stwórz usługę kafli wektorowych z liniami konturowymi przy użyciu funkcji Azure i PostgreSQL

W tej sekcji opisano kroki tworzenia i przetwarzania linii konturowych w QGIS, przesyłania ich do bazy danych PostgreSQL, a następnie tworzenia funkcji Azure do wykonywania zapytań w PostgreSQL w celu zwrócenia płytek wektorowych.

W QGIS otwórz scalone kafelki rastrowe w projekcji EPSG:4326, utworzone w kroku 3 z Utwórz wektorowe kafelki linii konturu i kafelki DEM zakodowane w RGB.
Wybierz pozycję Wyodrębnianie —> Kontur z menu Raster, aby otworzyć narzędzie Kontur.

Wybranie pozycji Uruchom powoduje utworzenie linii konturowych i dodanie ich jako warstwy do mapy. niektóre krawędzie linii konturowej mogą wydawać się nieco szorstkie. Zostanie to rozwiązane w następnym kroku.
Wybierz Przybornik z menu Przetwarzanie, aby wyświetlić Przybornik Przetwarzania.
Następnie wybierz Wygładzanie w sekcji Geometria Wektorowa w Przyborniku Przetwarzania.

Uwaga

Wygładzenie linii konturowej można znacznie poprawić, ale kosztem zwiększonego rozmiaru pliku.
Załaduj linie konturowe do bazy danych. W tym przewodniku jest używana bezpłatna wersja bazy danych PostgreSQL uruchomiona na hoście lokalnym. Można je również załadować do usługi Azure Database for PostgreSQL.

Następny krok wymaga bazy danych PostgreSQL z rozszerzeniem PostGIS .
Aby utworzyć połączenie z platformy QGIS do bazy danych PostgreSQL, wybierz pozycję Dodaj warstwę ->Dodaj warstwy PostGIS z menu Warstwa , a następnie wybierz przycisk Nowy .
Następnie załaduj dane z zestawu QGIS do bazy danych PostgreSQL przy użyciu Menedżera baz danych w systemie QGIS. W tym celu wybierz pozycję Menedżer bazy danych z menu Baza danych .
Połącz się z bazą danych PostGIS i wybierz pozycję Importuj warstwę/plik... w celu zaimportowania linii konturowych do bazy danych.

Teraz możesz użyć funkcji platformy Azure do wykonywania zapytań z bazy danych PostgreSQL i zwracania kafelków wektorowych dla linii konturowych. Serwer kafelków może służyć do tworzenia aplikacji internetowej usługi Azure Maps, która wyświetla linie konturowe na mapie.

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
        )