Catatan
Akses ke halaman ini memerlukan otorisasi. Anda dapat mencoba masuk atau mengubah direktori.
Akses ke halaman ini memerlukan otorisasi. Anda dapat mencoba mengubah direktori.
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.
Navigasi ke USGS EarthExplorer.
Di tab Kriteria Pencarian, pilih Poligon lalu pilih di mana saja di peta untuk membuat batas.
Himpunan data
Hasil
Pilih Hasil >> untuk menampilkan petak peta untuk wilayah dan himpunan data yang dipilih.
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.
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.
Tambahkan petak peta raster ke QGIS dengan menyeret file ke tab lapisan QGIS atau memilih Tambahkan Lapisan di menu Lapisan .
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...
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.
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.
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/
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)
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.
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.
Unggah petak peta ke Azure Blob Storage. Azure Storage Explorer adalah alat yang berguna untuk tujuan ini.
Mengunggah petak peta ke Azure Blob Storage dapat memakan waktu beberapa menit untuk diselesaikan.
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.
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.
Pilih Ekstraksi -> Kontur dari menu Raster untuk membuka alat Contour.
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.
Pilih Kotak Alat dari menu Pemrosesan untuk memunculkan Kotak Alat Pemrosesan.
Kemudian pilih Smooth di bagian geometri Vektor dari Kotak Alat Pemrosesan.
Catatan
Penghalusan garis kontur dapat ditingkatkan secara substansial tetapi dengan biaya peningkatan ukuran file.
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 .
Untuk membuat koneksi dari QGIS ke PostgreSQL, pilih Tambahkan Lapisan ->Tambahkan Lapisan PostGIS dari menu Lapisan, lalu pilih tombol Baru.
Selanjutnya, muat Data dari QGIS ke PostgreSQL menggunakan Database Manager di QGIS. Untuk melakukan ini, pilih Manajer DB dari menu Database .
Sambungkan ke database PostGIS dan pilih Impor Lapisan/File... untuk Mengimpor baris kontur ke database.
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}