Bagikan melalui

Membuat data & layanan elevasi

Panduan ini menjelaskan cara menggunakan data DEM usgs di seluruh dunia dari misi SRTM mereka dengan akurasi 30m untuk membangun layanan Elevation di Microsoft Azure Cloud.

Artikel ini menjelaskan cara:

  • Buat petak peta vektor garis Contour dan petak peta DEM yang dikodekan RGB.
  • Buat API Elevasi menggunakan Azure Function dan petak peta DEM yang dikodekan RGB dari Azure Blob Storage.
  • Buat layanan petak peta vektor baris Contour menggunakan Azure Function dan PostgreSQL.

Prasyarat

Panduan ini memerlukan penggunaan perangkat lunak dan data pihak ketiga berikut:

  • Data USGS. Data DEM dapat diunduh sebagai GeoTiff dengan cakupan kedua 1 busur per petak peta melalui USGS EarthExplorer. Ini memerlukan akun EarthExplorer, tetapi data dapat diunduh secara gratis.
  • Aplikasi GIS desktop QGIS digunakan untuk memproses dan menghaluskan petak peta Raster. QGIS gratis untuk diunduh dan digunakan. Panduan ini menggunakan QGIS versi 3.26.2-Buenos Aires.
  • Paket rio-rgbify Python, yang dikembangkan oleh MapBox, digunakan untuk mengodekan GeoTIFF sebagai RGB.
  • Database PostgreSQL dengan ekstensi spasial PostGIS .

Membuat petak peta vektor garis Contour dan petak peta DEM yang dikodekan RGB

Panduan ini menggunakan 36 petak peta yang mencakup negara bagian Washington, tersedia dari USGS EarthExplorer.

Unduh petak peta raster dari USGS EarthExplorer

Kriteria pencarian

Pilih wilayah yang Anda inginkan untuk petak peta raster. Untuk tujuan demonstrasi, panduan ini menggunakan metode "Poligon" untuk memilih wilayah di peta.

  1. Navigasi ke USGS EarthExplorer.

  2. Di tab Kriteria Pencarian, pilih Poligon lalu pilih di mana saja di peta untuk membuat batas.

    Cuplikan layar memperlihatkan tab kriteria pencarian di situs web penjelajah bumi USGS.

Himpunan data

  1. Pilih tab Himpunan Data.

  2. Pilih SRTM 1 Arc-Second Global dari bagian Elevasi Digital.

    Cuplikan layar memperlihatkan tab himpunan data di situs web penjelajah bumi USGS.

Hasil

  1. Pilih Hasil >> untuk menampilkan petak peta untuk wilayah dan himpunan data yang dipilih.

  2. Daftar petak peta yang dapat diunduh muncul di halaman hasil. Untuk mengunduh hanya petak peta yang Anda inginkan, pilih tombol Opsi Unduh pada kartu hasil untuk setiap petak peta, memilih opsi GeoTIFF 1 Arc-Second dan ulangi langkah ini untuk petak peta yang tersisa.

    Cuplikan layar memperlihatkan tab hasil di situs web usgs earth explorer.

  3. Atau, gunakan opsi unduh massal dan pilih GeoTIFF 1 Arc-second.

Menambahkan petak peta raster ke QGIS

Setelah Anda memiliki petak peta raster yang Anda butuhkan, Anda dapat mengimpornya di QGIS.

  1. Tambahkan petak peta raster ke QGIS dengan menyeret file ke tab lapisan QGIS atau memilih Tambahkan Lapisan di menu Lapisan .

    Cuplikan layar memperlihatkan petak peta raster di QGIS.

  2. Ketika lapisan raster dimuat ke dalam QGIS, mungkin ada nuansa petak peta yang berbeda. Perbaiki ini dengan menggabungkan lapisan raster, yang menghasilkan satu gambar raster halus dalam format GeoTIFF. Untuk melakukan ini, pilih Lain-lain dari menu Raster , lalu Gabungkan...

    Cuplikan layar memperlihatkan menu gabungkan raster di QGIS.

  3. Proyek ulang lapisan raster gabungan ke EPSG:3857 (WGS84 / Pseudo-Mercator) menggunakan Simpan Lapisan Raster sebagai klik kanan untuk mengakses lapisan raster gabungan dalam daftar isi ->Ekspor ->Simpan Sebagai opsi. EPSG:3857 diperlukan untuk menggunakannya dengan Azure Maps Web SDK.

    Cuplikan layar yang menunjukkan bagaimana menu gabungkan lapisan raster di QGIS.

  4. Jika Anda hanya ingin membuat petak peta vektor garis kontur, Anda dapat melewati langkah-langkah berikut dan masuk ke Membuat layanan petak peta vektor baris Contour menggunakan Azure Function dan PostgreSQL.

  5. Untuk membuat API Elevasi, langkah selanjutnya adalah RGB-Encode GeoTIFF. Ini dapat dilakukan menggunakan rio-rgbify, yang dikembangkan oleh MapBox. Ada beberapa tantangan menjalankan alat ini langsung di Windows, sehingga lebih mudah dijalankan dari WSL. Di bawah ini adalah langkah-langkah dalam Ubuntu di 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/

    Cuplikan layar memperlihatkan GeoTIFF yang dikodekan RGB di QGIS.

    GeoTIFF yang dikodekan RGB memungkinkan Anda mengambil nilai R, G, dan B untuk piksel dan menghitung elevasi dari nilai-nilai ini:

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

  6. Selanjutnya, buat set petak peta yang akan digunakan dengan kontrol peta dan/atau gunakan untuk mendapatkan Elevasi untuk koordinat geografis apa pun dalam tingkat peta set petak peta. Set petak peta dapat dibuat di QGIS menggunakan alat Hasilkan petak peta (Direktori) XYZ.

    Cuplikan layar yang menunjukkan alat Hasilkan petak peta (Direktori) XYZ di QGIS.

  7. Simpan lokasi kumpulan petak peta, Anda akan menggunakannya di Bagian berikutnya.

Membuat API Elevasi menggunakan Azure Function dan petak peta DEM yang dikodekan RGB dari Azure Blob Storage

Ubin DEM yang dikodekan RGB perlu diunggah ke penyimpanan database sebelum dapat digunakan dengan Azure Functions untuk membuat API.

  1. Unggah petak peta ke Azure Blob Storage. Azure Storage Explorer adalah alat yang berguna untuk tujuan ini.

    Cuplikan layar memperlihatkan Microsoft Azure Storage Explorer.

    Mengunggah petak peta ke Azure Blob Storage dapat memakan waktu beberapa menit untuk diselesaikan.

  2. Setelah unggahan selesai, Anda dapat membuat Azure Function untuk membangun API yang mengembalikan elevasi untuk koordinat geografis tertentu.

    Fungsi ini menerima pasangan koordinat, menentukan petak peta yang mencakupnya pada tingkat zoom 14, lalu menentukan koordinat piksel dalam petak peta yang cocok dengan koordinat geografis. Kemudian mengambil petak peta, mendapatkan nilai RGB untuk piksel tersebut, lalu menggunakan rumus berikut untuk menentukan elevasi:

    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
        )

Untuk melihat hasil sampel kode, jalankan secara lokal:

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

Membuat layanan petak peta vektor baris Contour menggunakan Azure Function dan PostgreSQL

Bagian ini menjelaskan langkah-langkah untuk membuat dan memproses garis kontur di QGIS, mengunggahnya ke PostgreSQL lalu membuat Fungsi Azure ke Query PostgreSQL untuk mengembalikan petak peta vektor.

  1. Di QGIS, buka petak peta raster gabungan dalam proyeksi EPSG:4326 yang dibuat di langkah 3 dari Buat petak peta vektor baris Contour dan petak peta DEM yang dikodekan RGB.

  2. Pilih Ekstraksi -> Kontur dari menu Raster untuk membuka alat Contour.

    Cuplikan layar memperlihatkan dialog kontur di QGIS.

    Memilih Jalankan membuat garis kontur dan menambahkannya sebagai lapisan ke peta. beberapa tepi garis kontur mungkin tampak sedikit kasar. Ini akan ditangani pada langkah berikutnya.

    Cuplikan layar memperlihatkan peta dengan kontur di QGIS.

  3. Pilih Kotak Alat dari menu Pemrosesan untuk memunculkan Kotak Alat Pemrosesan.

  4. Kemudian pilih Smooth di bagian geometri Vektor dari Kotak Alat Pemrosesan.

    Cuplikan layar memperlihatkan dialog lancar di QGIS.

    Catatan

    Penghalusan garis kontur dapat ditingkatkan secara substansial tetapi dengan biaya peningkatan ukuran file.

  5. Muat baris kontur ke database. Panduan ini menggunakan versi gratis database PostgreSQL yang berjalan di localhost. Anda juga dapat memuatnya ke Azure Database for PostgreSQL.

    Langkah berikutnya memerlukan database PostgreSQL dengan ekstensi PostGIS .

  6. Untuk membuat koneksi dari QGIS ke PostgreSQL, pilih Tambahkan Lapisan ->Tambahkan Lapisan PostGIS dari menu Lapisan, lalu pilih tombol Baru.

    Cuplikan layar memperlihatkan dialog buat koneksi PostGIG baru di QGIS.

  7. Selanjutnya, muat Data dari QGIS ke PostgreSQL menggunakan Database Manager di QGIS. Untuk melakukan ini, pilih Manajer DB dari menu Database .

    Cuplikan layar memperlihatkan Manajer DB di QGIS.

  8. Sambungkan ke database PostGIS dan pilih Impor Lapisan/File... untuk Mengimpor baris kontur ke database.

    Cuplikan layar memperlihatkan dialog impor vektor di QGIS.

  9. Sekarang Anda dapat menggunakan Fungsi Azure untuk Mengkueri PostgreSQL dan mengembalikan petak peta vektor untuk garis kontur. Server petak peta dapat digunakan dengan SDK web Azure Maps untuk membuat aplikasi web yang menampilkan garis kontur di peta.

    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
        )

Untuk melihat hasil sampel kode, jalankan secara lokal:

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