登山好きプログラミング素人が衛星データで登りたい山を探してみた！〜②衛星データで周囲で一番高い山を探してみる編〜

はじめに

・登山好きプログラミング初心者が衛星データを利用して山探しする企画です。
・こう書いた方がより良いなどありましたら、勉強になりますのでぜひコメントください。
・本企画は以下の４部構成を予定しており、今回は第２弾です。
①日本の山リスト作成編
②衛星データで周囲で一番高い山を探してみる編
③衛星データで山周辺の地域の晴天率を調べてみる編
④衛星データで眺望の良い、山頂がひらけている山を探してみる編
前回記事(①日本の山リスト作成編)で作成した日本の山リストに対して、標高の衛星データを使って自分たちの登りたい山を絞っていきます。

今回は以下の二つの条件で、登りたい山を絞っていきます。
・標高2000m以上
・周辺で最も高い山

標高２０００m以上の山を絞る

標高2000m以上の山を絞ります。
衛星データを用いて絞ることもできますが、前回記事で作成したcsvに標高の値があるため、今回はそちらを使って絞ります。

まずは前回記事で作成した日本の山リストのcsvを読み込みます。
その後、標高の列の値が2000以上のものを抽出します。
抽出後は行数の値が飛び飛びになってしまっているため、reset_index()で番号を振り直します。
最後にcsvで保存します。

df = pd.read_csv("/content/drive/My Drive/mountain_list_fin.csv", index_col=0)

# 標高2000m以上を抽出
df2000 = df[df['標高'] >= 2000]
df2000.reset_index(drop=True, inplace=True)
df2000.to_csv("/content/drive/My Drive/mountain_list_2000.csv")

この結果を前回と同じようにGoogleマップに表示してみます。
前回の結果から山が絞られた様子が分かります。
驚いたことに2000m以上の山は西日本にはないんですね…

周辺で最も高い山の絞り方

山頂に着いたときに、近くに自分がいる位置より高い山があるのはなんか悔しいですよね。
そのために自分が周辺で最も高い山を衛星データを用いて絞っていきます。
以下の図のように山頂を中心に大小二つの正方形の領域で標高データを取得し、その最大値を比較します。最大値が同じだった場合には青い領域でも、中心の山頂が最も標高が高かったことになるため、この山は周囲で最も高い山であると言えます。逆に最大値が異なっていた場合は、青い領域に中心の山頂以上に高い標高の場所があると言えるため、この山は周囲で最も高い山ではないことがわかります。
この方法を用いて、周囲で最も高い山を絞ります。

標高データの紹介

今回使用する衛星データは「ALOS DSM Global 30m」です。
https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_AW3D30_V3_2
こちらのデータのDSMというバンドで標高データを取得することができます。各衛星には複数のバンドがあり、バンド名は上記のGEEのデータカタログにて調査が可能です。
まず試しに富士山を中心として画像を取得し表示してみます。

GEEのデータを使用するため、GEEの認証とモジュールのインポートを行います。
GEE認証を実行すると認証に必要なURLが表示されるので、URLへアクセスしてGoogleアカウントを指定し、表示された認証コードをGoogle Colabのボックスへコピペしてください。

!pip install rasterio #初めて実行するときのみ
import ee
import numpy as np
import matplotlib.pyplot as plt
import rasterio

try:
  ee.Initialize()
except Exception as e:
  ee.Authenticate()
  ee.Initialize()

衛星データをロードする関数を定義します。
衛星名(snippet)、場所(geometry)、バンド(band)を指定して、衛星データをロードします。通常はどの期間のデータをロードするのかも指定しますが、今回使用する標高データはデータ数が一つのため、不要です。

# GEEのデータロード
def load_data(snippet, geometry, band):
    # パラメータの条件にしたがってデータを抽出
    dataset = ee.ImageCollection(snippet).filter(
    ee.Filter.geometry(geometry)).select(band)
    # リスト型へ変換
    data_list = dataset.toList(dataset.size().getInfo())
    # 対象データ数とデータリストを出力
    return dataset.size().getInfo(), data_list

次に衛星名(snippet)、バンド(band)、関心領域(aoi)を指定します。
場所には富士山を中心として、緯度経度0.2[°]の四角形を指定しています。

## パラメーターの指定
# 衛星を指定
snippet = 'JAXA/ALOS/AW3D30/V3_2'
# バンド名を指定
band = 'DSM'

# 富士山周辺を指定
lon = 138.727363
lat = 35.360626
aoi = ee.Geometry.Rectangle([lon - 0.1, lat - 0.1, lon + 0.1, lat + 0.1])

データリストを取得し、データ数を確認します。
今回はデータ数は”１”と表示されます。

##　データ数の確認
num_data, data_list = load_data(snippet=snippet, geometry=aoi, band=band)
print('Datasets; ', num_data)

指定したデータをGoogle Drive内にTIFF(.tif)ファイルで保存します。
下記のコードを実行すると、My Drive > fijisan_test > にfuji_muntain.tifというファイルが保存されます。

image = ee.Image(data_list.get(0))

# Gdriveへ保存
task = ee.batch.Export.image.toDrive(**{
    'image': image,# 対象データの指定
    'description': 'fuji_mountain',# ファイル名の指定
    'folder':'fijisan_test',# Google Drive(MyDrive)のフォルダ名
    'scale': 30,# 解像度の指定
    'region': aoi.getInfo()['coordinates']# 上記で指定した対象エリア
})

# 処理の実行
task.start()

データの取得には時間がかかるので、下記コードを実行して、Falseと表示されるまで待ちます。

# Falseと表示されれば処理終了
task.active()

無事データが取得されたら、最大標高を取得し、画像を表示します。
No such file or directoryというエラーが出た場合には、まだデータのダウンロード中か、Google Driveのマウントができていないかと思われます。確認してみてください。

# データの読み込み
with rasterio.open('/content/drive/My Drive/fijisan_test/fuji_mountain.tif') as src:
    arr = src.read()

print("最大標高: ", arr.max())

# 可視化
plt.imshow(arr[0])

以下のように表示されました。
富士山の標高3776mと比べて多少誤差があるものの、無事富士山山頂を中心とする画像の取得ができました。

reducerを用いて画像の最大値を取得する

富士山の例のように、全ての山に対して画像ファイルをダウンロードして、それを読み取る方法で画像の最大値を取得することができます。しかし今回のように、複数の画像を調査したい、代表値のみ取得できれば良い、といった際には”reducer”という機能が便利です。
reducerを使えば、画像をローカルにダウンロードせずに、画像の最大値を取得することができます。reducerを使って富士山の最大標高を取得してみましょう。

先程取得したimageに対してreducerを適用します。
reduceRegion()に、reducer、geometry、scaleを与えることで、画像の代表値を取得することができます。reducerにはmax()、min()、meam()、median()などがあります。

max = image.reduceRegion(reducer=ee.Reducer.max(), geometry=aoi, scale=30)
print(max.values().getInfo()[0])

実行結果は以下のようになりました。
上と同じように”3771”という富士山の最高標高が取得できています。
reducerを用いることで簡単に最大値を取得することができました。
一枚の画像ではあまり効果を感じないかもしれませんが、たくさんの画像に対して処理を行いたい場合にはとても便利な機能です。

周囲で最も高い山を探す

いよいよ衛星データを使って登りたい山を絞っていきます。
まずはGEEにログインします。

import ee
import pandas as pd
import csv

try:
 ee.Initialize()
except Exception as e:
 ee.Authenticate()
 ee.Initialize()

上と同じように衛星データをロードする関数を定義します。

def load_data(snippet, geometry, band):
   # パラメータの条件にしたがってデータを抽出
   dataset = ee.ImageCollection(snippet).filter(
   ee.Filter.geometry(geometry)).select(band)
   # リスト型へ変換
   data_list = dataset.toList(dataset.size().getInfo())
   # データリストを出力
   return data_list

衛星名とバンドを指定します。
ここまでは上と同じです。

## パラメーターの指定
# 衛星を指定
snippet = 'JAXA/ALOS/AW3D30/V3_2'
# バンド名を指定
band = 'DSM'

画像から関心領域の最大値を取得する関数を定義します。
reducerを用いて画像の最大値を取得し、返しています。

def get_max(aoi):
 data_list = load_data(snippet=snippet, geometry=aoi, band=band)
 image = ee.Image(data_list.get(0))
 max = image.reduceRegion(reducer=ee.Reducer.max(), geometry=aoi, scale=30)
 return max.values().getInfo()[0]

先ほど作成した標高2000m以上の山リストを基に周辺で最も高い山を絞っていきます。
リストにある山に対して、ループ処理で以下の①〜④を適用していきます。
①リストから山頂の緯度経度を取得します。
②緯度経度0.02[°]の正方形の標高データから最大値を取得
③緯度経度0.2[°]の正方形の標高データから最大値を取得
④②と③の最大値を比較し、異なっていればリストから削除する

df_top = df = pd.read_csv("/content/drive/My Drive/mountain_list_2000.csv", index_col=0)
for i in range(len(df_top)):
 print(i)
 lat = df_top.at[i, "緯度"]
 lon = df_top.at[i, "経度"]

 aoi_inner = ee.Geometry.Rectangle([lon - 0.01, lat - 0.01, lon + 0.01, lat + 0.01])
 max1 = get_max(aoi_inner)

 aoi_outer = ee.Geometry.Rectangle([lon - 0.1, lat - 0.1, lon + 0.1, lat + 0.1])
 max2 = get_max(aoi_outer)

 if max1 != max2:
   df_top = df_top.drop(i, axis=0)

最後に結果を保存します。

df_top.reset_index(drop=True, inplace=True)
df_top.to_csv("/content/drive/My Drive/mountain_list_top.csv")

結果をGoogleマップに表示します。
山は32個に絞られました。

高い山の連なる北アルプス周辺を見てみます。
北アルプス最高峰は奥穂高岳なので、その周辺の常念岳、焼岳、槍ヶ岳などは周囲で一番高い山ではないと判断されていることがわかります。
しかし涸沢岳のようにかなり近くに2000m以上の山がある場合には、内側の小さな長方形内に二つの山が入ってしまい、どちらも残ってしまっています。厳密に絞りたいのであればここは要改善です。
しかし今回は概ね周囲で一番高い山を絞ることができているため、このまま次のステップに進みたいと思います。

まとめ

今回はpandasを用いて2000m以上の山、標高の衛星データを用いて周囲で最も高い山を絞りました。一応良好な結果を得られたのではないでしょうか。
次回からは気象の衛星データを用いて晴天率の高い山、植生指数という衛星データを用いて山頂地点での見晴らしが良い山を探していきます。
③衛星データで山周辺の地域の晴天率を調べてみる編

参考文献

※１　 GEEデータカタログ：ALOS DSM Global 30m
https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_AW3D30_V3_2
※２　ImageCollection Reductions
https://developers.google.com/earth-engine/guides/reducers_image_collection