Pythonで国土地理院DEMデータを利用するときのお役立ちCode & Function（その１）

日付:

Thema: Code and Functions of processing DEM data of GSI
(Geospatial Information Authority) in Japan){Python)

国土地理院(GSI)のDigital Elevtation Modelは地理院Tileという概念を用いているため、この地理院Tile座標やTile内部の個々のDEM値の位置を世界座標系で表現するPixel座標を、他のGIS(Geographic Information System)で利用されることが多い緯度経度(Longitude, Latitude)に変換しないといけない場面に遭遇します。

The case like thatのために使用するFuctionなどをIntroduceしたいと思います。

１．Tile座標を緯度経度に変換するFunction

最初に、あるTile座標とZoomがgiveされたときに、そのTileのNorth-WestEdgeのLon,Lat(Longitude,Latitude)をReturnするFunction 'tile_to_lolat'です。 


# Function of exchange to lon, lat of the tile's north-west edge from the tile zoom and coodinetes / TileのZoomとx,y座標からタイル北西端の緯度経度を計算する関数
# z:zoom, x:タイルのx座標, y:タイルのy座標
def tile_to_lonlat(z, x, y):
    lon = (x / 2.0**z) * 360 - 180 # 経度（東経）
    mapy = (y / 2.0**z) * 2 * pi - pi
    lat = 2 * atan(e ** (- mapy)) * 180 / pi - 90 # 緯度（北緯）
    return lon,lat

Google Geminiさんに聞くと、Latitudeの計算方法は以下の２種類あるようですが、どちらの式を使っても正しい緯度を計算できるそうです。

①lat = math.atan(math.sinh(math.pi * (1 - 2 * y / 2**z))) * 180 / math.pi

②mapy = (y / 2.0**z) * 2 * pi - pi, lat = 2 * atan(e ** (- mapy)) * 180 / pi - 90 

Hyperbaric関数を解くのがめんどくさいので、ここはGoogleGeminiさんを信じることにします。

２．Tile内の標高値Meshの平均的なSizeをCalculateするFunction

Tile内には東西256×南北256個のElevationValueが格納されており、その１つ１つのElevationValueMeshの東西・南北の間隔(距離)とOneMeshの面積をCaluculateしてReturnするFunctionでございます。


# Function that calculatete the avarege mesh size and area(m2) from zoom, tile coordinate north-west edge and south-east edge
# / タイル内のメッシュ平均サイズ、平均メッシュ面積(m2単位)を計算する関数
def mesh_size_inner_tile(z, wn_x, wn_y, es_x, es_y):
    lon1, lat1 = tile_to_lonlat(z, wn_x, wn_y)
    lon2, lat2 = tile_to_lonlat(z, es_x, es_y)
    lon_size = (lon2-lon1) / (256 * (es_x - wn_x)) # 角度刻み
    lat_size = (lat2-lat1) / (256 * (es_y - wn_y)) # 角度刻み
    lon_pich = lonlat_to_dist(lon1, lat1, lon2, lat1) / (256 * (es_x - wn_x)) # 距離刻み
    lat_pich = lonlat_to_dist(lon1, lat1, lon1, lat2) / (256 * (wn_y - es_y)) # 距離刻み
    area = abs(lon_pich * lat_pich)
    return lon_size, lat_size, lon_pich, lat_pich, area

引数zはZoomLevel。wn_x, wn_y, es_x, es_yは今回対象としている検討領域のWest-North EdgeおよびはEast-SouthののTile座標でございます。　

lon_size, lat_sizeは緯度・経度の角度単位で、lon_pich, lat_pichは距離(m)単位の表示、areaはOneMeshの面積(m2)でございます。

３．緯度・経度とZoomLevelからPixel座標およびTile内部(局所)座標(Inner(Local) Tile Coordinate)をCalculateするFunction

緯度・経度(Lon, Lat)からあるZoomLevelにおけるPixel座標を、また、複数のTileを結合した結合TileModelの北西端を(0,0)とした内部(局所)Pixel座標(わたしは勝手にCol(Column:列), Row(行)と名付けました)

当初、１つのTileに対してTile北西端緯度・経度－Tile南東端緯度・経度を求めて、それを緯度差分÷256、経度差分÷256みたいな単純計算で行けるのかなぁと思って、Col, Rowといった内部(局所)座標を試行しましたが、そんな単純な方法では緯度経度から変換した位置が全く異なるということがわかり、緯度・経度とZoomから地球全体のPixel座標を求める変換式を探しました。

ChatGPTが公開される前だったのでWeb上を一生懸命探しました。ヒットしたのは「Trail Note」さまの以下のURLでございます。

https://www.trail-note.net/tech/coordinate/

緯度経度と地理院TileのPixel座標の関係を説明していただいている貴重な日本語のWebSiteでございます。本当にありがとうございます。このWebSiteを参考にして、Functionを以下のように定義いたします。

引数はz=ZoomLevel. lonlat=(Longitude, Latitude), edge = (wn_edge_lon, wn_edge_lat)[検討領域北西端の緯度、経度]でございます。


# Function that exchange to column and row number in the conbined elevation model from zoom and north-west edge lon,kat
# / 任意の緯度・経度について、結合した標高値モデルのzoomと北西端緯度・経度から結合標高値モデルにおける行列番号に変換する関数
def func_lonlat_to_colrow(z, lonlat, nw_edge_lonlat):

    parameter_l = 85.05112878
    lon, lat = lonlat
    wn_edge_lon, wn_edge_lat = nw_edge_lonlat

# Calculate inner(local) tile column, row and the Pixel Coordinate(Global Coordinate) from zoom and lon, lat
# / zoomと緯度経度から全地球系のPixcel座標を計算
    global_column = int(2 ** (z + 7) * (lon / 180 + 1))
    global_row = int(2 ** (z + 7) / pi * (-1 * np.arctanh(np.sin(pi * lat / 180))+np.arctanh(np.sin(pi / 180 * parameter_l))))

# Calculate the Pixel coordinate of west and north edge of the combined elevation model
    w_edge_column =  int(2 ** (z + 7) * (wn_edge_lon / 180 + 1))
    n_edge_row = int(2 ** (z + 7) / pi * (-1 * np.arctanh(np.sin(pi * wn_edge_lat / 180)) + np.arctanh(np.sin(pi / 180 * parameter_l))))

# Calculate the row and column number in the local coodinate system(at conbined elevation model)
    column = global_column - w_edge_column
    row = global_row - n_edge_row

    return column ,row

4．Tile内部(局所)座標(Inner(Local) Tile Coordinate)とZoomLevelから緯度経度(LonLat)をCalculateするFunction

前述の３．のReverseのTransformでございます。

引数はzoom=ZoomLevel. colrow=(column, Row)[対象領域北西端座標を[0, 0]としたときの任意の標高点の行列番号], nw_edge_lonlat = (wn_edge_lon, wn_edge_lat)[検討領域北西端の緯度、経度]でございます。


# Calculate lon, lat from the inner(local) tile column, row based on zoom, north-west edge lon, lat of the consideration area
# 任意のLocalColumnRowNumber'colrow'からLonLat(経度緯度)をCalculateするFunction'colrow_to_lonlat()'
def func_colrow_to_lonlat_v02(zoom, colrow, nw_edge_lonlat):

    parameter_l = 85.05112878
    col, row = colrow
    nw_edge_lon, nw_edge_lat = nw_edge_lonlat

# LocalTile北西端のPixelCoordをCalculate
    nw_edge_col = int(2 ** (zoom + 7) * (nw_edge_lon / 180 + 1))
    nw_edge_row = int(2 ** (zoom + 7) / pi * (-1 * np.arctanh(np.sin(pi * nw_edge_lat / 180)) + np.arctanh(np.sin(pi / 180 * parameter_l))))

# LocalTileのColumnRowNumberをPixelCoordにExchange
    pixcel_col = col + nw_edge_col
    pixcel_row = row + nw_edge_row

    lon = 180 * (pixcel_col / 2 ** (zoom + 7) - 1)
    lat = 180/pi * (np.arcsin(np.tanh(-1 * pi / (2**(zoom+7)) * pixcel_row + np.arctanh(np.sin(pi/180 *parameter_l)))))

    return lon, lat


５．LonLan(緯度・経度)が既知の2点間のDistanceをCalcするFunction’lonlat_to_dist'

　緯度・経度が既知の２点間の距離をCalcする関数です。２点の緯度経度(lon1, lat1), (lon2, lat2)を引数として与えて、２点間の距離distanceをReturnします。

import pyproj

# Calculate the distance of two ponts lon, lat
# 緯度経度が与えられた2点間距離を計算する関数'lonlat_to_dist()'
def lonlat_to_dist(lon1, lat1, lon2, lat2):
    grs80 = pyproj.Geod(ellps='GRS80')
    az, bkw_az, dist = grs80.inv(lon1, lat1, lon2, lat2)
    return dist

　かなりPastにCreateした関数なので中身を忘却してしまいました。中身についてGemini先生に聞いてみると以下の通りの回答をいただきました。

  • pyproj.Geod を使用して、GRS80 楕円体に基づく地理座標計算オブジェクトを作成します。

  • grs80.inv() メソッドを使用して、2点間の距離を計算します。このメソッドは、方位角、逆方位角、距離を返しますが、ここでは距離のみを使用します。

    • 　このコードのポイントは'grs80 = pyproj.Geod(ellps='GRS80')'です。Gemini先生に解説していただきます。

    grs80 = pyproj.Geod(ellps='GRS80') というコードは、Python の pyproj ライブラリを使用して、GRS80 楕円体に基づく地理座標計算オブジェクトを作成するものです。それぞれの部分について説明します。

    pyproj ライブラリ

    • pyproj は、Python で地理座標変換や測地線計算を行うためのライブラリです。
    • 地球の形状を考慮した高精度な計算が可能になります。

    pyproj.Geod() クラス

    • pyproj.Geod() は、測地線計算を行うためのクラスです。
    • 地球の楕円体モデルに基づいて、2点間の距離や方位角などを計算できます。

    ellps='GRS80' 引数

    • ellps は、楕円体（ellipsoid）モデルを指定する引数です。
    • 'GRS80' は、GRS80（Geodetic Reference System 1980）楕円体モデルを指定しています。
    • GRS80 楕円体は、地球の形状を近似する回転楕円体であり、測地学における基準として広く使用されています。
    • この引数を指定することで、GRS80 楕円体に基づいた計算が行われるようになります。

    grs80 変数

    • grs80 は、作成された pyproj.Geod オブジェクトを格納する変数です。
    • このオブジェクトを使用して、距離計算や方位角計算などの測地線計算を行うことができます。

    まとめ

    このコードは、GRS80 楕円体モデルを使用して、地球上の2点間の距離や方位角などを計算するための準備を行っています。作成された grs80 オブジェクトを後続のコードで使用することで、高精度な地理座標計算が可能になります。

    　このgrs80オブジェクトを使用してaz, bkw_az, dist = grs80.inv(lon1, lat1, lon2, lat2)''というコードで以下の数値を計算しています。pyprojライブラリはおそるべしでございます。

  • az: 1点目から2点目への方位角（Azimuth）。角度単位（度）。
  • bkw_az: 2点目から1点目への逆方位角（Back Azimuth）。角度単位（度）。

  • dist: 2点間の距離（Distance）。メートル単位。

    • 　GRS80楕円体についても聞いてみました。

    GRS80楕円体（Geodetic Reference System 1980）は、地球の形状を近似する回転楕円体であり、測地学における基準として広く使用されています。以下に、GRS80楕円体に関する詳細な情報を提供します。

    GRS80楕円体の特徴

    • 形状の定義:
      • 地球を回転楕円体としてモデル化し、その形状を数学的に定義しています。
      • 長半径（赤道半径）と扁平率によって、楕円体の大きさと形状が決定されます。
    • 主なパラメータ:
      • 長半径 (a): 6,378,137 メートル
      • 扁平率 (f): 1/298.257222101
    • 測地学における役割:
      • 地球上の位置を正確に測定し、地図を作成するための基準となります。
      • GPSなどの衛星測位システムや、地理情報システム（GIS）においても、位置情報の基準として使用されます。
    • WGS84との関係:
      • GPSで使用されているWGS84（World Geodetic System 1984）楕円体は、GRS80楕円体と非常に近い値を持っています。
      • 実用上、両者はほぼ同一のものとして扱われることが多いです。

    GRS80楕円体の重要性

    • 地球の形状を正確にモデル化することで、高精度な測量や位置測定が可能になります。
    • 国際的な測地学の基準として、様々な分野で広く利用されています。
    • GPSなどの衛星測位システムを支える重要な基盤となっています。

    GRS80楕円体とWGS84楕円体の違い

    GRS80楕円体とWGS84楕円体は非常に近い値を持っていますが、厳密にはわずかな違いがあります。

    • WGS84は、GPSの運用に合わせて、より精密な測定に基づいて更新されています。
    • そのため、最新の測量においては、WGS84の方がより正確な値を提供します。

    しかし、一般的な用途においては、両者の違いは無視できるほど小さいため、どちらを使用しても実用上の問題はほとんどありません。

    参考資料

    ５．Digital Elevation Modelの標高値を平均化するFunction 'tile_elev_sm1'

    　使用するかどうかは別にして、Tile内部の行列数を変えずに、与えられた引数'number'だけの上下左右の範囲の標高値を平均して、当該行列Pointの標高値を置き換える関数です。

    　地理院DEMは既に平均化されているように思えるので、このFunctionの要否については微妙化もしれませんが、ゆくゆくドローン測量した標高値など誤差を含む可能性があるModelを利用する際には有効かもしれません。


    # Smoothing the DEM data
    ### DEM(標高)データの平滑化
    # DEMを平滑化(1回目)する関数'tile_elev_sm1st()'
    # 当該内部タイルに隣接するnumber個分の標高値を平均
    def tile_elev_sm1(number):
        x_num = tile_num_x * 256 #東西方向の標高値のデータ数
        y_num = tile_num_y * 256 #南北方向の標高値のデータ数
        for y in range(0 , y_num - number):
            if y <= number - 1:
                continue
            for x in range(0, x_num - number):
                if x <= number -1:
                    continue
                total = 0.0
                for n in range(-1 * number, number + 1):
                    for m in range(-1 * number, number + 1):
                        total = total + float(tile[y + n][x + m])
                sm1_tile[y][x] = total / ((2 * number + 1)**2)


    ６．n×n個のタイルの標高値を平均化して、１個の標高値にまとめるコード

    　このコードについてはFunctionにすることを失念しており、MainCodeのままでお届けします。何かの時に使うことがあるかもしれないとCreatしたコードで、VariousなCalculateをするときに、DEMのDataが多すぎると処理時間が膨大になる可能性があると思って、n×n個のタイルの標高値を平均化したうえで、１個の標高値にまとめてしまおうというCodeです。

    　このコードを通すことで、例えば、256×256の標高値を持つTileを8×8の領域で読み込んだら、大変な量の標高値を持つModelになりますが、8×8個のタイルの標高値を平均化して１つにまとめることで256×256の標高値モデルに簡略化できるということになります。


    #集水面積計算用　集約タイル作成
    #各種変数の初期化
    ca_tile = []; temp_tile = []
    x = 0; y = 0; x_num = 0; y_num = 0
    ca_x_num = 0; ca_y_num = 0
    x_num = tile_num_x * 256 #東西方向の標高値のデータ数
    y_num = tile_num_y * 256 #南北方向の標高値のデータ数
    ca_x_num = int(x_num / ca_tile_pitch) #東西方向の集水面積計算用標高値のデータ数
    ca_y_num = int(y_num / ca_tile_pitch) #南北方向の集水面積計算用標高値のデータ数
    print("ca_x_num, ca_y_num", ca_x_num, ca_y_num)

    for y in range(0, y_num, ca_tile_pitch):
        ca_tile_y = int(y / ca_tile_pitch)
        temp_tile = []
        for x in range(0, x_num, ca_tile_pitch):

            ca_tile_x = int(x / ca_tile_pitch)
            total = 0.0
            for n in range(0, ca_tile_pitch):
                for m in range(0, ca_tile_pitch):
                    total = total + float(sm1_tile[y + n][x + m])
    # タイル内部座標に対して東西・南北方向に'ca_tile_pitch'分の標高値を平均した数値を持つ集水面積計算用のリストを生成
            temp_tile.append(total / ca_tile_pitch **2)
        ca_tile.append(temp_tile)

    （補足）前述５．、６．のCodeをWorkさせるためには、以下のようなInitial Inputが必要になります。

    #################################
    ##### DEMデータの平滑化条件入力 #####
    #### Code No.2 use this data ####
    #################################
    ## 平滑化の対象とする隣接タイル数の入力
    # 20250208_Not adopt smoothing elevation data
    # 1回目平準化に用いる平準化対象内部タイル数'sm1_range'
    sm1_range = 0
    # 集水面積計算用のタイル内部座標ピッチ'ca_tile_pitch'
    ca_tile_pitch = 1 # Default = 8


    ７．Confirmation(確認)

    　四国地方のMajor MountainsのLongitude、LatitudeをGoogleMapでGetし、そのLongitude、Latitudeを入力値として、DEM(GSITile)のElevationValueをObtainして、OfficialなElavation of the Top of MountainとCompereするCodeです。

    　さらに、ComareResultについては、Foliumを使ってGSIがPublishしているStandardMap上にMarkerをPlotして、DEM標高値とOfficialElevationをPopupでShowupしています。

    import folium

    # 地図の中心座標（四国の中央付近）
    center_lat = 33.7
    center_lon = 133.3
    # 地理院地図（標準地図）のタイルURL
    tile_url = "https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png"

    # 地図の作成
    map = folium.Map(location=[center_lat, center_lon], zoom_start=9, tiles=tile_url,
        attr='<a href="https://maps.gsi.go.jp/development/ichiran.html" target="_blank">地理院タイル</a>'
    )

    #point_list 山の名称、公式標高値、経度緯度
    point_list = [
        ["石鎚山天狗岳", 1982, [133.11516, 33.76771]],
        ["剣山山頂", 1955, [134.09408, 33.85358]],
        ["稲村山山頂", 1506, [133.35912, 33.74710]],
        ["大滝山山頂", 946, [134.12712, 34.12331]],
        ["伊予富士山頂", 1756, [133.24812, 33.78806]],
        ["飯野山(讃岐富士)山頂", 421.9, [133.84571, 34.27428]]
    ]

    # point_listを順番に読み込むforループ
    for point in point_list:
     
        lonlat = point[2]
        col, row = func_lonlat_to_colrow(zoom, lonlat, nw_edge_lonlat)
        print(point[0], "経度・緯度: ", lonlat, "Tile内部座標(列・行番号): ", col, row)
        print("公式標高値: ", point[1], "m, DEM標高値: ", int(tile[row][col]), \
              "m, 集水面積計算用標高値: ", int(ca_tile[int(row/ca_tile_pitch)][int(col/ca_tile_pitch)]), "m")
        print("")
       
        #Creat the Marker on the Map
        popup=folium.Popup(f"{point[0]}<br> 公式標高値: {point[1]}m<br> DEM標高値: {int(tile[row][col])}m", max_width=300)

        folium.Marker(
            location = [lonlat[1], lonlat[0]],
            popup = popup,
            icon = folium.Icon(color="green"),
        ).add_to(map)

    map

    まずは、数値のOutputです。

    石鎚山天狗岳 経度・緯度: [133.11516, 33.76771] Tile内部座標(列・行番号): 1351 948 公式標高値: 1982 m, DEM標高値: 1956 m, 集水面積計算用標高値: 1956 m 剣山山頂 経度・緯度: [134.09408, 33.85358] Tile内部座標(列・行番号): 2776 797 公式標高値: 1955 m, DEM標高値: 1942 m, 集水面積計算用標高値: 1942 m 稲村山山頂 経度・緯度: [133.35912, 33.7471] Tile内部座標(列・行番号): 1706 984 公式標高値: 1506 m, DEM標高値: 1493 m, 集水面積計算用標高値: 1493 m 大滝山山頂 経度・緯度: [134.12712, 34.12331] Tile内部座標(列・行番号): 2824 323 公式標高値: 946 m, DEM標高値: 933 m, 集水面積計算用標高値: 933 m 伊予富士山頂 経度・緯度: [133.24812, 33.78806] Tile内部座標(列・行番号): 1544 912 公式標高値: 1756 m, DEM標高値: 1736 m, 集水面積計算用標高値: 1736 m 飯野山(讃岐富士)山頂 経度・緯度: [133.84571, 34.27428] Tile内部座標(列・行番号): 2414 57 公式標高値: 421.9 m, DEM標高値: 414 m, 集水面積計算用標高値: 414 m

    全ての結果において、Offcial Elevation＞DEM Elevationになっています。これはわたしの想像ですが、OfficialElevationは山頂の点の数値ですが、DEM標高値は山頂周辺の平均的な標高になっているのではないかなぁと思います。

    FoliumによるPlotのResultです。最初は四国全域が表示され、指定したお山付近にマーカーが表示されていると思います。



    石鎚山天狗岳付近をZoomUPしてみるとこんな感じです。マーカーはほぼ天狗岳山頂付近に位置しており、経度・緯度はほぼ正確にGetできているようです。



    コメント

    コメントを投稿