你当前正在访问 Microsoft Azure Global Edition 技术文档网站。如果需要访问由世纪互联运营的 Microsoft Azure 中国技术文档网站，请访问 https://docs.azure.cn。

创建提升数据和服务

2023-03-24

本指南介绍了如何使用 SRTM 任务中的 USGS 全球 DEM 数据以 30 米的准确度在 Microsoft Azure 云上生成高程服务。

本文介绍如何：

创建等高线矢量图块和 RGB 编码的 DEM 图块。
从 Azure Blob 存储使用 Azure 函数和 RGB 编码的 DEM 图块创建高程 API。
使用 Azure 函数和 PostgreSQL 创建等高线矢量图块服务。

先决条件

本指南要求使用以下第三方软件和数据：

USGS 数据。 DEM 数据可作为 GeoTiff 下载，每个图块通过 USGS EarthExplorer 进行 1 弧秒覆盖。这需要 EarthExplorer 帐户，但可以免费下载数据。
QGIS 桌面 GIS 应用程序用于处理光栅图块并使其平滑。可以免费下载和使用 QGIS。本指南使用 QGIS 版本 3.26.2-布宜诺斯艾利斯。
MapBox 开发的 rio-rgbify Python 包用于将 GeoTIFF 编码启用为 RGB。
具有 PostGIS 空间扩展的 PostgreSQL 数据库。

创建等高线矢量图块和 RGB 编码的 DEM 图块

本指南使用涵盖华盛顿州的 36 个图块，可从 USGS EarthExplorer 获取。

从 USGS EarthExplorer 下载光栅图块

搜索条件

选择要为其使用光栅图块的区域。出于演示目的，本指南使用“多边形”法在地图上选择区域。

导航到 USGS EarthExplorer。
在“搜索条件”选项卡中选择“多边形”，然后选择地图上的任何位置以创建边界。

数据集

选择“数据集”选项卡。
从“数字高程”部分选择“SRTM 1 弧秒全局”。

结果

选择“结果 >>”以查看所选区域和数据集的图块。
结果页面会显示可下载图块列表。要仅下载所需图块，请选择每个图块结果卡上的“下载选项”按钮，选择“GeoTIFF 1 弧秒”选项，然后对其余图块重复此步骤。
或者，使用批量下载选项并选择“GeoTIFF 1 弧秒”。

将光栅图块添加到 QGIS

获得所需的光栅图块后，可以在 QGIS 中导入它们。

将文件拖动到“QGIS 层”选项卡或选择“层”菜单中的“添加层”，从而将光栅图块添加到 QGIS。
将光栅层加载到 QGIS 时，可能会有不同的图块阴影。合并光栅层以解决此问题，这会形成 GeoTIFF 格式的单个平滑光栅图像。为此，请从“光栅”菜单中选择“杂项”，然后选择“合并...”
使用“将光栅层另存为”右键菜单访问“目录”->“导出”->“另存为”选项，将合并的光栅层重新投影到 EPSG:3857（WGS84/伪墨卡托投影）。需要 EPSG:3857 以将其与 Azure Maps Web SDK 结合使用。
如果只想创建等高线矢量图块，可以跳过以下步骤并转到“使用 Azure 函数和 PostgreSQL 创建等高线矢量图块服务”。

要创建高程 API，下一步是对 GeoTIFF 进行 RGB 编码。这可以使用 MapBox 开发的 rio-rgbify 来完成。直接在 Windows 中运行此工具存在一些挑战，因此从 WSL 运行更容易。下面是在 Ubuntu on 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/

显示 QGIS 中 RGB 编码的 GeoTIFF 的屏幕截图。

RGB 编码的 GeoTIFF 允许检索像素的 R、G 和 B 值，并从以下值计算高程：

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

接下来，创建图块集以与地图控件结合使用和/或用于获取图块集地图范围内任何地理坐标的高程。可以使用 生成 XYZ 图块（目录） 工具在 QGIS 中创建图块集。
保存图块集的位置，其将用于下一部分。

从 Azure Blob 存储使用 Azure 函数和 RGB 编码的 DEM 图块创建高程 API

需要先将 RGB 编码的 DEM 图块上传到数据库存储，然后再将其与 Azure Functions 结合使用以创建 API。

将图块上传到 Azure Blob 存储。 Azure 存储资源管理器是用于此目的的有用工具。

将图块上传到 Azure Blob 存储可能需要几分钟才能完成。
上传完成后，可以创建 Azure 函数以生成 API，该 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 中，在创建等高线矢量图块和 RGB 编码的 DEM 图块的步骤 3 中创建的 EPSG:4326 投影中打开合并的光栅图块。
从“光栅”菜单中选择“提取 - >等高线”，从而打开“等高线”工具。

选择“运行”以创建等高线并将其作为层添加到地图中。一些等高线边缘可能看起来有点粗糙。下一步将解决此问题。
从“处理”菜单中选择“工具箱”，从而打开“处理工具箱”。
然后，在“处理工具箱”的“矢量几何”部分，选择“平滑”。

注意

等高线平滑可以大幅改进，但代价是文件大小增加。
将等高线加载到数据库。本指南使用在 localhost 上运行的免费 PostgreSQL 数据库版本。还可以将它们加载到 Azure Database for PostgreSQL。

下一步需要具有 PostGIS 扩展的 PostgreSQL 数据库。
要创建从 QGIS 到 PostgreSQL 的连接，请从“层”菜单中选择“添加层”-“>添加 PostGIS 层”，然后选择“新建”按钮。
接下来，使用 QGIS 中的数据库管理器将数据从 QGIS 加载到 PostgreSQL。为此，请从“数据库”菜单中选择“数据库管理器”。
连接到 PostGIS 数据库，并选择“导入层/文件...”以将等高线导入数据库。

现在可以使用 Azure 函数查询 PostgreSQL 并返回等高线的矢量图块。图块服务器可与 Azure Maps Web SDK 结合使用，从而创建在地图上显示等高线的 Web 应用。

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
        )