🛰

衛星データを使って関東近辺の地形を3D表示させる

2023/04/08に公開

Python

衛星データ

tech

無料で取得できるデータと Python を使って、関東近辺の地形を 3D 表示させます。

ここでは、NASA の SRTM (Shuttle Radar Topography Mission) データを使用し、Python で 3D モデルを作成します。

実行環境

Ubuntu 22.04
Python 3.10.4

データの取得

SRTM データは、NASA が提供する無料の地形データです。

USGS EarthExplorer から SRTM データを取得する

USGS EarthExplorer にアクセス

以下のリンクから USGS EarthExplorer にアクセスします。

https://earthexplorer.usgs.gov/
アカウント作成

データをダウンロードするには、USGS EarthExplorer のアカウントが必要です。
登録ページ（https://ers.cr.usgs.gov/register/ ）でアカウントを作成してください。
地域とデータセットの選択

EarthExplorer のメインページで、検索ボックスに緯度経度や地名を入力して検索し、取得したい地域を選択することができます。
今回はマップをクリックしてエリアを指定し、データを検索します。

エリアを指定した後は「Data Sets」ボタンをクリックし、「Digital Elevation」カテゴリを展開し、「SRTM」を選択します。ここで SRTM 1 Arc-Second Global を選択します。
データのダウンロード

「Results」タブをクリックし、データセットのリストから適切なものを見つけます。

今回は比較のために以下３つのデータをダウンロードします。
- 東京を含むデータ
- 関東南側の海のデータ
- 富士山を含むデータ
ダウンロードアイコンをクリックして、データファイル（GeoTIFF ファイル）をダウンロードしてください。

画像の赤い枠が検索のために手動でプロットしたものです。茶色や紫色がダウンロードできるデータの範囲です。

Python ライブラリのインストール

以下の Python ライブラリを使います。

rasterio: 地形データの読み込みに使用
numpy: データ処理に使用
matplotlib: 3D プロットに利用

インストールは以下のように行います。

pip install rasterio numpy matplotlib

実装

ダウンロードした SRTM データから 3D モデルを作成し表示する関数を実装します。

なお、Jupyter を利用することを想定しています。

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# SRTMデータを読み込む関数
def read_elevation_data(file_path):
    with rasterio.open(file_path) as src:
        elevation_data = src.read(1)
        elevation_data = elevation_data.astype(np.float32)
        elevation_data[elevation_data < 0] = np.nan
    return elevation_data

# SRTMデータを3Dプロットする関数
def plot_elevation_data(elevation_data_list):
    fig = plt.figure(figsize=(20, 8))

    for i, elevation_data in enumerate(elevation_data_list):
        ax = fig.add_subplot(1, len(elevation_data_list), i+1, projection="3d")
        x = np.arange(0, elevation_data.shape[1], 1)
        y = np.arange(0, elevation_data.shape[0], 1)
        x, y = np.meshgrid(x, y)
        ax.plot_surface(x, y, elevation_data, cmap="terrain")

        # 東西南北の軸名を表示
        ax.set_xlabel('west --- east')
        ax.set_ylabel('south --- north')
        ax.set_zlabel('hight [m]')

        # 高さ軸の最大値を3000に固定
        ax.set_zlim(0, 3000)

    plt.show()

描画

SRTM データは python スクリプトと同階層に保存しておいてください。

では定義した関数に、ダウンロードしたファイルを指定して、3D 表示させてみます。

file_paths = ["./n35_e139_1arc_v3.tif",
              "./n34_e139_1arc_v3.tif",
              "./n35_e138_1arc_v3.tif"]

elevation_data_list = [read_elevation_data(file_path) for file_path in file_paths]
plot_elevation_data(elevation_data_list)

結果は以下のようになります。

一番左は、東京を含む関東平野の結果です。東京や神奈川は青色でほぼ平坦になっています。西の方だけ隆起しており、関東の西の山があることがわかります。

真ん中は、関東南側の近海の結果です。ポツポツと島の隆起があることがわかりますが、それ以外は海なので青く真っ平らになっています。

一番右は、富士山近辺です。ほとんど山の範囲のデータなので激しく隆起していることがわかります。

まとめ

このように衛星データと Python を使って、地球の高さ情報を可視化することができます。

今回使った SRTM データはデジタル標高モデル（DEM: Digital Elevation Model）と呼ばれるデータの一種です。DEM データは地球上の地形や標高を表現するデータで、標高の数値が格子状に配置されたデータセットです。

今回は無料のデータを使ったため、解像度が粗く、詳細な情報は得られませんでしたが、解像度が細かい衛星データ（おそらく有償）を使うことでより詳細な高さ情報を得ることができます。

この DEM データは、以下のような様々な分野で利用されます。

地形解析: 地形の特徴や地形変化を把握し、傾斜、起伏などの地形要素を計算するために使用する。
洪水リスク評価: 洪水が起こりやすい場所や浸水範囲をシミュレーションし、防災対策や避難計画を立てるために使用する。
農業・土地利用: 土壌や水資源の分布を解析し、農地や森林などの土地利用計画や土地管理を行う。
交通・インフラ計画: 道路や鉄道、ダムなどのインフラ建設計画において、最適なルートや位置を決定するために使用される。

また、今回は matplotlib を用いて DEM データを可視化しましたが、もっとインタラクティブに可視化したい場合は、plotly や PyVista といったライブラリもあります。ただし、サイズの大きい衛星データを扱うにはある程度のマシンスペックが必要になるのでご注意ください。

https://plotly.com/python/

実行環境

データの取得

USGS EarthExplorer から SRTM データを取得する

Python ライブラリのインストール

実装

描画

まとめ

Discussion