GIS、リモセンで地理問題を解析しよう！

こんにちわ、今日は記事にしようとおもっていたものの、３か月近く放置していたGISとGOOGLE EARTH ENGINEの衛星データを使って特定地域の水の流域面積変化を測定しようという内容です。
これは名古屋大学の衛星データの授業の一環としてのグループワークで内容の一部です。
内容はメコン川の流域の人々が上流のダムの建設によってどのような影響受けたかという壮大なテーマではあったのですが、私はひたすら流域面積の測定に特化しました。その防備録です。

準備１　メコン川の場所の特定と測定範囲を絞る

メコン川の流域ですが以下からとてきました。
https://data.opendevelopmentmekong.net/dataset/greater-mekong-subregion-rivers
私はGOOGLE COLABを使いましたが直接おいてもいいですし、いったんGDRIVEの保存してそこのパスから読んであげてもいいかもしれません

import folium
import geopandas as gpd
import json

# 1) Shapefile 読み込み
mekong_shp_path = input_path + "mekong4/Mekong_River_Full.shp"
mekong_gdf = gpd.read_file(mekong_shp_path)
wrs2_shp_path = input_path + "WRS2/WRS2_descending.shp"
wrs2_gdf = gpd.read_file(wrs2_shp_path)

# 2) 表示用マップの作成（中心座標はメコン川のジオメトリから自動取得）
mekong_centroid = mekong_gdf.geometry.unary_union.centroid
map_center = [mekong_centroid.y, mekong_centroid.x]
m = folium.Map(location=map_center, zoom_start=6, tiles="CartoDB Positron")

# 3) GeoDataFrame → GeoJSON 変換
mekong_geojson = json.loads(mekong_gdf.to_json())
wrs2_geojson   = json.loads(wrs2_gdf.to_json())

# 4) Folium へ GeoJSON レイヤーとして追加
folium.GeoJson(
    mekong_geojson,
    name="Mekong River",
    style_function=lambda feat: {
        'color': 'blue',
        'weight': 2,
        'fillOpacity': 0
    }
).add_to(m)
folium.GeoJson(
    wrs2_geojson,
    name="WRS-2 Descending",
    style_function=lambda feat: {
        'color': 'red',
        'weight': 1,
        'fillOpacity': 0
    }
).add_to(m)
# 5) レイヤーコントロールを追加して ON/OFF 切り替え可能にする
folium.LayerControl().add_to(m)
# 6) 地図を表示（Jupyter Notebook 等では m を最後に置くか display(m)）
m

上記を実行すると最後で地図が表示されます。またこのSHPファイルのベクターデータも表示されます。
でも結構いい加減なデータです。川と大幅にずれている。これはあるあるなのですが、研究対象とする地域でこういったデータがきちんと整備されていないことはよくあることです。今回はそれでもこういったシェープファイルが見つかっただけ幸運とは言えます。

ちなみに上の図は、WRSというグリッドを非表示にしてます。表示すると以下のような見栄えに。

この赤い四角はWRS（Worldwide Reference System）というもので、Landsat 衛星の撮像範囲を「Path（軌道）」と「Row（列）」の格子状グリッドで管理する参照システムです。
Landsat衛星3以前ではWRS-1 をLandsat衛星4以降は WRS-2 を使ってます。
今回はLandsat衛星５の画像を使うためWRS-2のグリッドを導入してます。
WRS-2は地球全体を 233 の Path（軌道） × 248 の Row（列）に分けています。これを入れている理由はあとで対象とする地域をGEOプロセスでいうところで衛星画像と先ほどのシェープファイルの場所を重ね合わせるために後でＣＬＩＰ（切り抜き）をするために使ってます。

準備2 対象地域を絞る

これから以下の操作をします。
** WRS2とメコン川のシェープファイルの座標系を合わせます。
** メコン川の小域を対象地域（WRS2で指定）で切り取ります（CLIP）
** メコン川の川の場所の精度が少し心配なのでシェープファイルからバッファ化します
** さらに生成されたバッファの中で必要なものだけ取り出します。

●WRS2とメコン川のシェープファイルの座標系を合わせておきます。
（地球を観測して地図にするにはいろいろな座標系という基準がありそれを合わせないと重なりません,この例ではEPSG:4326を利用してます。EPSG:4326 は、地球上の位置を「緯度／経度（latitude／longitude）」で表す座標系（CRS）で、基準とする測地系（ジオイドや楕円体）は WGS84（World Geodetic System 1984）です）
以下の例では普段使われているGPSと同じ基準での座標系に合わせてます。

# メコン川データの座標系を確認し、必要に応じて変換
if mekong_gdf.crs.to_string() != "EPSG:4326":
    mekong_gdf = mekong_gdf.to_crs(epsg=4326)
# WRS-2データの座標系を確認し、必要に応じて変換
if wrs2_gdf.crs.to_string() != "EPSG:4326":
    wrs2_gdf = wrs2_gdf.to_crs(epsg=4326)

●次にメコン川の小域を対象地域（WRS2で指定）で切り取ります（CLIP）

# 対象とするPATH/ROWを選択
selected_area_paths_row = wrs2_gdf[
    (wrs2_gdf['PATH'] == 130) & (wrs2_gdf['ROW'] == 46)
]
# メコン川のデータを選択した領域でクリップ
mekong_area_clip = gpd.clip(mekong_gdf, selected_area_paths_row)

●次にシェープファイルより5000mでバッファをとります。

# バッファ距離を設定（単位はメートル）
buffer_distance = 5000  # 例: 5000メートル
# 適切な投影法に変換（例: UTMゾーン48N）
mekong_area_clip = mekong_area_clip.to_crs(epsg=32648)
# バッファポリゴンを作成
mekong_buffer = mekong_area_clip.buffer(buffer_distance)
# バッファポリゴンをGeoDataFrameに変換
buffer_gdf = gpd.GeoDataFrame(geometry=mekong_buffer, crs=mekong_area_clip.crs)
# 元の座標系（EPSG:4326）に再投影
buffer_gdf = buffer_gdf.to_crs(epsg=4326)

一旦ここで地図でみてみます。

selected_buffer = buffer_gdf
selected_buffer_geojson = json.loads(selected_buffer.to_json())
# Foliumマップの作成（選択した領域の中心を指定）
map_center = [selected_buffer.geometry.centroid.y.mean(), selected_buffer.geometry.centroid.x.mean()]
m = folium.Map(location=map_center, zoom_start=10, tiles='CartoDB Positron')
# 選択したバッファ領域のみ表示（青色）
folium.GeoJson(
    selected_buffer_geojson,
    name="Selected Buffer Geometry",
    style_function=lambda feature: {'color': 'blue', 'fillOpacity': 0.3}
).add_to(m)
# WRS2 の選択領域（PATH=130, ROW=46）の GeoDataFrame を GeoJSON に変換
selected_area_geojson = json.loads(selected_area_paths_row.to_json())
# WRS2 の選択領域を境界線のみ（赤）で追加
folium.GeoJson(
    selected_area_geojson,
    name="Selected WRS2 Boundary",
    style_function=lambda feature: {'color': 'red', 'fillOpacity': 0}
).add_to(m)
# レイヤーコントロールを追加
folium.LayerControl().add_to(m)
# マップの表示
m

バッファの対象地域配列が最終的にだされます。私のテーマではもっと絞るため、と規定の場所だけにしました。上のコードで

# 例: buffer_gdf の3番目のフィーチャー（インデックス 2）だけを選択
#selected_buffer = buffer_gdf
selected_buffer = buffer_gdf.iloc[[3]]

としてください。複数のバッファのうち対象地域だけにさらに絞ります。

準備3 GOOGLE EARTH　ENGINEに登録してPROJECTを作成

https://earthengine.google.com/
から登録します。非商用であれば無料で使えます。商用では月額課金がされます。通常衛星の生データはそのままでは使えないことが多く、大気補正や雲の補正など様々な補正が必要で衛星画像研磨の技術がないと極めて敷居が高い代物ですが、GOOGLE ERTH　ENGINEはすでにいろいろな指標がそのまま使える形でデータを得られます。例えば光合成指標などいろいろな直接観測値から指標計算済みのデータも場所と年代によっては得られます。
レベル０ 生データ Raw Data
レベル１ 輝度データ+軌道情報のキャリブレーション
レベル２ 正規化輝度+大気補正+バンド間演算
レベル３ 　地球規模幾何補正と複合化
レベル４ 他センサー情報で指標計算
それでは登録をしてPROJECTIDを取っておきます。そしてそのPROJECTIDを使って以下でEARTH ENGINEぼライブラリをPYTHONでインポートします。

import ee  # Earth Engineライブラリ
ee.Authenticate()
ee.Initialize(project='ee-YourProjectID')  # ご自身のプロジェクトIDを入力してください

まず、興味のある年代でほしい衛星の画像がどの程度あるかをチェックしておきましょう

for year in range(1987, 2000):
    start_date = ee.Date.fromYMD(year, 1, 1)
    end_date = ee.Date.fromYMD(year, 12, 31)
    # 指定した期間の画像数をカウント
    count = landsat5.filterDate(start_date, end_date).size().getInfo()
    # 結果を表示
    print(f"{year}: {count}")

1987: 1
1988: 10
1989: 11
1990: 16
1991: 10
1992: 11
1993: 12
1994: 17
1995: 15
1996: 11
1997: 13
1998: 13
1999: 14
とおおむね1988~1999年には１月１回弱くらいで撮影はされているようです。

では先ほど特定した場所の画像を人間の視覚（RGB）で表示してみましょう。
ここでLANDSAT/LT05/C02/T1_TOA　はGOOGLE EARTH ENGINEのカタログIDです
また'1987-01-01', '1999-12-31'と日にちを絞ってます。返される画像はこれらの平均をとったものになります。

まず再度、目的の場所を特定して

# buffer_gdfをGeoJSON形式に変換
buffer_geojson = json.loads(buffer_gdf.to_json())
# GeoJSONからEarth Engineのジオメトリに変換
buffer_geometry = ee.Geometry(buffer_geojson['features'][3]['geometry'])

表示します。

landsat5 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_TOA')
    .filterBounds(buffer_geometry)
    .filterDate('1987-01-01', '1999-12-31')
    .map(lambda img: img.clip(buffer_geometry))
)
# 4. たとえば 1987 年の median composite を作成し、RGB 用に可視化パラメータを設定
year = 1987
start_date = ee.Date.fromYMD(year, 1, 1)
end_date   = ee.Date.fromYMD(year, 12, 31)
composite  = landsat5.filterDate(start_date, end_date).median()
# Landsat 5 でのバンド 3 (Red), 2 (Green), 1 (Blue) を使って表示
vis_params = {
    'bands': ['B3', 'B2', 'B1'],
    'min': 0,      # TOA リフレクタンス値の下限（適宜変更可）
    'max': 0.3,    # TOA リフレクタンス値の上限（適宜変更可）
    'gamma': 1.2   # コントラスト調整（任意）
}
# 5. Folium マップを作成し、カスタム関数で EE 画像レイヤーを追加できるよう登録
def add_ee_layer(self, ee_image_object, vis_params, name):
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
        tiles=map_id_dict['tile_fetcher'].url_format,
        attr="Map Data © Google Earth Engine",
        name=name,
        overlay=True,
        control=True
    ).add_to(self)
folium.Map.add_ee_layer = add_ee_layer
# 6. マップの中心を設定（buffer_geometry の重心を中心にする例）
centroid = buffer_geometry.centroid().coordinates().getInfo()
map_center = [centroid[1], centroid[0]]  # folium は [lat, lon]
m = folium.Map(location=map_center, zoom_start=10, tiles='CartoDB Positron')
# 7. 先ほどの可視化パラメータを使って RGB レイヤーを追加
m.add_ee_layer(composite, vis_params, f'Landsat5_{year}_RGB')
# 8. レイヤーコントロールを表示
folium.LayerControl().add_to(m)
# 9. マップを表示（Jupyter Notebook 等ではこれでレンダリングされます）
m

きれいに映っていると思います。ただし今回の目的は水の領域を知りたいのです。
RGB画像表示が目的ではないので、水を特定するための指標で計算して表示してます。
ここではNDWIという指標を使ってます。
Landsat 5 TM を使う場合、NDWI は緑色バンド（B2）と近赤外（NIR）バンド（B4）の輝度値を使って次の式で算出します：

B2（Green, 緑バンド）：水面は可視光の緑色領域で比較的高い反射を示す
B4（NIR, 近赤外バンド）：水面は近赤外領域ではほぼ吸収してしまい、非常に低い反射を示す
この２つのバンドの差（Green － NIR）が大きいほど、水面の特徴が強調されます。

# Landsat 5のデータを取得
landsat5 = (ee.ImageCollection('LANDSAT/LT05/C02/T1_TOA')
                    .filterBounds(buffer_geometry)
                    .filterDate('1987-01-01', '1999-12-31')
                    .map(lambda image: image.clip(buffer_geometry)))

# まず results_ndwi をタプル形式で再作成
results_ndwi = []
for year in range(1987, 2000):
    start_date = ee.Date.fromYMD(year, 1, 1)
    end_date = ee.Date.fromYMD(year, 12, 31)
    composite = landsat5.filterDate(start_date, end_date).median()
    ndwi = composite.normalizedDifference(['B2', 'B4'])
    ndwi = ndwi.updateMask(ndwi.gte(0))
    results_ndwi.append({str(year): ndwi})

それでは表示します。

# Folium 地図の中心は対象エリアの centroid から計算
centroid = buffer_geometry.centroid().coordinates().getInfo()
map_center = [centroid[1], centroid[0]]  # folium は [lat, lon]
m = folium.Map(location=map_center, zoom_start=10, tiles='CartoDB Positron')
# WRS2 の選択領域を境界線のみ（赤）で追加
folium.GeoJson(
    selected_area_geojson,
    name="Selected WRS2 Boundary",
    style_function=lambda feature: {'color': 'red', 'fillOpacity': 0}
).add_to(m)
# 対象領域（バッファ）の表示（青）
folium.GeoJson(
    selected_buffer_geojson,
    name="Selected Buffer Geometry",
    style_function=lambda feature: {'color': 'yellow', 'fillOpacity': 0}
).add_to(m)
# ===================================================================
# Folium に Earth Engine のレイヤーを追加するためのカスタム関数を定義
def add_ee_layer(self, ee_image_object, vis_params, name):
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
        tiles = map_id_dict['tile_fetcher'].url_format,
        attr = "Map Data &copy; <a href='https://earthengine.google.com/'>Google Earth Engine</a>",
        name = name,
        overlay = True,
        control = True
    ).add_to(self)
# この関数を folium.Map に追加
folium.Map.add_ee_layer = add_ee_layer
# ===================================================================
# NDWI レイヤーを Folium 地図に追加
# すべてのレイヤーで同じ青系統の色を使用（例: "#0000FF"）
blue_color = "#0000FF"
for item in results_ndwi:
    for key, ndwi in item.items():
        vis_params = {
            "min": -1,
            "max": 1,
            "palette": ["white", blue_color]
        }
        m.add_ee_layer(ndwi, vis_params, key)
# レイヤーコントロールを追加して、各レイヤーの表示/非表示を切り替え可能にする
folium.LayerControl().add_to(m)
# 地図の表示
m

以下のようにかなり正しくメコン川の対象地域の（バッファーで計算した地域）の水面がとられられていそうです！

定量化する！

地図に出すだけだとまあ、先ほどのメコン川と水領域などの定性的チェックには役立ちますが、それだけではないも言えません。たたおえば水域の変化が必要であればもちろん地図を目で見比べるなんてこともできますが大変ですし数値がないとなにも言えません。なので地図だけでは不十分なのです。
そこで、水域の面積を計算してみました！NDVIを指標をもとめ、そして水である基準の０以上の部分でマスクします。

#対象地域の水指標を出す（陸地はマイナス、水域はプラス）
def calcNDWI(image):
    """
    NDWI = (Green - NIR) / (Green + NIR)
    Landsat 5 TM では Green=B2, NIR=B4
    """
    ndwiImage = image.normalizedDifference(['B2', 'B4']).rename('NDWI')
    # もとの画像プロパティ（特に撮影日時）を引き継ぐ
    return ndwiImage.set('system:time_start', image.get('system:time_start'))
ndwi = landsat5.map(calcNDWI)

# NDWI画像に対して、値が0以上の部分（水域）だけを残すマスクを適用する関数
def maskNDWI(image):
    return image.updateMask(image.gte(0))
# 各画像に対してマスクを適用（ImageCollection全体に map() で適用）
ndwi_masked = ndwi.map(maskNDWI)

上記ndwi_maskedの情報を使って水域面積を計算する関数を作成します。

# 各画像から水域面積と撮影日を抽出する関数
def computeWaterArea(image):
    # すでに ndwi_masked は値が0以上の部分のみ残っているので、
    # 再度 gte(0) を使ってバイナリマスクを作成（あるいはそのまま利用してもよい）
    water_mask = image.gte(0).rename('water_mask')
    # 各ピクセルの面積（平方メートル）を計算
    pixel_area = ee.Image.pixelArea()
    # 水域ピクセルの面積を算出する画像を作成
    water_area_image = water_mask.multiply(pixel_area).rename('water_area')
    # 指定領域（buffer_geometry）内の水域面積の総和を計算（scaleは30m）
    water_area_dict = water_area_image.reduceRegion(
        reducer=ee.Reducer.sum(),
        geometry=buffer_geometry,
        scale=30,
        maxPixels=1e9
    )
    # 水域面積（平方メートル）を取得。画像によっては値が存在しない場合があるので ee.Number() によるキャスト
    water_area = ee.Number(water_area_dict.get('water_area'))
    # 撮影日の取得（'system:time_start' プロパティがセットされている前提）
    date = ee.Date(image.get('system:time_start')).format('YYYY-MM-dd')
    # Feature として返す（プロパティに日付と水域面積をセット）
    return ee.Feature(None, {'date': date, 'water_area': water_area})

では関数を読んで実際に日付と水域面積を表データ化します。

# ndwi_masked は日付付きの NDWI のマスク済み ImageCollection
water_features_ndwi = ndwi_masked.map(computeWaterArea)
# FeatureCollection を Python 側にリストとして取得
water_list_ndwi = water_features_ndwi.getInfo()['features']
data = []
for f in water_list_ndwi:
    props = f['properties']
    # None チェック：水域面積が取得できない場合があるので
    if props['water_area'] is not None:
        # 水域面積は平方メートルなので、平方キロメートルに変換
        water_area_sqkm = props['water_area'] / 1e6
    else:
        water_area_sqkm = None
    data.append({'date': props['date'], 'water_area_sqkm': water_area_sqkm})

# Pandas DataFrame に変換
df_ndwi = pd.DataFrame(data)
print(df_ndwi)

以下のように撮影日と水域面積が得られました！
date water_area_sqkm
0 1987-12-16 40.842750
1 1988-01-17 40.371686
2 1988-02-02 40.595547
3 1988-04-06 55.946846
4 1988-05-24 0.054914
.. ... ...
149 1999-10-30 5.588682
150 1999-11-15 13.071911
151 1999-12-01 49.904220
152 1999-12-17 48.088657
153 1996-09-10 0.000000
[154 rows x 2 columns]

せっかくなので、グラフ化します。そのままのプロットと移動平均の折れ線を付けます。

# df_ndwi に 'date' 列と 'water_area_sqkm' 列がある前提
df = df_ndwi.copy()
df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d')
df.set_index('date', inplace=True)
df.sort_index(inplace=True)

# 移動平均（例：3シーンのローリング平均）を計算してプロット
df['rolling_mean'] = df['water_area_sqkm'].rolling(window=3, min_periods=1).mean()
plt.figure(figsize=(24,6))
plt.plot(df.index, df['water_area_sqkm'], marker='o', linestyle='', label='Water Area')
plt.plot(df.index, df['rolling_mean'], color='red', label='Rolling Mean (window=3)')
plt.title('Water Area with Rolling Mean')
plt.xlabel('Date')
plt.ylabel('Water Area (sq km)')
plt.legend()
plt.tight_layout()
plt.show()

以下のように対象地域の水域の面積変化のグラフができました！