إنشاء خدمات وبيانات الارتفاع

يصف هذا الدليل كيفية استخدام بيانات USGS حول العالم DEM من مهمة SRTM الخاصة بهم بدقة 30m لإنشاء خدمة رفع على Microsoft Azure Cloud.

توضح هذه المقالة كيفية:

إنشاء تجانبات خط خط إحاطة ولوحات DEM المرمزة ب RGB.
إنشاء واجهة برمجة تطبيقات رفع باستخدام Azure Function ولوحات DEM المشفرة RGB من Azure Blob Storage.
إنشاء خدمة تجانب خط خط Contour باستخدام Azure Function وPostgreSQL.

المتطلبات الأساسية

يتطلب هذا الدليل استخدام برامج وبيانات الجهات الخارجية التالية:

بيانات USGS. يمكن تنزيل بيانات DEM ك GeoTiff مع تغطية ثانية واحدة للقوس لكل تجانب من خلال USGS EarthExplorer. يتطلب هذا حساب EarthExplorer، ولكن يمكن تنزيل البيانات مجانا.
يستخدم تطبيق QGIS لسطح المكتب GIS لمعالجة تجانبات Raster وتنعيمها. QGIS مجاني للتنزيل والاستخدام. يستخدم هذا الدليل الإصدار 3.26.2 من QGIS- بوينس آيرس.
يتم استخدام حزمة ريو-rgbify Python، التي تم تطويرها بواسطة MapBox، لترميز GeoTIFF ك RGB.
قاعدة بيانات PostgreSQL مع ملحق PostGIS المكاني.

إنشاء تجانبات خط خط إحاطة ولوحات DEM المرمزة ب RGB

يستخدم هذا الدليل الإطارات المتجانبة ال 36 التي تغطي ولاية واشنطن، والمتوفرة من USGS EarthExplorer.

تنزيل اللوحات النقطية من USGS EarthExplorer

معايير البحث

حدد المنطقة التي تريد إنشاء لوحات نقطية لها. لأغراض العرض التوضيحي، يستخدم هذا الدليل طريقة "المضلع" لتحديد المنطقة على الخريطة.

انتقل إلى USGS EarthExplorer.
في علامة التبويب معايير البحث، حدد مضلع ثم حدد أي مكان على الخريطة لإنشاء الحدود.

مجموعات البيانات

حدد علامة التبويب مجموعات البيانات.
حدد SRTM 1 Arc-Second Global من قسم Digital Elevations .

النتائج

حدد النتائج >> لعرض الإطارات المتجانبة للمنطقة المحددة ومجموعة البيانات.
تظهر قائمة اللوحات القابلة للتنزيل في صفحة النتائج. لتنزيل اللوحات التي تريدها فقط، حدد الزر خيارات التنزيل على بطاقة النتائج لكل لوحة، وحدد الخيار GeoTIFF 1 Arc-Second وكرر هذه الخطوة لللوحات المتبقية.
بدلا من ذلك، استخدم خيار التنزيل المجمع وحدد GeoTIFF 1 Arc-second.

إضافة تجانبات نقطية إلى QGIS

بمجرد أن يكون لديك الإطارات المتجانبة النقطية التي تحتاجها، يمكنك استيرادها في QGIS.

أضف تجانبات نقطية إلى QGIS عن طريق سحب الملفات إلى علامة تبويب طبقة QGIS أو تحديد إضافة طبقة في قائمة الطبقة .
عند تحميل طبقات النقط في QGIS، يمكن أن تكون هناك ظلال مختلفة من الإطارات المتجانبة. قم بإصلاح ذلك عن طريق دمج طبقات النقط، مما يؤدي إلى صورة نقطية واحدة متجانسة بتنسيق GeoTIFF. للقيام بذلك، حدد Miscellaneous من قائمة Raster ، ثم Merge...
أعد مشروع طبقة النقط المدمجة إلى EPSG:3857 (WGS84 / Pseudo-Mercator) باستخدام حفظ طبقة نقطية بالنقر بزر الماوس الأيمن للوصول إلى طبقة النقط المدمجة في جدول المحتوى ->خيار التصدير ->حفظ باسم . مطلوب EPSG:3857 لاستخدامه مع خرائط Azure Web SDK.
إذا كنت تريد فقط إنشاء إطارات متجانبة لناقل خط إحاطة، يمكنك تخطي الخطوات التالية والانتقال إلى إنشاء خدمة تجانب خط Contour باستخدام Azure Function وPostgreSQL.
لإنشاء واجهة برمجة تطبيقات رفع، الخطوة التالية هي ترميز RGB-Encode في GeoTIFF. يمكن القيام بذلك باستخدام rio-rgbify، الذي تم تطويره بواسطة MapBox. هناك بعض التحديات التي تواجه تشغيل هذه الأداة مباشرة في Windows، لذلك من الأسهل تشغيلها من WSL. فيما يلي الخطوات الواردة في Ubuntu على 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 المشفرة RGB باسترداد قيم R وG وB لبكسل وحساب الارتفاع من هذه القيم:

elevation (m) = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
بعد ذلك، أنشئ مجموعة تجانب لاستخدامها مع عنصر تحكم الخريطة و/أو استخدمها للحصول على رفع لأي إحداثيات جغرافية ضمن نطاق الخريطة لمجموعة الإطارات المتجانبة. يمكن إنشاء مجموعة التجانب في QGIS باستخدام أداة إنشاء تجانبات XYZ (الدليل ).
احفظ موقع مجموعة الإطارات المتجانبة، وستستخدمها في المقطع التالي.

إنشاء واجهة برمجة تطبيقات رفع باستخدام Azure Function ولوحات DEM المشفرة RGB من Azure Blob Storage

يجب تحميل تجانبات DEM المشفرة RGB إلى تخزين قاعدة بيانات قبل أن يمكن استخدامها مع Azure Functions لإنشاء واجهة برمجة تطبيقات.

تحميل الإطارات المتجانبة إلى Azure Blob Storage. Azure Storage Explorer هو أداة مفيدة لهذا الغرض.

قد يستغرق تحميل الإطارات المتجانبة إلى Azure Blob Storage عدة دقائق حتى يكتمل.
بمجرد اكتمال التحميل، يمكنك إنشاء Azure Function لإنشاء واجهة برمجة تطبيقات تقوم بإرجاع الارتفاع لإحداثيات جغرافية معينة.

تتلقى هذه الدالة زوجا من الإحداثيات، وتحدد اللوحة التي تغطيها في مستوى التكبير/التصغير 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`

إنشاء خدمة تجانب خط خط Contour باستخدام Azure Function وPostgreSQL

يصف هذا القسم خطوات إنشاء خطوط إحاطة ومعالجتها في QGIS، وتحميلها إلى PostgreSQL ثم إنشاء وظيفة Azure للاستعلام عن PostgreSQL لإرجاع تجانبات المتجهات.

في QGIS، افتح الإطارات المتجانبة النقطية المدمجة في إسقاط EPSG:4326 الذي تم إنشاؤه في الخطوة 3 من إنشاء تجانبات خط خط Contour ولوحات DEM المشفرة RGB.
حدد Extraction -> Contour من قائمة Raster لفتح أداة Contour.

يؤدي تحديد Run إلى إنشاء خطوط إحاطة وإضافتها كطبقة إلى الخريطة. قد تظهر بعض حواف خط إحاطة تقريبية بعض الشيء. ستتم معالجة هذا في الخطوة التالية.
حدد Toolbox من قائمة Processing لإظهار مربع أدوات المعالجة.
ثم حدد Smooth في قسم Vector geometry في مربع أدوات المعالجة.

إشعار

يمكن تحسين تجانس خط إحاطة بشكل كبير ولكن على حساب زيادة حجم الملف.
تحميل خطوط إحاطة إلى قاعدة البيانات. يستخدم هذا الدليل الإصدار المجاني من قاعدة بيانات PostgreSQL التي تعمل على المضيف المحلي. يمكنك أيضا تحميلها إلى قاعدة بيانات Azure ل PostgreSQL.

تتطلب الخطوة التالية قاعدة بيانات PostgreSQL مع ملحق PostGIS .
لإنشاء اتصال من QGIS إلى PostgreSQL، حدد إضافة طبقة ->إضافة طبقات PostGIS من قائمة الطبقة ، ثم حدد الزر جديد .
بعد ذلك، قم بتحميل البيانات من QGIS إلى PostgreSQL باستخدام Database Manager في QGIS. للقيام بذلك، حدد DB Manager من قائمة قاعدة البيانات .
اتصل بقاعدة بيانات PostGIS وحدد Import Layer/File... لاستيراد خطوط إحاطة إلى قاعدة البيانات.

يمكنك الآن استخدام دالة Azure للاستعلام عن PostgreSQL وإرجاع الإطارات المتجانبة المتجهة لخطوط إحاطة. يمكن استخدام خادم التجانب مع Azure Maps Web SDK لإنشاء تطبيق ويب يعرض خطوط إحاطة على الخريطة.

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
        )