For Not only Japan domein, but also global dataでTerrain AnalysisをActするべく、Geospatial Information Authority of Japan（国土地理院）のDEM_Dataから離れて、JAXAがDevelop&OperateしたArtificial Steliteのだいち（ALOS）がGetした高精度の全天球ElevationModelの活用についてLearningでございます。
（参考）JAXAのだいち（ALOS）に関する解説Page：https://www.satnavi.jaxa.jp/ja/project/alos/index.html
だいち（ALOS）にはたくさんの機能があると、改めて思いますが、ここはSimpleにElevationDataの活用を目指します。このWebSiteはALOSをはじめとするJAXAのArtificialSteliteのIntroductionで、ALOSのVarius kind of DataをGettingは次のURLになります。
（参考）ALOS利用推進研究プロジェクトのURL：https://www.eorc.jaxa.jp/ALOS/jp/index_j.htm
このPageにもいろいろなDataがexistしますが、今回は「全球高精度デジタル３Ｄ地図(ALOS_World_3D)」を活用させていただきます。このDataの概要は同Pageによると次のように解説されています。「「だいち」搭載光学センサPRISMを用いた全球陸域を対象とした高精度デジタル3D地図。水平解像度30m相当の全球数値標高モデル第3.2版を公開中。」
ALOS_World_3DのURL（https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d_j.htm）に移動すると、ProfessinalSpecの平面解像度5mVersionと、PublicSpec&FreeOfChargeの平面解像度30mVersionがIntroduceされています。今回はLearning＆Practiceなので、PriceFreeの30mPressisionVersionをMakeUseします。このPageに「全球数値地表モデル (DSM) の30m解像度版データセット」へのリンク（https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d30/aw3d30_j.htm）があるので、こちらに移動すると、Like BelowのGlobalMapが表示されます。
LandRegionのAllOfAreaがColorDrawnなので、ver4.1までに全世界のElevationDataがBeAbleToCompleteのRecognitionです。すごいのが、Commercial or Non-Commercial問わずにFreeでAvailableというPointでございます。JAXAさまに大感謝でございます。
わたくしはAmatureUserなので、早速、「４．データのダウンロード方法」と「ユーザー登録」しました。ユーザー登録は指示に従ってProcessすれば問題なく完了するはずです。
UserRegisterationが完了したら、早速データのダウンロード用のリンク先（https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm）に移動します。
最初はPythonのコードから直接Download Siteにアクセスして、圧縮ファイルを取得しようとおもったのですが、UserIDとpassの入力が必要なためWeb scripingなどの操作が煩雑になるため、一度、必要な圧縮ファイルをダウンロードして、それをPythonCodeで解凍した方がよいとGemini先生に助言をいただいたので、その方法を採用しました。
まずは手始めに中国、四国、九州地方を含む５✕５のTileをDownLoadします。
話題が前後しますが、ALOSのTile座標の考え方をGemini先生に確認すると、おおむね以下のような回答をいただきました。
- ALOSのTileの大きさは緯度１度、経度１度単位で設定
- Tileの中の標高値は東西3600×南北3600個を格納
- Tile中のPixel座標は東西・南北ともに緯度・経度と線形関係で設定（標高値の平面的感覚は緯度・経度ともに１秒ごとに設定されている）
Pixel座標と緯度・経度の変換式は、具体的には別途学習したいと思いますが、Geospatial Information Authority of Japanのいわゆる地理院Tileと比較すると、複雑ではなさそうです。
話題を戻して、先ほどのDataDownloadSiteに進み、ID/Passを入力すると次の画面が表示されました。
このMap上でDownloadしたいRegionをLeftCrickしていくと、MapがEnlargeしていきます。数回LeftClickすると下のような画面になります。これが５×５のTileの一括Download画面になり、オレンジ色のDownloadButtonをClickすると、DataのDownloadが始まります。この画面（四国地方全域と、中国・九州地方の一部）のDataは「N030E130_N035E135.zip」というNameのZipFileでしたBottomLeftのTileの左上角部の座標が北緯30度・東経130度を、UpperRightの左上角部の座標が北緯35度・東経135度であることを表してるNameであると予想しました。
ちなみに、この５×５のTileのいずれか１つのTileをClickすると、１つのTileのみをDownloadする以下のような画面にTransferします。
おおむね準備が整いました。「N030E130_N035E135.zip」を解凍して、格納されている標高DataをShowPlotするPythonCodeをGemini先生のAdviceをいただきながら、Creatしていきます。
下準備として、Downloadした「N030E130_N035E135.zip」をGoogleColabのTemporaryFolderにDrag&Dropして、Uploadします。
そして、GoogleColabのPythonCodeで「N030E130_N035E135.zip」を読み込んで、解凍します。Gemini先生のOpusです。Gemini先生はtry-exceptを多用してCodeのVolumeを増やすくせがあるようで、最近少々困っています。
import zipfile
import os
# アップロードしたzipファイルのパス
zip_filepath = '/content/N030E130_N035E135.zip'
# 解凍先のベースディレクトリ
# この中に 'N030E130_N035E135' フォルダが作成され、GeoTIFFファイルが格納されます
extracted_base_dir = '/content/extracted_elevation_data'
# 解凍先の最終パスを生成 (zipファイル内のフォルダ構造を考慮)
# zipファイルが N030E130_N035E135 というフォルダを含んでいると仮定
final_extracted_data_dir = os.path.join(extracted_base_dir, 'N030E130_N035E135')
# 解凍先のディレクトリが存在しない場合は作成
os.makedirs(final_extracted_data_dir, exist_ok=True)
print(f"Zipファイル '{zip_filepath}' の解凍を開始します...")
try:
with zipfile.ZipFile(zip_filepath, 'r') as zip_ref:
# extractall() は、zipファイル内のトップレベルのフォルダ構造を維持して解凍します。
# そのため、extracted_base_dir の中に N030E130_N035E135 フォルダが作成されます。
zip_ref.extractall(extracted_base_dir)
print(f"Zipファイルを '{extracted_base_dir}' に正常に解凍しました。")
print(f"GeoTIFFファイルは '{final_extracted_data_dir}' の中に格納されています。")
except FileNotFoundError:
print(f"エラー: Zipファイル '{zip_filepath}' が見つかりません。ファイル名とパスを確認してください。")
except zipfile.BadZipFile:
print(f"エラー: '{zip_filepath}' は無効なzipファイルです。再ダウンロードを検討してください。")
except Exception as e:
print(f"Zipファイルの解凍中に予期せぬエラーが発生しました: {e}")
ZipFileの回答が成功すると、以下のようなMassegeが出力されます。
Zipファイル '/content/N030E130_N035E135.zip' の解凍を開始します...
Zipファイルを '/content/extracted_elevation_data' に正常に解凍しました。
GeoTIFFファイルは '/content/extracted_elevation_data/N030E130_N035E135' の中に格納されています。
スムーズにFileを解凍できた模様です。FileのA‗Part Of Contentsはこのようになっています。標高データは「…_DSM.tif」に格納されており、FileNameの後半のN030E130は当該TileのUpperRightCornerが北緯30度東経130度であることを表しています。
最初に「ALPSMLC30_N031E130_DSM.tif」という１つのTileのElevationDataをMatplotlibを使ってContourをCreatしてみます。Gemini先生にCodeをお願いしました。
import rasterio
import matplotlib.pyplot as plt
import numpy as np
# TIFFファイルのパス
tif_file_path = '/content/extracted_elevation_data/N030E130_N035E135/ALPSMLC30_N031E130_DSM.tif'
print(f"TIFFファイル '{tif_file_path}' の読み込みを開始します...")
try:
# rasterio を使ってTIFFファイルを読み込む
with rasterio.open(tif_file_path) as src:
# 標高データを読み込む (通常は1番目のバンド)
elevation_data = src.read(1)
# ファイルのメタデータから座標参照系 (CRS) や変換情報 (transform) を取得
crs = src.crs
transform = src.transform
bounds = src.bounds
print(f"ファイル名: {src.name}")
print(f"幅 (pixels): {src.width}")
print(f"高さ (pixels): {src.height}")
print(f"バンド数: {src.count}")
print(f"CRS (座標参照系): {crs}")
print(f"バウンディングボックス (minx, miny, maxx, maxy): {bounds}")
print(f"データ型: {elevation_data.dtype}")
print(f"データ範囲 (min, max): {np.nanmin(elevation_data)}, {np.nanmax(elevation_data)}")
print("データ読み込み完了。図化を開始します...")
# Matplotlibで図化
plt.figure(figsize=(10, 8)) # 図のサイズを設定
# imshowを使ってデータを画像として表示
# extent: データの地理的な範囲 (left, right, bottom, top)
# cmap: カラーマップ (標高データには'terrain', 'viridis', 'gray'などがよく使われます)
# origin='upper': 画像の原点を左上にする (地理空間データでは一般的)
# interpolation='none': ピクセル補間なし (データがそのまま表示される)
img_plot = plt.imshow(
elevation_data,
cmap='terrain', # 'viridis', 'gray', 'jet' なども試してみてください
origin='upper',
extent=[bounds.left, bounds.right, bounds.bottom, bounds.top]
)
# カラーバーを追加 (標高の凡例)
plt.colorbar(img_plot, label='Elevation (meters, assumed)')
# タイトルと軸ラベルを設定
plt.title(f'Elevation Data from {tif_file_path.split("/")[-1]}')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
# グリッドを追加 (必要に応じて)
plt.grid(True)
# 図を表示
plt.show()
print("\n図化が完了しました。")
except rasterio.errors.RasterioIOError as e:
print(f"エラー: TIFFファイルの読み込み中に問題が発生しました。ファイルが存在するか、破損していないか確認してください。\n詳細: {e}")
except Exception as e:
print(f"予期せぬエラーが発生しました: {e}")
BelowはOutputです。このShapeからThinkすると、おそらくKagoshimaPref.で左上の高標高が韓国岳や開聞岳、霧島連山で、中央付近のIslandは桜島ではないかとImageしました。ここまではSimpleなOperationです。
次は、Some Tileを使って広域なDEMDataBaseをCreatしていきます。
コメント
コメントを投稿