Создание данных высот и услуг

2025-04-21

В этом руководстве описывается, как использовать глобальные данные DEM от миссии SRTM USGS с точностью 30м для создания высотной службы на платформе Microsoft Azure Cloud.

В этой статье описывается, как выполнить следующее:

Создайте векторные плитки контурной линии и плитки DEM в кодировке RGB.
Создание API высот с использованием функции Azure и плиток ЦМР, закодированных в RGB, из хранилища блобов Azure.
Создайте векторный тайловый сервис линий контура, используя Azure Functions и PostgreSQL.

Предварительные условия

В этом руководстве требуется использовать следующее стороннее программное обеспечение и данные:

Данные USGS. Данные DEM можно скачать как GeoTiff с разрешением 1 угловая секунда на тайл через USGS EarthExplorer. Для этого требуется учетная запись EarthExplorer, но данные можно скачать бесплатно.
Настольное приложение QGIS используется для обработки и сглаживания растровых плиток. QGIS предоставляется бесплатно для скачивания и использования. В этом руководстве используется QGIS версии 3.26.2-Буэнос-Айрес.
Пакет Python rio-rgbify, разработанный MapBox, используется для кодирования GeoTIFF как RGB.
База данных PostgreSQL с пространственным расширением PostGIS .

Создание векторных плиток линии контура и плитки DEM в кодировке RGB

В этом руководстве используются 36 плиток, покрывающих штат Вашингтон, доступные через USGS EarthExplorer.

Скачивание растровых плиток из USGS EarthExplorer

Критерии поиска

Выберите регион, для которого требуется растровые плитки. Для демонстрационных целей в этом руководстве используется метод Polygon для выбора региона на карте.

Перейдите к USGS EarthExplorer.
На вкладке "Критерии поиска" выберите Polygon и выберите в любом месте карты, чтобы создать границу.

Наборы данных

Перейдите на вкладку "Наборы данных".
Выберите SRTM 1 Arc-Second Global в разделе "Цифровые высоты".

Результаты

Выберите "Результаты >> ", чтобы просмотреть плитки для выбранного региона и набора данных.
Список загружаемых плиток отображается на странице результатов. Чтобы скачать только нужные плитки, нажмите кнопку "Параметры загрузки" на карточке результатов для каждой плитки, выбрав параметр GeoTIFF 1 Arc-Second и повторите этот шаг для оставшихся плиток.
Кроме того, используйте параметр массового скачивания и выберите GeoTIFF 1 Arc-second.

Добавление растровых плиток в QGIS

После получения нужных растровых плиток их можно импортировать в QGIS.

Добавьте растровые плитки в QGIS, перетащив файлы на вкладку уровня QGIS или выбрав "Добавить слой" в меню "Слой".
Когда слои растра загружаются в QGIS, могут быть разные оттенки плиток. Исправьте это путем объединения слоев растров, что приводит к одному гладкому растровому изображению в формате GeoTIFF. Для этого выберите Прочее из меню Растр, а затем Слияние...
Перепроектируйте объединенный слой растра в EPSG:3857 (WGS84 / Псевдо-Меркатор) с помощью Сохранить растровый слой как, щелкнув правой кнопкой мыши, чтобы получить доступ к объединенному слою растра в таблице содержимого -Экспорт -Сохранить как. EPSG:3857 требуется для использования с Azure Maps Web SDK.
Если вы хотите создать только векторы контурной линии, можно пропустить следующие шаги и перейти к службе создания векторной плитки контурной линии с помощью функции Azure и PostgreSQL.
Чтобы создать API высоты, следующим шагом необходимо выполнить RGB-кодирование GeoTIFF. Это можно сделать с помощью rio-rgbify, разработанного MapBox. Есть несколько проблем, связанных с выполнением этого средства непосредственно в Windows, поэтому проще запускать его из WSL. Ниже приведены шаги в Ubuntu на 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/
```
В RGB-кодированном GeoTIFF можно получить значения R, G и B пикселя и вычислить высоту из этих значений.

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Затем создайте набор плиток, используемый с элементом управления картой и (или) используйте его для получения высоты для любых географических координат в пределах области карты набора плиток. Набор плиток можно создать в QGIS с помощью средства Создания плиток XYZ (Directory).
Сохраните расположение набора плиток, используйте его в следующем разделе.

Создание API для определения высоты с использованием функции Azure и плит DEM, закодированных в RGB, из хранилища BLOB-объектов Azure.

Плитки DEM в кодировке RGB необходимо передать в хранилище базы данных, прежде чем их можно будет использовать с Azure Functions для создания API.

Отправьте плитки в Хранилище BLOB-объектов Azure. служба хранилища Azure Explorer — это полезное средство для этой цели.

Отправка плиток в Хранилище BLOB-объектов Azure может занять несколько минут.
После завершения загрузки вы можете создать функцию Azure для построения API, возвращающего высоту для заданной географической координаты.

Эта функция получает пару координат, определяет плитку, которая охватывает ее на уровне масштабирования 14, а затем определяет координаты пикселей внутри этой плитки, которая соответствует географическим координатам. Затем он извлекает плитку, получает значения RGB для этого пикселя, а затем использует следующую формулу для определения повышения:

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
        )

Чтобы просмотреть результаты примера кода, запустите его локально:

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

Создать службу векторных тайлов линий контура с помощью функции Azure и PostgreSQL

В этом разделе описаны шаги по созданию и обработке линий контура в QGIS, их отправке в PostgreSQL, а затем создание функции Azure для запроса PostgreSQL для возврата векторных плиток.

В QGIS откройте объединенные растровые плитки в проекции EPSG:4326, созданные на шаге 3 задачи Создание векторных плиток с контурными линиями и плиток ЦМР с кодировкой RGB.
Выберите извлечение —> контур из меню Растра, чтобы открыть средство "Контур".

При выборе Выполнить создаются контурные линии и добавляются как слой на карту. некоторые края линии контура могут показаться немного грубыми. Это будет решено на следующем шаге.
Выберите Инструменты в меню Обработка, чтобы открыть Инструменты обработки.
Затем выберите "Smooth" в разделе Геометрия вектораИнструменты обработки.

Примечание.

Сглаживание линии контура может быть значительно улучшено, но за счет увеличения размера файла.
Загрузите линии контура в базу данных. В этом руководстве используется бесплатная версия базы данных PostgreSQL , которая выполняется в localhost. Их также можно загрузить в База данных Azure для PostgreSQL.

Для следующего шага требуется база данных PostgreSQL с расширением PostGIS .
Чтобы создать подключение из QGIS к PostgreSQL, выберите Add Layer ->Добавить слои PostGIS из меню Слой, а затем нажмите кнопку "Создать".
Затем загрузите данные из QGIS в PostgreSQL с помощью диспетчера баз данных в QGIS. Для этого выберите db Manager в меню "База данных ".
Подключитесь к базе данных PostGIS и выберите "Импорт слоя или файла", чтобы импортировать линии контура в базу данных.

Теперь можно использовать функцию Azure для запроса PostgreSQL и возврата векторных плиток для линий контура. Сервер плиток можно использовать с веб-пакетом SDK Azure Maps для создания веб-приложения, отображающего контурные линии на карте.

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
        )