前回：https://shikuuk.blogspot.com/2025/07/jaxaalospython.html
第２部はさきほどDownloadした「N030E130_N035E135.zip」に含まれている５×５のTileを１枚の大きなModelに結合する処理（モザイク処理と呼ぶそうです）のCodeをGemini先生にお願いしました。
※JAXAさま AW3D30 DSM data map Download Page参照
第２部はさきほど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を作成する。
続く。
コメント
コメントを投稿