영어로 읽기

다음을 통해 공유


상승 데이터 및 서비스 만들기

이 가이드에서는 30m 정확도의 SRTM 임무에 맞춰 USGS 전 세계 DEM 데이터를 사용하여 Microsoft Azure Cloud에서 고도 서비스를 빌드하는 방법을 설명합니다.

이 문서에서는 다음과 같이 방법을 설명합니다.

  • 윤곽선 벡터 타일 및 RGB로 인코딩된 DEM 타일을 만듭니다.
  • Azure Function을 사용하여 고도 API를 만들고 Azure Blob Storage에서 RGB로 인코딩된 DEM 타일을 만듭니다.
  • Azure Function 및 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로 인코딩된 DEM 타일 만들기

이 가이드에서는 USGS EarthExplorer에서 사용할 수 있는 워싱턴 주를 포괄하는 36개 타일을 사용합니다.

USGS EarthExplorer에서 래스터 타일 다운로드

검색 조건

래스터 타일을 사용할 지역을 선택합니다. 데모를 위해 이 가이드에서는 "Polygon" 메서드를 사용하여 지도에서 지역을 선택합니다.

  1. USGS EarthExplorer로 이동합니다.

  2. 검색 기준 탭에서 다각형을 선택한 다음 맵의 아무 곳이나 선택하여 경계를 만듭니다.

    USGS 지구 탐색기 웹 사이트의 검색 기준 탭을 보여 주는 스크린샷.

데이터 집합

  1. 데이터 세트 탭을 선택합니다.

  2. 디지털 고도 섹션에서 SRTM 1 Arc-Second Global을 선택합니다.

    USGS 지구 탐색기 웹 사이트의 데이터 세트 탭을 보여 주는 스크린샷.

결과

  1. 결과>>를 선택하여 선택한 지역 및 데이터 세트의 타일을 봅니다.

  2. 다운로드 가능한 타일 목록이 결과 페이지에 표시됩니다. 원하는 타일만 다운로드하려면 각 타일에 대한 결과 카드에서 다운로드 옵션 단추를 선택하고 GeoTIFF 1 Arc-Second 옵션을 선택한 후 나머지 타일에 대해 이 단계를 반복합니다.

    USGS 지구 탐색기 웹 사이트의 결과 탭을 보여 주는 스크린샷.

  3. 또는 대량 다운로드 옵션을 사용하고 GeoTIFF 1 Arc-second를 선택합니다.

QGIS에 래스터 타일 추가

필요한 래스터 타일이 있으면 QGIS에서 가져올 수 있습니다.

  1. 파일을 QGIS 계층 탭으로 끌어오거나 계층 메뉴에서 계층 추가를 선택하여 래스터 타일을 QGIS에 추가합니다.

    QGIS에서 래스터 타일을 보여 주는 스크린샷.

  2. 래스터 계층이 QGIS에 로드되면 타일의 음영이 다를 수 있습니다. 래스터 계층을 병합하여 이 문제를 해결하면 하나의 매끄러운 래스터 이미지가 GeoTIFF 형식으로 생성됩니다. 이렇게 하려면 래스터 메뉴에서 기타를 선택한 다음, 병합...을 선택합니다.

    QGIS의 래스터 병합 메뉴를 보여 주는 스크린샷.

  3. 목차 ->내보내기 ->다른 이름으로 저장 옵션의 병합된 래스터 계층을 마우스 오른쪽 단추로 클릭하여 액세스되는 래스터 계층으로 저장을 사용하여 병합된 래스터 계층을 EPSG:3857(WGS84/Pseudo-Mercator)로 재투영합니다. Azure Maps Web SDK에서 사용하려면 EPSG:3857이 필요합니다.

    QGIS에서 래스터 계층 병합 메뉴를 보여 주는 스크린샷.

  4. 윤곽선 벡터 타일만 만들려면 다음 단계를 건너뛰고 Azure Function 및 PostgreSQL을 사용하여 윤곽선 벡터 타일 서비스 만들기로 이동합니다.

  5. 고도 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/
    

    QGIS에서 RGB로 인코딩된 GeoTIFF를 보여 주는 스크린샷.

    RGB로 인코딩된 GeoTIFF를 사용하면 픽셀의 R, G 및 B 값을 검색하고 이러한 값에서 고도를 계산할 수 있습니다.

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

  6. 다음으로, 지도 컨트롤과 함께 사용할 타일 집합을 만들거나, 타일 집합의 지도 범위 내에서 지리적 좌표의 고도를 가져오는 데 사용합니다. 타일 집합은 XYZ 타일 생성(디렉터리) 도구를 사용하여 QGIS에서 만들 수 있습니다.

    QGIS의 XYZ 타일 생성(디렉터리) 도구를 보여 주는 스크린샷.

  7. 타일 집합의 위치를 저장하면 다음 섹션에서 사용합니다.

Azure Function을 사용하여 고도 API를 만들고 Azure Blob Storage에서 RGB로 인코딩된 DEM 타일 만들기

RGB로 인코딩된 DEM 타일을 데이터베이스 스토리지에 업로드해야 Azure Functions에서 API를 만드는 데 사용할 수 있습니다.

  1. 타일을 Azure Blob Storage에 업로드합니다. Azure Storage Explorer는 이 용도에 유용한 도구입니다.

    Microsoft Azure Storage Explorer를 보여 주는 스크린샷.

    Azure Blob Storage로의 타일 업로드를 완료하는 데 몇 분 정도 걸릴 수 있습니다.

  2. 업로드가 완료되면 Azure Function을 만들어 지정된 지리적 좌표에 대한 고도를 반환하는 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 Function 및 PostgreSQL을 사용하여 윤곽선 벡터 타일 서비스 만들기

이 섹션에서는 QGIS에서 윤곽선을 만들어 처리하고, PostgreSQL에 업로드한 다음, PostgreSQL을 쿼리하여 벡터 타일을 반환하는 Azure 함수를 만드는 단계를 설명합니다.

  1. QGIS에서 윤곽선 벡터 타일 및 RGB로 인코딩된 DEM 타일 만들기의 3단계에서 만든 EPSG:4326 도법에서 병합된 래스터 타일을 엽니다.

  2. 래스터 메뉴에서 추출 -> 윤곽선을 선택하여 윤곽선 도구를 엽니다.

    QGIS의 등고선 대화 상자를 보여 주는 스크린샷.

    실행을 선택하면 윤곽선이 만들어지고 지도에 계층으로 추가됩니다. 일부 윤곽선 가장자리는 약간 거칠게 나타날 수 있습니다. 이 문제는 다음 단계에서 해결됩니다.

    QGIS에서 등고선이 있는 맵을 보여 주는 스크린샷.

  3. 처리 메뉴에서 도구 상자를 선택하여 처리 도구 상자를 표시합니다.

  4. 그런 다음, 처리 도구 상자벡터 기하 도형 섹션에서 다듬기를 선택합니다.

    QGIS의 원활한 대화 상자를 보여 주는 스크린샷.

    참고

    윤곽선 다듬기를 크게 개선할 수 있지만 파일 크기가 증가합니다.

  5. 데이터베이스에 윤곽선을 로드합니다. 이 가이드에서는 localhost에서 실행되는 PostgreSQL 데이터베이스의 무료 버전을 사용합니다. Azure Database for PostgreSQL에 로드할 수도 있습니다.

    다음 단계에서는 PostGIS 확장이 있는 PostgreSQL 데이터베이스가 필요합니다.

  6. QGIS에서 PostgreSQL로의 연결을 만들려면 계층 메뉴에서 계층 추가 ->PostGIS 계층 추가를 선택한 다음, 새로 만들기 단추를 선택합니다.

    QGIS에서 새로운 PostGIG 연결 만들기 대화 상자를 보여 주는 스크린샷.

  7. 다음으로, QGIS의 데이터베이스 관리자를 사용하여 QGIS에서 PostgreSQL로 데이터를 로드합니다. 이렇게 하려면 데이터베이스 메뉴에서 DB 관리자를 선택합니다.

    QGIS의 DB 관리자를 보여 주는 스크린샷.

  8. PostGIS 데이터베이스에 연결하고 계층/파일 가져오기...를 선택하여 데이터베이스에 윤곽선을 가져옵니다.

    QGIS에서 벡터 가져오기 대화 상자를 보여 주는 스크린샷.

  9. 이제 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
            )
    

코드 샘플의 결과를 보려면 로컬로 실행합니다.

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