Créer des services et des données d’élévation

2025-04-30

Ce guide explique comment utiliser les données DEM mondiales de l’USGS à partir de sa mission SRTM avec une précision de 30 m pour créer un service d’élévation sur le cloud Microsoft Azure.

Cet article explique comment :

Créer des vignettes vectorielles de ligne de contour et des vignettes DEM encodées en RVB.
Créer une API d'élévation en utilisant une fonction Azure et des tuiles DEM encodées en RVB depuis Azure Blob Storage.
Créer un service de vignettes vectorielles de ligne de contour à l’aide d’une fonction Azure et de PostgreSQL.

Prérequis

Ce guide nécessite l’utilisation des données et logiciels tiers suivants :

Données USGS. Les données DEM peuvent être téléchargées au format GeoTiff avec une couverture d’1 seconde d’arc par vignette par le biais d’USGS EarthExplorer. Cela nécessite un compte EarthExplorer, mais les données peuvent être téléchargées gratuitement.
L’application GIS de bureau QGIS est utilisée pour traiter et lisser les vignettes raster. Vous pouvez télécharger et utiliser gratuitement QGIS. Ce guide utilise QGIS version 3.26.2 Buenos Aires.
Le package Python rio-rgbify, développé par MapBox, est utilisé pour encoder le format GeoTIFF au format RVB.
Base de données PostgreSQL avec l’extension spatiale PostGIS.

Créer des vignettes vectorielles de ligne de contour et des vignettes DEM encodées en RVB

Ce guide utilise les 36 vignettes couvrant l’État de Washington, disponibles auprès d’USGS EarthExplorer.

Télécharger les vignettes raster à partir d’USGS EarthExplorer

Critères de recherche

Sélectionnez la région dont vous souhaitez obtenir les vignettes raster. À des fins de démonstration, ce guide utilise la méthode « Polygon » (Polygone) pour sélectionner la région sur la carte.

Accédez à USGS EarthExplorer.
Dans l’onglet Critères de recherche, sélectionnez Polygone, puis sélectionnez n’importe où sur la carte pour créer la limite.

Jeux de données

Sélectionnez l’onglet Jeux de données.
Sélectionnez SRTM 1 Arc-Second Global (SRTM mondial 1 seconde d’arc) dans la section Digital Elevations (Élévations numériques).

Résultats

Sélectionnez Results>> (Résultats) pour afficher les vignettes de la région et du jeu de données sélectionnés.
La liste des vignettes téléchargeables s’affiche sur la page de résultats. Pour télécharger les vignettes souhaitées uniquement, sélectionnez le bouton Download Options (Options de téléchargement) dans la fiche de résultat de chaque vignette, en sélectionnant l’option GeoTIFF 1 Arc-Second (GeoTIFF 1 seconde d’arc), puis répétez cette étape pour les vignettes restantes.
Vous pouvez également utiliser l’option de téléchargement en bloc et sélectionner GeoTIFF 1 Arc-Second (GeoTIFF 1 seconde d’arc).

Ajouter des vignettes raster à QGIS

Une fois que vous disposez des vignettes raster dont vous avez besoin, vous pouvez les importer dans QGIS.

Ajoutez les vignettes raster à QGIS en faisant glisser les fichiers vers l’onglet QGIS layer (Couche QGIS) ou en sélectionnant Add Layer (Ajouter une couche) dans le menu Couche.
Lorsque les couches raster sont chargées dans QGIS, il peut y avoir différentes nuances de vignettes. Corrigez ce problème en fusionnant les couches raster, ce qui entraîne une image raster lisse unique au format GeoTIFF. Pour ce faire, sélectionnez Miscellaneous dans le menu Raster, puis Merge....
Projetez à nouveau la couche raster fusionnée sur EPSG:3857 (WGS84 / Pseudo-Mercator) en faisant un clic droit sur Enregistrer la couche raster comme pour accéder à la couche raster fusionnée dans l’option table des matières ->Exporter ->Enregistrer sous. EPSG:3857 est requis pour utiliser le SDK Web Azure Maps.
Si vous souhaitez créer uniquement des vignettes vectorielles de ligne de contour, vous pouvez ignorer les étapes suivantes et accéder à l’étape Créer un service de vignettes vectorielles de ligne de contour à l’aide d’une fonction Azure et de PostgreSQL.
Pour créer une API Elevation, l’étape suivante consiste à encoder en RVB le format GeoTIFF. Pour ce faire, vous pouvez utiliser rio-rgbify, développé par MapBox. L’exécution de cet outil directement dans Windows pose quelques problèmes ; il est dès lors plus facile de l’utiliser à partir de WSL. Voici la procédure à suivre dans Ubuntu sur 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/
```
Le format GeoTIFF encodé en RVB vous permet de récupérer les valeurs R, V et B d’un pixel et de calculer l’élévation à partir de ces valeurs :

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Ensuite, créez un jeu de vignettes à utiliser avec le contrôle de carte et/ou utilisez-le pour obtenir l’élévation pour n’importe quelles coordonnées géographiques dans l’étendue de la carte du jeu de vignettes. Le jeu de vignettes peut être créé dans QGIS à l’aide de l’outil Generate XYZ tiles (Directory) (Générer des vignettes XYZ (répertoire)).
Enregistrez l’emplacement du jeu de vignettes. Vous l’utiliserez dans la section suivante.

Créer une API Elevation à l’aide d’une fonction Azure et de vignettes DEM encodées en RVB à partir de Stockage Blob Azure

Vous devez charger les vignettes DEM encodées en RVB dans un stockage de base de données avant de pouvoir utiliser ces vignettes avec les fonctions Azure pour créer une API.

Chargez les vignettes dans Stockage Blob Azure. Pour cela, l’Explorateur Stockage Azure est un outil utile.

Le chargement de tuiles dans Stockage Blob Azure peut prendre plusieurs minutes.
Une fois le chargement terminé, vous pouvez créer une fonction Azure pour générer une API qui retourne l’élévation de coordonnées géographiques données.

Cette fonction reçoit une paire de coordonnées, détermine la vignette qui la couvre au niveau de zoom 14, puis détermine les coordonnées du pixel dans cette vignette qui correspondent aux coordonnées géographiques. Elle récupère ensuite la vignette, obtient les valeurs RVB de ce pixel, puis utilise la formule suivante pour déterminer l’élévation :

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
        )

Pour afficher les résultats de l’exemple de code, exécutez-le localement :

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

Créer un service de vignettes vectorielles de ligne de contour à l’aide d’une fonction Azure et de PostgreSQL

Cette section décrit la procédure à suivre pour créer et traiter des lignes de contour dans QGIS, les charger dans PostgreSQL, puis créer une fonction Azure pour interroger PostgreSQL pour retourner des vignettes vectorielles.

Dans QGIS, ouvrez les vignettes raster fusionnées dans la projection EPSG:4326 créée à l’étape 3 de la section Créer des vignettes vectorielles de ligne de contour et des vignettes DEM encodées en RVB.
Sélectionnez Extraction -> Contour dans le menu Raster pour ouvrir l’outil Contour.

Le bouton Run (Exécuter) permet de créer des lignes de contour et de les ajouter en tant que couche à la carte. Certains bords de ligne de contour peuvent sembler un peu grossiers. Ce problème sera corrigé à l’étape suivante.
Sélectionnez Toolbox (Boîte à outils) dans le menu Processing (Traitement) pour afficher la boîte à outils de traitements.
Sélectionnez ensuite Smooth (Lisser) dans la section Vector geometry (Géométrie vectorielle) de la boîte à outils de traitements.

Remarque

Le lissage des lignes de contour peut être considérablement amélioré, mais au détriment d’une taille de fichier accrue.
Chargez les lignes de contour dans la base de données. Ce guide utilise la version gratuite de la base de données PostgreSQL qui s’exécute sur localhost. Vous pouvez également les charger dans Azure Database pour PostgreSQL.

L’étape suivante nécessite une base de données PostgreSQL avec l’extension PostGIS.
Pour créer une connexion entre QGIS et PostgreSQL, sélectionnez Add Layer (Ajouter une couche) ->Add PostGIS Layers (Ajouter des couches PostGIS) dans le menu Layer (Couche), puis sélectionnez le bouton New (Nouveau).
Ensuite, chargez les données de QGIS vers PostgreSQL à l’aide du gestionnaire de base de données dans QGIS. Pour cela, sélectionnez DB Manager dans le menu Database.
Connectez-vous à la base de données PostGIS, puis sélectionnez Import Layer/File (Import de couche/fichier) pour importer les lignes de contour dans la base de données.

Vous pouvez à présent utiliser une fonction Azure pour interroger PostgreSQL et retourner des vignettes vectorielles pour les lignes de contour. Le serveur de vignettes peut être utilisé avec le Kit de développement logiciel (SDK) web Azure Maps pour créer une application web qui affiche des lignes de contour sur la carte.

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
        )