前回：https://shikuuk.blogspot.com/2025/07/jaxaalospython.html
Print文でDiscribeしたMosicedModelのContentsは以下のとおりです。
In relatively larger scale, The summits of mailly mountain locate almost exact location.
続く。
第２部はさきほどDownloadした「N030E130_N035E135.zip」に含まれている５×５のTileを１枚の大きなModelに結合する処理（モザイク処理と呼ぶそうです）のCodeをGemini先生にお願いしました。
import rasterio
import numpy as np
import os
from rasterio.enums import Resampling
from rasterio.merge import merge # Ensure merge is imported
def merge_elevation_tiles(tile_dir, output_filepath):
"""
Reads all GeoTIFF elevation tiles in the specified directory,
merges them into a single GeoTIFF file, and saves it.
Args:
tile_dir (str): Path to the directory containing GeoTIFF files.
output_filepath (str): Output path for the merged GeoTIFF file.
Returns:
str: Path to the merged GeoTIFF file (on success), or None (on failure).
"""
src_files_to_mosaic = []
# Iterate through all files in the directory to find GeoTIFFs
for filename in os.listdir(tile_dir):
if filename.endswith('.tif'): # Only process .tif files
filepath = os.path.join(tile_dir, filename)
try:
src_files_to_mosaic.append(rasterio.open(filepath))
except rasterio.errors.RasterioIOError as e:
print(f"Warning: Could not open {filepath} (possibly corrupted or not a valid GeoTIFF): {e}. Skipping.")
except Exception as e:
print(f"Warning: Unexpected error opening {filepath}: {e}. Skipping.")
if not src_files_to_mosaic:
print("No valid GeoTIFF tile files found to merge in the specified directory. Exiting.")
return None
print(f"Merging {len(src_files_to_mosaic)} tiles...")
# Merge the datasets; resampling='nearest' for discrete elevation values
mosaic, out_transform = merge(src_files_to_mosaic, resampling=Resampling.nearest)
# Convert to float32 if needed to save memory (from float64)
if mosaic.dtype == np.float64:
mosaic = mosaic.astype(np.float32)
print("Converted mosaic to float32 to save memory.")
# Copy metadata from the first source file and update with merged properties
out_meta = src_files_to_mosaic[0].meta.copy()
out_meta.update({
"driver": "GTiff",
"height": mosaic.shape[1], # New height of the merged mosaic
"width": mosaic.shape[2], # New width of the merged mosaic
"transform": out_transform, # New georeferencing transform
"crs": src_files_to_mosaic[0].crs, # Coordinate Reference System (assumed consistent)
"nodata": src_files_to_mosaic[0].nodata # NoData value
})
# Close all source files to release resources
for src in src_files_to_mosaic:
src.close()
# Save the merged mosaic to a new GeoTIFF file
try:
with rasterio.open(output_filepath, "w", **out_meta) as dest:
dest.write(mosaic)
print(f"Successfully merged tiles to '{output_filepath}'.")
return output_filepath
except Exception as e:
print(f"Error saving merged GeoTIFF: {e}")
return None
# --- Main Execution Block ---
if __name__ == "__main__":
# Define the directory where GeoTIFF files are extracted
extracted_data_dir = '/content/extracted_elevation_data/N030E130_N035E135'
# Define the output path for the merged GeoTIFF file
output_merged_filepath = os.path.join(extracted_data_dir, 'merged_5x5_elevation.tif')
print("Starting GeoTIFF tile merging process...")
merged_file_path = merge_elevation_tiles(extracted_data_dir, output_merged_filepath)
if merged_file_path:
print(f"\nSuccessfully merged all tiles. Now analyzing statistics for: {merged_file_path}")
try:
with rasterio.open(merged_file_path) as src:
# Read the entire merged data
merged_elevation_data = src.read(1)
# Convert NoData values to NaN
if src.nodata is not None:
merged_elevation_data[merged_elevation_data == src.nodata] = np.nan
# Get valid (non-NaN) elevation data for statistics
valid_merged_elevation_data = merged_elevation_data[~np.isnan(merged_elevation_data)]
if valid_merged_elevation_data.size > 0:
print(f"Merged data original pixel size: {src.height}x{src.width} pixels")
print(f"Merged data geographical bounds (Lon/Lat): Left={src.bounds.left:.2f}, Bottom={src.bounds.bottom:.2f}, Right={src.bounds.right:.2f}, Top={src.bounds.top:.2f}")
print(f"Merged data elevation range: Min={np.min(valid_merged_elevation_data):.2f}m, Max={np.max(valid_merged_elevation_data):.2f}m")
print(f"Merged data mean elevation: {np.mean(valid_merged_elevation_data):.2f}m")
print(f"Merged data standard deviation: {np.std(valid_merged_elevation_data):.2f}m")
print(f"Merged data 2nd percentile: {np.percentile(valid_merged_elevation_data, 2):.2f}m")
print(f"Merged data 98th percentile: {np.percentile(valid_merged_elevation_data, 98):.2f}m")
else:
print("No valid elevation data found in the merged file for statistics.")
except Exception as e:
print(f"An error occurred during statistics analysis of the merged file: {e}")
else:
print("GeoTIFF file merging failed, cannot proceed with statistics analysis.")
何度かGemini先生にCodeを書き直してもらっていますが、それほど苦労せずにモザイク処理ができました。このコードの大まかな流れは以下の通りと解釈しました。
- mainで5*5のTileが格納されているDirectoryと、Marge（結合＝モザイク処理）したModelを格納するFileNameをassignしたあと、Function‘merge_elevation_tiles’を呼び出す。
- Function‘merge_elevation_tiles’では、Tileが格納されているDirectory内のFileNameを全て読み込み、その中から拡張子が‘.tif’（Geotiffファイル）をSelectし、DirectoryとFileNameを合成したStringを作成し、それらのStringをList‘src_files_to_mosaic’にAppendする。
- ‘src_files_to_mosaic’に格納されたFileListのすべてを使って、‘mosaic, out_transform = merge(src_files_to_mosaic, resampling=Resampling.nearest)’でmosaic画像と地理参照情報out_transformを作成する。このOneArticleでTileの結合ができるというのはすさまじいConvenientですね。
このCodeのOutputはLikeBelowです。PixelSizeが18000*18000なので3600*3600のTileが5*5個分結合されているらしいのですが、ElavationRangeがFrom0mTo255.0mというのは西日本最高峰の石鎚山が1900m級であることをThinkするとCouldNotUnderstandなResultです。
Starting GeoTIFF tile merging process...
Merging 54 tiles...
Successfully merged tiles to '/content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif'.
Successfully merged all tiles. Now analyzing statistics for: /content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif
Merged data original pixel size: 18000x18000 pixels
Merged data geographical bounds (Lon/Lat): Left=130.00, Bottom=30.00, Right=135.00, Top=35.00
Merged data elevation range: Min=0.00m, Max=255.00m
Merged data mean elevation: 13.94m
Merged data standard deviation: 41.39m
Merged data 2nd percentile: 0.00m
Merged data 98th percentile: 192.00m
というわけで、結合したFile‘merged_5x5_elevation.tif’をMatplotlibでDrawImageするCodeをGemini先生にAskToCreatしました。
### 20250728_Latest Version
import rasterio
import matplotlib.pyplot as plt
import numpy as np
import os
def plot_elevation_tile(filepath):
"""
指定されたGeoTIFF標高ファイルをMatplotlibで表示します。
Args:
filepath (str): 表示するGeoTIFFファイルのパス。
"""
try:
with rasterio.open(filepath) as src:
# 標高データを読み込む（通常はバンド1）
elevation_data = src.read(1)
# NoData値が存在する場合、それをNaN（非数）に変換してプロット時に無視されるようにする
if src.nodata is not None:
elevation_data[elevation_data == src.nodata] = np.nan
# プロットのセットアップ
plt.figure(figsize=(10, 8)) # 図のサイズを設定
# 画像として標高データを表示
# extent: データの地理的な範囲を設定し、軸のラベルを地理座標にする
# cmap: カラーマップ（標高データに適した'viridis'を使用）
# origin: 画像の原点（GeoTIFFは通常左上が原点）
# vmin, vmax: カラーマップの最小値と最大値を設定し、色の範囲を固定する
# ここでは、データの2パーセンタイルと98パーセンタイルを使用
# これにより、外れ値がカラーマップの表示範囲を歪めるのを防ぐ
valid_elevation_data = elevation_data[~np.isnan(elevation_data)]
if valid_elevation_data.size > 0:
vmin = np.percentile(valid_elevation_data, 2)
vmax = np.percentile(valid_elevation_data, 98)
else:
vmin = np.min(elevation_data[~np.isnan(elevation_data)]) if not np.all(np.isnan(elevation_data)) else 0
vmax = np.max(elevation_data[~np.isnan(elevation_data)]) if not np.all(np.isnan(elevation_data)) else 100
# 画像の表示範囲を設定 (左経度, 右経度, 下緯度, 上緯度)
# rasterioのboundsは (left, bottom, right, top) の順
bounds = src.bounds
extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]
im = plt.imshow(
elevation_data,
cmap='viridis', # 標高データに適したカラーマップ
origin='upper', # 画像の原点を左上に設定
extent=extent, # 地理的な範囲を設定
vmin=vmin, # カラーマップの最小値
vmax=vmax # カラーマップの最大値
)
plt.title(f'Elevation Data: {os.path.basename(filepath)}')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.colorbar(im, label='Elevation (meters)') # カラーバーを追加
plt.grid(True, linestyle=':', alpha=0.7) # グリッドを追加
plt.tight_layout() # レイアウトを調整
plt.show()
except rasterio.errors.RasterioIOError as e:
print(f"エラー: GeoTIFFファイル '{filepath}' を開けませんでした: {e}")
except Exception as e:
print(f"予期せぬエラーが発生しました: {e}")
# --- メイン実行ブロック ---
if __name__ == "__main__":
# 前のコードで結合されたGeoTIFFファイルのパス
# このパスは、前のコードの `output_merged_filepath` と一致している必要があります
merged_tile_filepath = '/content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif'
# 結合されたタイルが存在するか確認
if os.path.exists(merged_tile_filepath):
print(f"GeoTIFFファイル '{merged_tile_filepath}' を表示します...")
plot_elevation_tile(merged_tile_filepath)
else:
print(f"エラー: GeoTIFFファイル '{merged_tile_filepath}' が見つかりません。")
print("先に結合処理を実行してファイルが生成されていることを確認してください。")
OutputはLikeBelowです
Actually Full Of DisconfortableなPlotでございます。このあと数日Gemini先生とDiscussionしながら、Cotinue to worryしました。例えば、ConbineしたMosicTileのDataCapacityが大きすぎて、MatplotlibでPlotできないとか、いろいろな理由を探しましたが、最終的に至った結論は、5*5のTileを読む込むときに、必要なTiffFileをLoadできていないことがCauseだとSeemされました。そこで、NecessaryなTiffFileをMyselfでSelectして読み込む方法にReviceします。
先ほどのCodeをGemini先生のPromptにCopy＆Pasteしたうえで、以下のRequestをWriteします。
>input<
このコードだとtifデータを読み込めないみたいなので、コード側で以下のファイル名を指定して、データを読み込むようにしてもらえませんか。
'ALPSMLC30_N034E134_DSM.tif',
'ALPSMLC30_N033E134_DSM.tif',
'ALPSMLC30_N034E133_DSM.tif',
'ALPSMLC30_N033E133_DSM.tif',
'ALPSMLC30_N032E133_DSM.tif',
'ALPSMLC30_N033E132_DSM.tif',
'ALPSMLC30_N034E132_DSM.tif',
'ALPSMLC30_N032E132_DSM.tif',
'ALPSMLC30_N033E131_DSM.tif',
'ALPSMLC30_N034E131_DSM.tif',
'ALPSMLC30_N032E131_DSM.tif',
'ALPSMLC30_N031E131_DSM.tif',
'ALPSMLC30_N034E130_DSM.tif',
'ALPSMLC30_N033E130_DSM.tif',
'ALPSMLC30_N032E130_DSM.tif',
'ALPSMLC30_N031E130_DSM.tif',
'ALPSMLC30_N030E130_DSM.tif'
なお、tifデータが保存されているディレクトリは従来と同じです。
>OutputCode<
###20250729_LatestVersion
###Get the Tiff Data by assigning the file name directly
import rasterio
import matplotlib.pyplot as plt
import numpy as np
import os
from rasterio.enums import Resampling
from rasterio.merge import merge # Ensure merge is imported
def merge_elevation_tiles(tile_dir, output_filepath):
"""
Reads specific GeoTIFF elevation tiles in the specified directory,
merges them into a single GeoTIFF file, and saves it.
Args:
tile_dir (str): Path to the directory containing GeoTIFF files.
output_filepath (str): Output path for the merged GeoTIFF file.
Returns:
str: Path to the merged GeoTIFF file (on success), or None (on failure).
"""
# ユーザーが指定したGeoTIFFファイル名のリスト
specific_tile_filenames = [
'ALPSMLC30_N034E134_DSM.tif',
'ALPSMLC30_N033E134_DSM.tif',
'ALPSMLC30_N034E133_DSM.tif',
'ALPSMLC30_N033E133_DSM.tif',
'ALPSMLC30_N032E133_DSM.tif',
'ALPSMLC30_N033E132_DSM.tif',
'ALPSMLC30_N034E132_DSM.tif',
'ALPSMLC30_N032E132_DSM.tif',
'ALPSMLC30_N033E131_DSM.tif',
'ALPSMLC30_N034E131_DSM.tif',
'ALPSMLC30_N032E131_DSM.tif',
'ALPSMLC30_N031E131_DSM.tif',
'ALPSMLC30_N034E130_DSM.tif',
'ALPSMLC30_N033E130_DSM.tif',
'ALPSMLC30_N032E130_DSM.tif',
'ALPSMLC30_N031E130_DSM.tif',
'ALPSMLC30_N030E130_DSM.tif'
]
src_files_to_mosaic = []
# 指定されたファイル名リストをイテレートしてGeoTIFFを開く
for filename in specific_tile_filenames:
filepath = os.path.join(tile_dir, filename)
if not os.path.exists(filepath):
print(f"Warning: Specified file '{filepath}' does not exist. Skipping.")
continue # ファイルが存在しない場合はスキップ
try:
src_files_to_mosaic.append(rasterio.open(filepath))
except rasterio.errors.RasterioIOError as e:
print(f"Warning: Could not open {filepath} (possibly corrupted or not a valid GeoTIFF): {e}. Skipping.")
except Exception as e:
print(f"Warning: Unexpected error opening {filepath}: {e}. Skipping.")
if not src_files_to_mosaic:
print("No valid GeoTIFF tile files found to merge from the specified list in the directory. Exiting.")
return None
print(f"Merging {len(src_files_to_mosaic)} tiles...")
# Merge the datasets; resampling='nearest' for discrete elevation values
mosaic, out_transform = merge(src_files_to_mosaic, resampling=Resampling.nearest)
# Convert to float32 if needed to save memory (from float64)
if mosaic.dtype == np.float64:
mosaic = mosaic.astype(np.float32)
print("Converted mosaic to float32 to save memory.")
# Copy metadata from the first source file and update with merged properties
out_meta = src_files_to_mosaic[0].meta.copy()
out_meta.update({
"driver": "GTiff",
"height": mosaic.shape[1], # New height of the merged mosaic
"width": mosaic.shape[2], # New width of the merged mosaic
"transform": out_transform, # New georeferencing transform
"crs": src_files_to_mosaic[0].crs, # Coordinate Reference System (assumed consistent)
"nodata": src_files_to_mosaic[0].nodata # NoData value
})
# Close all source files to release resources
for src in src_files_to_mosaic:
src.close()
# Save the merged mosaic to a new GeoTIFF file
try:
with rasterio.open(output_filepath, "w", **out_meta) as dest:
dest.write(mosaic)
print(f"Successfully merged tiles to '{output_filepath}'.")
return output_filepath
except Exception as e:
print(f"Error saving merged GeoTIFF: {e}")
return None
def plot_elevation_tile(filepath):
"""
指定されたGeoTIFF標高ファイルをMatplotlibで表示します。
Args:
filepath (str): 表示するGeoTIFFファイルのパス。
"""
try:
with rasterio.open(filepath) as src:
# 標高データを読み込む（通常はバンド1）
elevation_data = src.read(1)
# NoData値が存在する場合、それをNaN（非数）に変換してプロット時に無視されるようにする
if src.nodata is not None:
elevation_data[elevation_data == src.nodata] = np.nan
# プロットのセットアップ
plt.figure(figsize=(10, 8)) # 図のサイズを設定
# 画像として標高データを表示
# extent: データの地理的な範囲を設定し、軸のラベルを地理座標にする
# cmap: カラーマップ（標高データに適した'viridis'を使用）
# origin: 画像の原点（GeoTIFFは通常左上が原点）
# vmin, vmax: カラーマップの最小値と最大値を設定し、色の範囲を固定する
# ここでは、データの2パーセンタイルと98パーセンタイルを使用
# これにより、外れ値がカラーマップの表示範囲を歪めるのを防ぐ
valid_elevation_data = elevation_data[~np.isnan(elevation_data)]
if valid_elevation_data.size > 0:
vmin = np.percentile(valid_elevation_data, 2)
vmax = np.percentile(valid_elevation_data, 98)
else:
vmin = np.min(elevation_data[~np.isnan(elevation_data)]) if not np.all(np.isnan(elevation_data)) else 0
vmax = np.max(elevation_data[~np.isnan(elevation_data)]) if not np.all(np.isnan(elevation_data)) else 100
# 画像の表示範囲を設定 (左経度, 右経度, 下緯度, 上緯度)
# rasterioのboundsは (left, bottom, right, top) の順
bounds = src.bounds
extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]
im = plt.imshow(
elevation_data,
cmap='viridis', # 標高データに適したカラーマップ
origin='upper', # 画像の原点を左上に設定
extent=extent, # 地理的な範囲を設定
vmin=vmin, # カラーマップの最小値
vmax=vmax # カラーマップの最大値
)
plt.title(f'Elevation Data: {os.path.basename(filepath)}')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.colorbar(im, label='Elevation (meters)') # カラーバーを追加
plt.grid(True, linestyle=':', alpha=0.7) # グリッドを追加
plt.tight_layout() # レイアウトを調整
plt.show()
except rasterio.errors.RasterioIOError as e:
print(f"エラー: GeoTIFFファイル '{filepath}' を開けませんでした: {e}")
except Exception as e:
print(f"予期せぬエラーが発生しました: {e}")
# --- メイン実行ブロック ---
if __name__ == "__main__":
# Define the directory where GeoTIFF files are extracted
extracted_data_dir = '/content/extracted_elevation_data/N030E130_N035E135'
# Define the output path for the merged GeoTIFF file
output_merged_filepath = os.path.join(extracted_data_dir, 'merged_5x5_elevation.tif')
print("Starting GeoTIFF tile merging process...")
merged_file_path = merge_elevation_tiles(extracted_data_dir, output_merged_filepath)
if merged_file_path:
print(f"\nSuccessfully merged all tiles. Now analyzing statistics for: {merged_file_path}")
try:
with rasterio.open(merged_file_path) as src:
# Read the entire merged data
merged_elevation_data = src.read(1)
# Convert NoData values to NaN
if src.nodata is not None:
merged_elevation_data[merged_elevation_data == src.nodata] = np.nan
# Get valid (non-NaN) elevation data for statistics
valid_merged_elevation_data = merged_elevation_data[~np.isnan(merged_elevation_data)]
if valid_merged_elevation_data.size > 0:
print(f"Merged data original pixel size: {src.height}x{src.width} pixels")
print(f"Merged data geographical bounds (Lon/Lat): Left={src.bounds.left:.2f}, Bottom={src.bounds.bottom:.2f}, Right={src.bounds.right:.2f}, Top={src.bounds.top:.2f}")
print(f"Merged data elevation range: Min={np.min(valid_merged_elevation_data):.2f}m, Max={np.max(valid_merged_elevation_data):.2f}m")
print(f"Merged data mean elevation: {np.mean(valid_merged_elevation_data):.2f}m")
print(f"Merged data standard deviation: {np.std(valid_merged_elevation_data):.2f}m")
print(f"Merged data 2nd percentile: {np.percentile(valid_merged_elevation_data, 2):.2f}m")
print(f"Merged data 98th percentile: {np.percentile(valid_merged_elevation_data, 98):.2f}m")
else:
print("No valid elevation data found in the merged file for statistics.")
except Exception as e:
print(f"An error occurred during statistics analysis of the merged file: {e}")
else:
print("GeoTIFF file merging failed, cannot proceed with statistics analysis.")
‘specific_tile_filenames’にListStyleでtifFileのnameをinputして、for_roopでin ordetにreadしてMargerProcessします。完成したMosicImageをMatplotlib（先ほどのCode）でPlotShowすると、このようにSubjectRegionのContourがPlotされました。
Starting GeoTIFF tile merging process...Merging 17 tiles...Successfully merged tiles to '/content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif'.Successfully merged all tiles. Now analyzing statistics for: /content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tifMerged data original pixel size: 18000x18000 pixelsMerged data geographical bounds (Lon/Lat): Left=130.00, Bottom=30.00, Right=135.00, Top=35.00Merged data elevation range: Min=-76.00m, Max=1970.00mMerged data mean elevation: 101.23mMerged data standard deviation: 222.09mMerged data 2nd percentile: 0.00mMerged data 98th percentile: 845.00m
MinimumElavationが-76mというViewが気になりますが、MaximumElavationが1970mというのは西日本最高峰の石鎚山が1982mであることを踏まえると、調和的かなぁと感じています。
ということで、以前、国土地理院Tileのときに実施した代表的なお山の山頂標高がMargeTileでどのような数値になっているかを確認してみます。
以下はGemini先生に作っていただいたCodeでございます。
### 20250728_Latest Version
import rasterio
import numpy as np
import os
def get_elevation_at_peaks(merged_elevation_filepath, mountain_peaks):
"""
結合された標高ファイルから、指定された山々の山頂の標高値を取得し表示します。
Args:
merged_elevation_filepath (str): 結合されたGeoTIFF標高ファイルのパス。
mountain_peaks (dict): 山のデータを含む辞書。各エントリは辞書で、
"lon", "lat", "actual_m" を含む。
"""
print(f"\nAnalyzing elevation at specified mountain peaks using '{merged_elevation_filepath}'...")
if not os.path.exists(merged_elevation_filepath):
print(f"エラー: 結合された標高ファイル '{merged_elevation_filepath}' が見つかりません。")
print("先に結合処理を実行してファイルが生成されていることを確認してください。")
return
try:
with rasterio.open(merged_elevation_filepath) as src:
# 標高データを読み込む（通常はバンド1）
elevation_data = src.read(1)
# NoData値が存在する場合、それをNaN（非数）に変換
if src.nodata is not None:
elevation_data[elevation_data == src.nodata] = np.nan
print(f"Merged data bounds: Left={src.bounds.left:.4f}, Bottom={src.bounds.bottom:.4f}, Right={src.bounds.right:.4f}, Top={src.bounds.top:.4f}")
results = []
for peak_name, peak_info in mountain_peaks.items():
lon = peak_info["lon"]
lat = peak_info["lat"]
actual_m = peak_info["actual_m"]
# 指定された緯度経度座標がデータの範囲内にあるか確認
if not (src.bounds.left <= lon <= src.bounds.right and
src.bounds.bottom <= lat <= src.bounds.top):
results.append({
"name": peak_name,
"lon": lon,
"lat": lat,
"actual_m": actual_m,
"sampled_m": "Out of bounds",
"diff": "N/A"
})
continue # 範囲外の場合は次の山へ
# 緯度経度座標からピクセル行/列に変換
# src.index(lon, lat) は (row, col) を返す
row, col = src.index(lon, lat)
# ピクセル座標が有効な範囲内か確認
if 0 <= row < src.height and 0 <= col < src.width:
sampled_elevation = elevation_data[row, col]
if np.isnan(sampled_elevation):
sampled_elevation_str = "NoData (or NaN)"
diff_str = "N/A"
else:
sampled_elevation_str = f"{sampled_elevation:.2f}m"
diff = sampled_elevation - actual_m
diff_str = f"{diff:.2f}m"
results.append({
"name": peak_name,
"lon": lon,
"lat": lat,
"actual_m": actual_m,
"sampled_m": sampled_elevation_str,
"diff": diff_str
})
else:
results.append({
"name": peak_name,
"lon": lon,
"lat": lat,
"actual_m": actual_m,
"sampled_m": "Pixel out of grid",
"diff": "N/A"
})
# 結果の表示
print("\n--- Mountain Peak Elevation Analysis Results ---")
print(f"{'Mountain Name':<25} | {'Longitude':<12} | {'Latitude':<12} | {'Actual (m)':<12} | {'Sampled (m)':<15} | {'Difference (m)':<15}")
print("-" * 100)
for res in results:
print(f"{res['name']:<25} | {res['lon']:<12.6f} | {res['lat']:<12.6f} | {res['actual_m']:<12} | {res['sampled_m']:<15} | {res['diff']:<15}")
print("-" * 100)
except rasterio.errors.RasterioIOError as e:
print(f"エラー: GeoTIFFファイル '{merged_elevation_filepath}' を開けませんでした: {e}")
except Exception as e:
print(f"予期せぬエラーが発生しました: {e}")
# --- メイン実行ブロック ---
if __name__ == "__main__":
# Define the directory where GeoTIFF files are extracted
extracted_data_dir = '/content/extracted_elevation_data/N030E130_N035E135'
# Define the output path for the merged GeoTIFF file
output_merged_filepath = os.path.join(extracted_data_dir, 'merged_5x5_elevation.tif')
print("Starting GeoTIFF tile merging process...")
# 結合処理を実行（ファイルが存在しない場合のみ）
if not os.path.exists(output_merged_filepath):
merged_file_path = merge_elevation_tiles(extracted_data_dir, output_merged_filepath)
else:
print(f"Merged file '{output_merged_filepath}' already exists. Skipping merge process.")
merged_file_path = output_merged_filepath # 既存のファイルパスを使用
if merged_file_path:
# 標高データの統計情報を表示（既存のコードからコピー）
print(f"\nSuccessfully merged/found tiles. Now analyzing statistics for: {merged_file_path}")
try:
with rasterio.open(merged_file_path) as src:
merged_elevation_data = src.read(1)
if src.nodata is not None:
merged_elevation_data[merged_elevation_data == src.nodata] = np.nan
valid_merged_elevation_data = merged_elevation_data[~np.isnan(merged_elevation_data)]
if valid_merged_elevation_data.size > 0:
print(f"Merged data original pixel size: {src.height}x{src.width} pixels")
print(f"Merged data geographical bounds (Lon/Lat): Left={src.bounds.left:.2f}, Bottom={src.bounds.bottom:.2f}, Right={src.bounds.right:.2f}, Top={src.bounds.top:.2f}")
print(f"Merged data elevation range: Min={np.min(valid_merged_elevation_data):.2f}m, Max={np.max(valid_merged_elevation_data):.2f}m")
print(f"Merged data mean elevation: {np.mean(valid_merged_elevation_data):.2f}m")
print(f"Merged data standard deviation: {np.std(valid_merged_elevation_data):.2f}m")
print(f"Merged data 2nd percentile: {np.percentile(valid_merged_elevation_data, 2):.2f}m")
print(f"Merged data 98th percentile: {np.percentile(valid_merged_elevation_data, 98):.2f}m")
else:
print("No valid elevation data found in the merged file for statistics.")
except Exception as e:
print(f"An error occurred during statistics analysis of the merged file: {e}")
# --- ここから山頂標高の確認・表示 ---
# Mountain data (updated with your requested changes)
mountain_peaks = {
"Mt. Ishizuchi (石鎚山)": {"lon": 133.11423824903136, "lat": 33.76840969420332, "actual_m": 1982},
"Mt. Tsurugi (剣山)": {"lon": 134.0940772370927, "lat": 33.85505389948576, "actual_m": 1955},
"Mt. Sanrei (三嶺)": {"lon": 133.9877611877695, "lat": 33.840358759990025, "actual_m": 1893},
"Mt. Yufu (由布岳)": {"lon": 131.3903, "lat": 33.2822, "actual_m": 1583},
"Mt. Karakuni (韓国岳)": {"lon": 130.86118851784443, "lat": 31.942399822041573, "actual_m": 1700},
"Mt. Takachiho (高千穂峰)": {"lon": 130.9189, "lat": 31.8860, "actual_m": 1574},
"Mt. Kaimondake (開聞岳)": {"lon": 130.5283, "lat": 31.1800, "actual_m": 924},
"Mt. Miyanoura (宮之浦岳)": {"lon": 130.5042, "lat": 30.3359, "actual_m": 1936},
"Mt. Daisen (大山)": {"lon": 133.5400, "lat": 35.3711, "actual_m": 1729},
"Mt. Aso (阿蘇山)": {"lon": 131.1000, "lat": 32.8800, "actual_m": 1592},
"Mt. Kuju (久住山)": {"lon": 131.2409422180058, "lat": 33.08243796768661, "actual_m": 1787},
"Mt. Hiko (英彦山)": {"lon": 130.9260117108225, "lat": 33.477174995714954, "actual_m": 1199}
}
get_elevation_at_peaks(merged_file_path, mountain_peaks)
else:
print("GeoTIFF file merging failed or file not found, cannot proceed with analysis.")
>>>>Output<<
Starting GeoTIFF tile merging process...
Merged file '/content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif' already exists. Skipping merge process.
Successfully merged/found tiles. Now analyzing statistics for: /content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif
Merged data original pixel size: 18000x18000 pixels
Merged data geographical bounds (Lon/Lat): Left=130.00, Bottom=30.00, Right=135.00, Top=35.00
Merged data elevation range: Min=-76.00m, Max=1970.00m
Merged data mean elevation: 101.23m
Merged data standard deviation: 222.09m
Merged data 2nd percentile: 0.00m
Merged data 98th percentile: 845.00m
Analyzing elevation at specified mountain peaks using '/content/extracted_elevation_data/N030E130_N035E135/merged_5x5_elevation.tif'...
Merged data bounds: Left=130.0000, Bottom=30.0000, Right=135.0000, Top=35.0000
--- Mountain Peak Elevation Analysis Results ---
Mountain Name | Longitude | Latitude | Actual (m) | Sampled (m) | Difference (m)
----------------------------------------------------------------------------------------------------
Mt. Ishizuchi (石鎚山) | 133.114238 | 33.768410 | 1982 | 1928.00m | -54.00m
Mt. Tsurugi (剣山) | 134.094077 | 33.855054 | 1955 | 1929.00m | -26.00m
Mt. Sanrei (三嶺) | 133.987761 | 33.840359 | 1893 | 1834.00m | -59.00m
Mt. Yufu (由布岳) | 131.390300 | 33.282200 | 1583 | 1557.00m | -26.00m
Mt. Karakuni (韓国岳) | 130.861189 | 31.942400 | 1700 | 1589.00m | -111.00m
Mt. Takachiho (高千穂峰) | 130.918900 | 31.886000 | 1574 | 1568.00m | -6.00m
Mt. Kaimondake (開聞岳) | 130.528300 | 31.180000 | 924 | 919.00m | -5.00m
Mt. Miyanoura (宮之浦岳) | 130.504200 | 30.335900 | 1936 | 1927.00m | -9.00m
Mt. Daisen (大山) | 133.540000 | 35.371100 | 1729 | Out of bounds | N/A
Mt. Aso (阿蘇山) | 131.100000 | 32.880000 | 1592 | 1419.00m | -173.00m
Mt. Kuju (久住山) | 131.240942 | 33.082438 | 1787 | 1772.00m | -15.00m
Mt. Hiko (英彦山) | 130.926012 | 33.477175 | 1199 | 1157.00m | -42.00m
----------------------------------------------------------------------------------------------------
阿蘇山で誤差が173m、韓国岳で111mとなっています。Errorがまあまあ大きいので、それぞれの山頂の緯度経度をFoliumを使って地図上に落としてみます。
In relatively larger scale, The summits of mailly mountain locate almost exact location.
But 詳細なlocationはAlosの解像度が30mMeshであることに起因して、国土地理院が示している平面的な山頂の位置とは少しずれるケースがあることがわかりました。
コメント
コメントを投稿