標高データとサービスを作成する

このガイドでは、SRTM ミッションの USGS ワールドワイド DEM データを 30 m の精度で使用して、Microsoft Azure Cloud に海抜サービスを構築する方法について説明します。

この記事では、次の方法について説明します。

等高線ベクタータイルとRGBでエンコードされたデジタル標高モデル（DEM）タイルを作成します。
Azure Function と RGB エンコードされた DEM タイルを使用して、Azure Blob Storage から Elevation API を作成します。
Azure 関数と PostgreSQL を使用して、等高線ベクタータイルサービスを作成します。

前提条件

このガイドでは、次のサードパーティ製ソフトウェアとデータを使用する必要があります。

USGS データ。 DEM データは、USGS EarthExplorer を介してタイルごとに 1 秒角のカバレッジを持つ GeoTiff としてダウンロードできます。これには EarthExplorer アカウントが必要ですが、データは無料でダウンロードできます。
QGIS デスクトップ GIS アプリケーションは、ラスタータイルを処理して滑らかにするために使用されます。 QGIS は無料でダウンロードして使用できます。このガイドでは QGIS バージョン 3.26.2-Buenos Aires を使用しています。
MapBox によって開発された rio-rgbify Python パッケージは、GeoTIFF を RGB としてエンコードするために使用されます。
PostGIS 空間拡張機能を持つ PostgreSQL データベース。

等高線ベクタータイルとRGBエンコードされたデジタル標高モデル(Digital Elevation Model)タイルを作成する

このガイドでは、USGS EarthExplorer から入手できる、ワシントン州をカバーする 36 タイルを使用します。

USGS EarthExplorer からラスタータイルをダウンロードする

検索条件

ラスタータイルの対象となる領域を選択します。デモする目的で、このガイドでは "Polygon" メソッドを使用してマップ上の領域を選択します。

USGS EarthExplorer に移動します。
[検索条件] タブで [ポリゴン] を選択し、マップの任意の場所を選択して境界を作成します。

データセット

[データセット] タブを選びます。
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 を使用して等高線ベクタータイルサービスを作成するに進むことができます。
Elevation API を作成するには、次の手順として GeoTIFF を RGB エンコードします。これは、MapBox によって開発された rio-rgbify を使用して行うことができます。 Windows でこのツールを直接実行する場合、いくつかの課題があるため、WSL から実行する方が簡単です。 WSL 上の Ubuntu の手順を次に示します。
```
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)
次に、タイルセットを作成して、マップコントロールで使用するか、タイルセットの範囲に含まれる地理座標の標高を取得するか、またはその両方に使用します。タイルセットは Generate XYZ tiles (Directory) ツールを使用して QGIS で作成することができます。
タイルセットの場所を保存します。次のセクションで使用します。

Azure Function と RGB エンコードされた DEM タイルを使用して、Azure Blob Storage から Elevation API を作成する

RGB でエンコードされた DEM タイルは、API を作成するために Azure Functions と共に使用する前に、データベースストレージにアップロードする必要があります。

Azure Blob Storage にタイルをアップロードします。 Azure Storage Explorer は、このために役立つツールです。

タイルを Azure Blob Storage にアップロードするには、数分かかる場合があります。
アップロードが完了したら、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 では、「等高線ベクトルタイルと RGB エンコード DEM タイルを作成する」の手順 3 で作成した EPSG:4326 プロジェクションでマージされたラスタータイルを開きます。
> メニューの [抽出] - [等高線] を選択して、等高線ツールを開きます。

[実行] を選択すると、等高線が作成され、マップにレイヤーとして追加されます。等高線のエッジの一部が少し粗く表示されることがあります。これは、次の手順で対処します。
[処理] メニューから [ツールボックス] を選択して、[処理ツールボックス] を表示します。
次に、[処理ツールボックス] の [ベクタージオメトリ] セクションで [スムーズ] を選択します。

注意

等高線のスムージングは大幅に改善できますが、ファイルサイズが大きくなります。
等高線をデータベースに読み込みます。このガイドでは、localhost で実行される PostgreSQL データベースの無料バージョンを使用します。 Azure Database for PostgreSQL に読み込むこともできます。

次の手順では、PostGIS 拡張機能を持つ PostgreSQL データベースが必要です。
QGIS から PostgreSQL への接続を作成するには、[レイヤー] メニューから > -[PostGIS レイヤーの追加] を選択し、[新規] ボタンを選択します。
次に、QGIS のデータベースマネージャーを使用して QGIS から PostgreSQL にデータを読み込みます。これを行うには、[データベース] メニューから [DB マネージャー] を選択します。
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
        )

コードサンプルの結果を表示するには、ローカルで実行します。

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

フィードバック

このページはお役に立ちましたか?

最終更新: 2023-03-23