Thema: Download the Digital Elevation Model(DEM) data from GSI(Geospatial Information Authority )in Japan (Python and GeoPandas)
GIS(Geographic Information System)の学習の一環として、国土地理院(GSI)のWeb SiteからDEM(Digital Elevation Model)をDownloadして、mutplotlibでPlotするCodeを作成してみました。
DEMの活用はGIS習得のうえでBaseになる技術だと考えています。
DEMは国土地理院(GSI)の電子国土Webのデータを利用させていただきます。国土地理院の電子国土Webデータの紹介はこちらです。
https://maps.gsi.go.jp/help/intro/ （国土地理院Webページ）
国土地理院の電子国土Webデータは国土地理院タイルという形式で整理されています。わたしの解釈では、ZoomLevelという概念を使って地球地図全体を正方形の複数個(Zoom×Zoom)のタイルで分割し、１つ１つのTileを256×256のMeshに分割して、そのMeshごとにElevationValueを格納しているという感じです。国土地理院タイルの説明はこちらです。
https://maps.gsi.go.jp/development/siyou.html（国土地理院Webページ）
国土地理院の電子国土Webを使って四国周辺の国土地理院標準地図を表示すると以下のような表示になります。これは「設定」で「タイル座標の表示」をonにSetUpした表示です。
地図内が赤い線で区切られているのが地理院タイルであり、「9/444/204」などの表示が地理院タイルのZoomLevelとTileCoordinateを表現しています。ちなみに、一番最初の数値(ここでは「9」）はZoomLevel、二番目の数値（ここでは「444」）はEast-West方向のTileCoordinate、三番目の数値（ここでは「204」）はSouth-West方向のTileCoordinateを表しています。
ZoomLevelが「1」ならば、地球全体を１つのタイルで、すなわち、256×256個の標高値で、また、ZoomLevelが「10」ならば、地球全体を100個のタイル、すなわち25600×25600個の標高値で表現することになり、ZoomLevelが大きくなるほど標高値の分解能が上昇し、これに伴ってデータ処理時間が大きくなっていきます。
ZoomLevelを'z'、East-West方向、North-South方向のTileCoordinateをそれぞれ’x'、’y'として、任意の１つのTile(z,x,y)を呼び出すためのPythonのFunction'func_fetch_tile'は以下のようになります。
import pandas as pd
# Fuction of get the DEM data of individial tile 'func_fetch_tile()' / 個別タイルのDEM(標高)データを取得する関数
# z:zoom, x: x coordinate(EW direction) of the tile, y: y coordinate(NS direction) of the tile
def func_fetch_tile(z, x, y):
# Enter the location on the WEB site of GIS Japan that has the DEM data / urlに国土地理院のHPからDEM(標高)データの場所を入力
url = "https://cyberjapandata.gsi.go.jp/xyz/dem/{z}/{x}/{y}.txt".format(z=z, x=x, y=y)
次に、当該Tileに格納されている256×256個の標高値csvをPadasDataFrameとしてreadします。同時に、国土地理院DEMでは海域の標高値としてSetUpされている'e'(empty?)について、数値処理上で不都合が発生するケースが多いため、'0'に置き換えます。
# Enter the Getting DEM data to PandasDataFrame'df', if directed URL and filename exist,/ 指定したURL,FileNameが存在する場合は、変数dfにpandasDataFrameとしてDEMデータを入力し、
# and replace 'e' values in original DEM data to '0' values / かつ、標高データに'e'が存在する場合は'0'に置き換え
try:
df = pd.read_csv(url, header=None).replace('e', 0)
次に、陸域が存在しないTile（例えば、AboveMapでいえば「9/446/206」や「9/447/206」など）にはDataそのものがExistしていないようなので、256×256個の'0'値を仮SetUpしたダミーのTileをCreateします。
from urllib.error import URLError, HTTPError
# As creating the tentitive PandasDataFrame for the tile that exist only ses area / 海上のみのタイルに仮のDEMデータを作成するため、
# Enter '0' values to all of df (column 256 * row 256) to PandasDataFrame'df', / 変数dfに全て'0'を入力(DEMデータは東西256×南北256のメッシュで構成)
# if the directed URL and filename don't exist(in the case of entire tile on the sea) / 指定したURLやFileNameが存在しない場合（例えば、指定したタイル座標が全て海上の場合）
except URLError as e:
temp_tile = []
temp_row = []
zero = 0.0
for i in range(0, 256):
temp_row.append(zero)
for i in range(0, 256):
temp_tile.append(temp_row)
df = pd.DataFrame(temp_tile)
except HTTPError as e:
temp_tile = []
temp_row = []
zero = 0.0
for i in range(0, 256):
temp_row.append(zero)
for i in range(0, 256):
temp_tile.append(temp_row)
df = pd.DataFrame(temp_tile)
return df.values
Functionの最後に作成した当該Tile１つ分のPandasDataFrame'df'を'df.value’でNumpy ndarrayの二次元配列として返します。
次に、複数のTileをCombineするFunction'func_fetch_all_tiles'をDefineします。Function’func_fetch_all_tiles’に渡す引数’north_west', 'south_east'は後述しますが、検討対象領域の北西端および南東端の(Zoom,X,Y)のことです。まず、引数’north_west', 'south_east'の一つ目すなわちZoomLevelがdifferentのCaseはErrorMassegeをOutputします。そして、Tile座標のx,y方向それぞれの個数をx_range,y_rangeとして計算し、Nampy ndarrayを結合するnp.concatenateを使用して、1つずつFunction’func_fetch_tile’をつかって１つずつTileを呼び出し、Combineして、検討対象領域の複数のTileを１つのModelにまとめます。
# function 'func_fetch_all_tiles' combine the multiple DEM tiles from Area of consideration that be difened by the tile numbers of North-West, South-East edge, North-East, South-West edges /
# 北西端・南東端タイル座標を指定した矩形領域の複数タイルのDEM(標高)データを結合する関数'func_fetch_all_tiles'
def func_fetch_all_tiles(north_west, south_east):
assert north_west[0] == south_east[0] , "タイル座標のzoomが一致していません"
x_range = range(north_west[1], south_east[1] + 1)
y_range = range(north_west[2], south_east[2] + 1)
return np.concatenate(
[
np.concatenate(
[func_fetch_tile(north_west[0], x, y) for y in y_range],
axis = 0
) for x in x_range
],
axis = 1
)
Two Functions Above I write, 検討対象領域の複数のDEMを１つのmodelにまとめるmainCodeは以下のとおりです。
# Calculate the number of South-East edge tile / 南東端のタイル番号の計算
south_east_x = north_west_x + tile_num_x -1 #x coordinate(number) South-East edge tile
south_east_y = north_west_y + tile_num_y -1 #y coordinate(number) South-East edge tile
# Combine a lot of tile in the area of consideration to 'tile'
tile = func_fetch_all_tiles((zoom,north_west_x, north_west_y), (zoom, south_east_x, south_east_y))
たったこれだけです。Function'func_fetch_all_tiles'に渡す引数は(zoom,north_west_x, north_west_y)のように3つのValueを持っています。
これだけでは何が起きているのかよくわからないので、作成したモデルをContour(等高線)にしてPlotします。
import numpy as np
import matplotlib.pyplot as plt
# サンプルデータ作成 (適宜変更してください)
grid_x = tile_num_x * 256
grid_y = tile_num_y * 256
X, Y = np.meshgrid(np.linspace(grid_x, 0, grid_x), np.linspace(grid_y, 0, grid_y))
# FigureとAxesを作成 (3行4列)
fig, axes = plt.subplots(3, 4, figsize=(16, 8), dpi=80, facecolor='w', edgecolor='k')
# 表示するcmapのリスト (お任せ)
cmaps = ['binary', 'viridis', 'spring_r', 'summer_r', 'autumn_r', 'winter_r', 'Reds', 'hot', 'cool', 'PuOr', 'RdYlGn', "twilight"]
GoogleGeminiの力も借りながら、CodeのOverviewをExplainしていきます。
◆np.linspace(grid_x, 0, grid_x)
np.linspace() は、指定された範囲で等間隔な数値を生成する関数です。この部分では、grid_x から 0 までの範囲を grid_x 個に分割した等間隔な数値を生成しています。生成された配列は、X軸方向の座標を表します。
◆np.meshgrid(...)
np.meshgrid() は、2つの1次元配列（ここではX軸方向とY軸方向の座標）を引数として受け取り、2次元のグリッド座標を生成します。具体的には、X軸方向の配列をY軸方向にコピーし、Y軸方向の配列をX軸方向にコピーすることで、すべてのグリッドポイントの座標値を持つ2つの2次元配列（XとY）を生成します。
◆fig, axes = plt.subplots(3, 4, figsize=(16, 8), dpi=80, facecolor='w', edgecolor='k')
このコードは、PythonのMatplotlibライブラリを使って、複数のグラフ（サブプロット）を1つのFigure内にグリッド状に配置するためのものです。
◆cmaps = ['binary', 'viridis', 'spring_r', 'summer_r', 'autumn_r', 'winter_r', 'Reds', 'hot', 'cool', 'PuOr', 'RdYlGn', "twilight"]
このコードは、PythonのMatplotlibライブラリで使用されるカラーマップ（colormap）の名前を格納したリストを作成するものです。
以下、Matplotlibライブラリを使って、複数のサブプロット（グラフ）を作成し、それぞれのサブプロットに異なるカラーマップを適用して等高線図を描画するものです。
# 各サブプロットでContourを描画
for i, cmap_name in enumerate(cmaps): # 4つのサブプロット
row = i // 4 # 行番号を修正
col = i % 4 # 列番号を修正
# tile変数がサポートされていないデータ型の場合、float型に変換 (各サブプロットで個別に行う)
if not np.issubdtype(tile.dtype, np.number):
tile = tile.astype(np.float64)
# 等高線を描画
elevation = range(0, 2000, 100)
cont = axes[row, col].contour(tile, levels=elevation, cmap=cmap_name) # axesを指定
# Y軸を反転
axes[row, col].invert_yaxis()
# カラーバーを追加 (各サブプロットに個別に追加)
cb = fig.colorbar(cont, ax=axes[row, col], shrink=0.5, aspect=10)
# グラフタイトルを設定
axes[row, col].set_title(cmap_name) # cmap名をタイトルに設定
plt.tight_layout() # レイアウト調整
plt.show()
このCodeに初期入力値として、四国地方全体をmodel化する前提で、以下を入力します。
zoom = 11 #zoomを入力
### 四国全域
north_west_x = 1776 #北西端tileのx座標を入力 基準点
north_west_y = 816 #北西端tileのy座標を入力 基準点
tile_num_x = 15 #東西方向のtile数を入力
tile_num_y = 12 #南北方向のtile数を入力
Varius Kind of Color Mapの四国地方のContour(等高線)が表示されます。
コメント
コメントを投稿