～～～作成中～～～
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 = lola_to_dist(lon1, lat1, lon2, lat1) / (256 * (es_x - wn_x)) # 距離刻み
lat_pich = lola_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
前述の３．のReverceの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
コメント
コメントを投稿