Yükseltme verileri ve hizmetleri oluşturma

2024-01-10

Bu kılavuzda, Microsoft Azure Bulut'ta yükseltme hizmeti oluşturmak için SRTM misyonundaki USGS dünya çapındaKI DEM verilerinin 30m doğrulukla nasıl kullanılacağı açıklanmaktadır.

Bu makalede, aşağıdakilerin nasıl yapılacağı açıklanmaktadır:

Kontur çizgi vektör kutucukları ve RGB ile kodlanmış DEM kutucukları oluşturun.
Azure İşlevi ve Azure Blob Depolama RGB ile kodlanmış DEM kutucuklarını kullanarak Yükseltme API'sini oluşturun.
Azure İşlevi ve PostgreSQL kullanarak Kontur çizgi vektör kutucuğu hizmeti oluşturun.

Önkoşullar

Bu kılavuz, aşağıdaki üçüncü taraf yazılım ve verilerin kullanılmasını gerektirir:

USGS Verileri. DEM verileri, USGS EarthExplorer aracılığıyla kutucuk başına 1 ark saniyelik kapsama sahip GeoTiff olarak indirilebilir. Bunun için bir EarthExplorer hesabı gerekir, ancak veriler ücretsiz olarak indirilebilir.
QGIS masaüstü GIS uygulaması, Raster kutucuklarını işlemek ve düzeltmek için kullanılır. QGIS'i indirip kullanabilirsiniz. Bu kılavuzda QGIS sürüm 3.26.2-Buenos Aires kullanılır.
MapBox tarafından geliştirilen rio-rgbify Python paketi, GeoTIFF'yi RGB olarak kodlamak için kullanılır.
PostGIS uzamsal uzantısına sahip PostgreSQL veritabanı.

Kontur çizgi vektör kutucukları ve RGB kodlanmış DEM kutucukları oluşturma

Bu kılavuzda, USGS EarthExplorer tarafından sağlanan Washington eyaletini kapsayan 36 kutucuk kullanılır.

USGS EarthExplorer'dan raster kutucuklarını indirin

Arama ölçütü

Tarama kutucuklarının olmasını istediğiniz bölgeyi seçin. Tanıtım amacıyla, bu kılavuz haritadaki bölgeyi seçmek için "Çokgen" yöntemini kullanır.

USGS EarthExplorer'a gidin.
Arama Ölçütleri sekmesinde Çokgen'i seçin ve ardından harita üzerinde herhangi bir yeri seçerek sınırı oluşturun.

Veri kümeleri

Veri Kümeleri sekmesini seçin.
Dijital Yükseltmeler bölümündeN SRTM 1 Arc-Second Global öğesini seçin.

Sonuçlar

Seçili bölge ve veri kümesinin kutucuklarını görüntülemek için Sonuçlar'ı >> seçin.
İndirilebilir kutucukların listesi sonuçlar sayfasında görünür. Yalnızca istediğiniz kutucukları indirmek için, her kutucuğun sonuç kartında İndirme Seçenekleri düğmesini seçip GeoTIFF 1 Arc-Second seçeneğini belirleyin ve kalan kutucuklar için bu adımı yineleyin.
Alternatif olarak toplu indirme seçeneğini kullanın ve GeoTIFF 1 Arc-second seçeneğini belirleyin.

QGIS'e raster kutucukları ekleme

İhtiyacınız olan raster kutucuklarına sahip olduktan sonra bunları QGIS'te içeri aktarabilirsiniz.

Dosyaları QGIS katmanı sekmesine sürükleyerek veya Katman menüsünde Katman Ekle'yi seçerek QGIS'e raster kutucukları ekleyin.
Raster katmanları QGIS'e yüklendiğinde, kutucukların farklı tonları olabilir. GeoTIFF biçiminde tek bir düzgün raster görüntüsüne neden olan raster katmanlarını birleştirerek bu sorunu düzeltin. Bunu yapmak için, Raster menüsünden Çeşitli'yi ve ardından Birleştir... seçeneğini belirleyin.
birleştirilmiş raster katmanını EPSG:3857 'ye (WGS84 / Pseudo-Mercator) yeniden oluşturarak Raster Katmanını Kaydet'i sağ tıklayarak içerik tablosundaki birleştirilmiş raster katmanına erişin ->Dışarı Aktar ->Farklı Kaydet seçeneği. AZURE HARITALAR Web SDK'sı ile kullanmak için EPSG:3857 gereklidir.
Yalnızca kontur çizgisi vektör kutucukları oluşturmak istiyorsanız aşağıdaki adımları atlayabilir ve Azure İşlevi ve PostgreSQL kullanarak Kontur çizgi vektör kutucuğu oluşturma hizmetine gidebilirsiniz.
Yükseltme API'sini oluşturmak için bir sonraki adım, GEOTIFF'yi RGB ile kodlamaktır. Bu işlem MapBox tarafından geliştirilen rio-rgbify kullanılarak yapılabilir. Bu aracı doğrudan Windows'ta çalıştırmanın bazı zorlukları vardır, bu nedenle WSL'den çalıştırmak daha kolaydır. WSL üzerinde Ubuntu'daki adımlar aşağıdadır:
```
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 kodlu GeoTIFF, bir piksel için R, G ve B değerlerini almanıza ve şu değerlerden yükseltmeyi hesaplamanıza olanak tanır:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
Ardından, harita denetimiyle kullanmak üzere bir kutucuk kümesi oluşturun ve/veya kutucuk kümesinin harita kapsamındaki coğrafi koordinatlar için Yükseltme elde etmek için bu kutucuğu kullanın. Kutucuk kümesi, XYZ kutucukları oluşturma (Dizin) aracı kullanılarak QGIS'te oluşturulabilir.
Kutucuk kümesinin konumunu kaydettiğinizde sonraki Bölümde kullanacaksınız.

Azure İşlevi ve Azure Blob Depolama RGB ile kodlanmış DEM kutucuklarını kullanarak Yükseltme API'sini oluşturma

RGB ile kodlanmış DEM Kutucuklarının bir API oluşturmak için Azure İşlevleri ile kullanılabilmesi için önce bir veritabanı depolama alanına yüklenmesi gerekir.

Kutucukları Azure Blob Depolama yükleyin. Azure Depolama Gezgini bu amaç için yararlı bir araçtır.

Kutucukların Azure Blob Depolama yüklenmesi birkaç dakika sürebilir.
Karşıya yükleme tamamlandıktan sonra, belirli bir coğrafi koordinat için yükseltme döndüren bir API oluşturmak için Azure İşlevi oluşturabilirsiniz.

Bu işlev bir koordinat çifti alır, yakınlaştırma düzeyi 14'te onu kapsayan kutucuğu belirler ve ardından bu kutucuk içinde coğrafi koordinatlarla eşleşen piksel koordinatlarını belirler. Ardından kutucuğu alır, bu pikselin RGB değerlerini alır ve ardından yükseltmeyi belirlemek için aşağıdaki formülü kullanır:

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
        )

Kod örneğinin sonuçlarını görmek için yerel olarak çalıştırın:

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

Azure İşlevi ve PostgreSQL kullanarak Kontur çizgi vektör kutucuğu hizmeti oluşturma

Bu bölümde QGIS'te kontur çizgileri oluşturma ve işleme, bunları PostgreSQL'e yükleme ve ardından vektör kutucuklarını döndürmek için PostgreSQL'i Sorgulamak için bir Azure İşlevi oluşturma adımları açıklanmaktadır.

QGIS'te, Kontur çizgi vektör kutucukları ve RGB ile kodlanmış DEM kutucukları oluşturma işleminin 3. adımında oluşturulan EPSG:4326 projeksiyonunda birleştirilmiş raster kutucuklarını açın.
Kontur aracını açmak için Raster menüsünden Ayıklama -> Kontur'a tıklayın.

Çalıştır seçildiğinde kontur çizgileri oluşturulur ve haritaya katman olarak eklenir. kontur çizgisi kenarlarından bazıları biraz pürüzlü görünebilir. Bu, sonraki adımda ele alınacaktır.
İşleme Araç Kutusu'nu açmak için İşleme menüsünden Araç Kutusu'nu seçin.
Ardından İşleme Araç Kutusu'nun Vektör geometrisi bölümünde Düzgün'e tıklayın.

Not

Kontur çizgisi düzeltme önemli ölçüde iyileştirilebilir, ancak dosya boyutunun artması karşılığında.
Dağılım çizgilerini veritabanına yükleyin. Bu kılavuzda, localhost üzerinde çalışan ücretsiz PostgreSQL veritabanı sürümü gerçekleştirilir. Ayrıca bunları PostgreSQL için Azure Veritabanı yükleyebilirsiniz.

Sonraki adım, PostGIS uzantısına sahip bir PostgreSQL veritabanı gerektirir.
QGIS'ten PostgreSQL'e bağlantı oluşturmak için Katman menüsünden Katman Ekle ->PostGIS Katmanları Ekle'yi ve ardından Yeni düğmesini seçin.
Ardından QGIS'teki Veritabanı Yöneticisi'ni kullanarak QGIS'ten PostgreSQL'e Veri yükleyin. Bunu yapmak için Veritabanı menüsünden DB Yöneticisi'ni seçin.
PostGIS veritabanına bağlanın ve Katman/Dosya İçeri Aktar... öğesini seçerek dağılım çizgilerini veritabanına aktarın.

Artık PostgreSQL'i sorgulamak ve kontur çizgileri için vektör kutucukları döndürmek için bir Azure İşlevi kullanabilirsiniz. Kutucuk sunucusu Azure Haritalar web SDK'sı ile birlikte kullanılarak harita üzerinde dağılım çizgileri görüntüleyen bir web uygulaması oluşturulabilir.

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
        )