衛星データ分析の基礎

1. はじめに

初めまして。データアナリティクスラボ株式会社データソリューション事業部の山本と申します。本記事では、社内活動の一環として学んできた衛星データの取得方法や取り扱い方、具体的な分析方法を衛星データ分析の基礎として整理した形でご紹介します。衛星データの取得方法や基本的な分析の仕方などを知りたい方はぜひご一読ください。

2. 衛星データについて

衛星データの取得方法や分析方法に入る前に、衛星データの概要について簡単に触れておきます。

2.1. 概要

2.1.1. 衛星データの種類

衛星データとは、地球の観測やモニタリングを目的として人工衛星から取得されるデータのことです。主に以下の2種類に分類されます。

1. ラスターデータ

• 内容: 衛星画像や地表面のデジタルデータで、ピクセル（セル）の集合として表現されます。
• 種類:
  • 光学衛星データ: 可視光、赤外線、マルチスペクトルデータなど。
  • レーダー衛星データ: マイクロ波を使った観測で、天候や光の影響を受けにくい。
• 用途: 土地利用、農業、環境モニタリング、災害対策など。
  
  　　　　　　　　　　　　　引用：国土地理院　地理院地図より

2. ベクターデータ

• 内容: 地物（建物、道路、河川など）を点、線、ポリゴンで表現したデータ。
• 用途: 地図作成、都市計画、GIS解析など。
  
  　　　　　　　　　　　　　引用：国土地理院　地理院地図より

2.1.2. 画像処理の必要性

ラスターデータ形式の衛星画像は基本的に使用する波長ごとに分かれています。ですので、例えば無加工のままではフルカラー画像として見ることができません。従って、衛星画像で分析をする際は、分析に適した画像に加工する技術が求められます。

衛星データについてより詳しく知りたい方は、是非こちらのサイトをご参照ください。

2.2. 公開中の衛星データ

2024年12月現在、インターネットでは様々な衛星データが公開されています。主に無料で公開されているオープンデータと、購入が必要な商用データの2つに分けられます。
ここでは代表的なサービスをいくつか紹介します。

1. 無料の衛星画像データ

無料で公開中の衛星画像データは研究や教育目的で使用されることが多く、特にNASAやESAといった宇宙機関が提供するデータがよく利用されます。

1.1. NASAのEarthdata

• プラットフォーム: NASA Earthdata（ https://www.earthdata.nasa.gov/ ）
• 特徴: NASAが運営するEarthdataでは、Landsatシリーズ、MODIS（Moderate Resolution Imaging Spectroradiometer）、VIIRS（Visible Infrared Imaging Radiometer Suite）などのデータを無料で提供しています。

1.2. USGS Earth Explorer

• プラットフォーム: USGS Earth Explorer（ https://earthexplorer.usgs.gov/ ）
• 特徴: Landsatシリーズの衛星画像データが中心です。Landsatデータは1970年代からの長期間のデータが利用可能です。

1.3. Sentinel Hub

• プラットフォーム: Sentinel Hub（ https://www.sentinel-hub.com/ ）
• 特徴: ヨーロッパ宇宙機関（ESA）のSentinel-1、Sentinel-2などの衛星データを提供しています。Sentinel-2は高解像度のマルチスペクトルデータを無料で利用できます。
  いずれのサービスも無料のアカウント登録が必要になります。検索ツールでエリアや時期を指定してダウンロードできます。なお、今回使用するこちらのSentinel HubはAPIを使用してデータをプログラムで取得することも可能です。

2. 商用の衛星画像データ

より高解像度のデータや最新の衛星画像が必要な場合は、商用のデータプロバイダーから購入する必要があります。

2.1. Maxar (DigitalGlobe)

• プラットフォーム: Maxar Technologies（ https://www.maxar.com/ ）
• 特徴: WorldView-3やWorldView-4などの高解像度衛星を運用しており、30cmクラスの超高解像度データを提供しています。地理空間分析、都市計画、災害対応など、商用利用で広く使われます。
  取得方法: 商用ライセンスが必要で、データの購入には問い合わせが必要です。

2.2. Planet Labs

• プラットフォーム: Planet Labs（ https://www.planet.com/ ）
• 特徴: Dove衛星群によって毎日地球全体をカバーする高頻度のデータを提供しています。農業、森林管理、環境監視など、迅速な変化を捉える用途で利用されています。
  取得方法: サブスクリプションモデルでデータを提供しており、利用目的に応じてライセンス契約を行います。

2.3. その他のプラットフォーム

• Google Earth Engine（ https://earthengine.google.com/ ）:
  大量の地理空間データをクラウド上で処理できるプラットフォーム。LandsatやSentinelなどの衛星データを検索、分析できます。
• Copernicus Open Access Hub（ https://scihub.copernicus.eu/ ）:
  ESAのCopernicusプログラムにより、Sentinelデータの無料アクセスを提供しています。

なお、今回はSentinel HubのAPIを利用する方法でSentinel-2の衛星画像を取得してみます。

3. Sentinel衛星画像の取得方法

3.1. Sentinel衛星シリーズについて

まずは、今回利用する衛星画像を撮影しているSentinel衛星シリーズについて簡単に触れておきます。
Sentinel衛星シリーズは、欧州宇宙機関（ESA）と欧州連合（EU）が主導する地球観測プログラム「Copernicus計画」の一部で、地球の環境監視や災害管理、気候変動の研究に広く利用されています。
Sentinel-1、Sentinel-2、Sentinel-3、Sentinel-5P、Sentinel-6 Michael Freilich（以下Sentinel-6）の5種類の衛星が現在運用中で、それぞれ観測方法や観測対象、用途に特徴があります。

	Sentinel-1	Sentinel-2
観測手法	C-SAR	光学
観測対象	陸域・海域	陸域
主用途	地表面形状　陸域画像取得 海面画像取得　海上風（水平方向の風向と風速）、 海面の流速、波高、海氷の分布・タイプを観測	地表面形状　陸域画像取得 土地被覆分類　植生変化 都市開発

	Sentinel-3	Sentinel-5P	Sentinel-6
観測手法	放射計	分光観測	放射計・高度計
観測対象	陸域・海域	空域	海域
主用途	海洋、陸地、河川、湖などの体系的計測 地球規模の気候現象のモニタリング	大気質の微量ガスやエアロゾルの計測	海面観測、気温や湿度などの気象情報の収集

3.2. Sentinel Hubのアカウント作成

ここからは、Sentinel衛星画像を取得する際に必要なアカウントの作成方法について説明します。
まず、こちらのサイトでユーザー登録を行います。

画面右上の緑色のボタンをクリック

画面右上の「SIGN IN」と書かれた緑色のボタンをクリック

画面中央下の「RESISTER」と書かれた緑色のボタンをクリック

必要情報を入力した後、登録したメールアドレスに送られてくるメールでユーザー認証を行えばアカウント登録は完了です。

次に、Sentinel HubをAPIで利用するために必要なOAuth認証の設定を行います。登録したアカウントでログインした後にメイン画面が表示されるので、以下の手順でOAuthの設定画面まで移動してください。

一番左の「Personal info」と書かれた箇所をクリック

画面左の「Sentinel Hub」と書かれた箇所をクリック

画面左下の「User settings」と書かれた箇所をクリック

画面右の「+ Create」と書かれた箇所をクリック

Client nameに任意の名前を入力する

入力が完了するとIDとKEYが発行されます。特にKEYは二度と見れなくなるので、安全な場所にメモしておきましょう。
なお、こちらのドキュメントにOAuthの詳しい設定方法が記載されていますので、不明な点がある際は適宜参照してください。

これで衛星画像データの取得に必要な準備が整いました。

4. Sentinel Hubの基本操作

ここからはSentinel Hubの基本操作について、最低限押さえるべき部分に絞って紹介します。

4.1. ライブラリのインポート

まず、Sentinel Hubを利用するために必要なライブラリをインポートします。
なお、この4章で記述されているコードはSentinel Hubの公式ドキュメントを参考にしています。より詳しい利用方法を知りたい方は参照ください。

# ライブラリのインポート
import cv2
import datetime
import glob
import japanize_matplotlib
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import requests
import shutil

import geopandas as gpd
from sentinelhub import (
SHConfig,
CRS,
BBox,
DataCollection,
DownloadRequest,
MimeType,
MosaickingOrder,
SentinelHubCatalog,
SentinelHubDownloadClient,
SentinelHubRequest,
bbox_to_dimensions,
Geometry
)
from rasterio.plot import show
from rasterio.mask import mask
import rasterio as rio

4.2. ユーザー認証

まず、登録したアカウントの資格情報を設定するユーザー認証というものを行う必要があります。
なお、OAuth認証はデフォルトだと90日間で期限を迎えるため、作成後90日以上経った場合は再度Oauth認証を作成してからユーザー認証を行ってください。

# 資格情報を設定（初回のみ）
config = SHConfig()
config.sh_client_id = "your client id"
config.sh_client_secret = "your secret key"
config.sh_token_url = "https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token"
config.sh_base_url = "https://sh.dataspace.copernicus.eu"
config.save("cdse")

# ２回目以降は下記のコードを実行するだけで良い
config = SHConfig("cdse")

# idもしくはkeyが設定されていない場合に警告文を表示
if not config.sh_client_id or not config.sh_client_secret:
    print("Warning! To use Process API, please provide the credentials (OAuth client ID and client secret).")

4.3. 読み込み範囲の指定

続いて取得したい領域情報を設定します。関心領域の座標を設定したのち、bboxの初期化と設定を行います。
例えば、日本の南関東あたりを設定したい場合は以下のように記述します。

# 関心領域の設定
aoi_coords_wgs84 = (139.589311,35.548802,139.948770,35.775579)
# 解像度を設定
resolution = 30
# 関心領域の初期化・設定
aoi_bbox = BBox(bbox=aoi_coords_wgs84, crs=CRS.WGS84)
aoi_size = bbox_to_dimensions(aoi_bbox, resolution=resolution)
# 解像度とピクセルサイズを表示
print(f'Image shape at {resolution} m resolution: {aoi_size} pixels')

4.4. データの取得

ここからいよいよ衛星画像データの取得になります。
衛星画像データを取得するには、sentinelhubライブラリのSentinelHubRequest関数を利用します。関数の利用にあたり、まず必須となるevalscriptを設定します。evalscriptとは、Sentinel Hubによる衛星データの処理方法とサービスが返す値を定義するJavascriptのコードになります。

# evalscriptの設定
evalscript_true_color = """
    //VERSION=3
    function setup() {
        return {
            input: [{
                bands: ["B02", "B03", "B04"]
            }],
            output: {
                bands: 3
            }
        };
    }
    function evaluatePixel(sample) {
        return [sample.B04, sample.B03, sample.B02];
    }
"""

このevalscriptでは、12あるバンドのうちB02とB03、B04の3種類のバンドを取得を要求することを表しています。このB02とB03、B04の３種類のバンドを加工することでカラー画像を出力することができます。

続いて、SentinelHubRequest関数を利用してリクエストを作成します。

# リクエストを作成
request_true_color = SentinelHubRequest(
    evalscript=evalscript_true_color,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A.define_from(
                name="s2", 
                service_url="https://sh.dataspace.copernicus.eu"
            ),
            time_interval=('2020-08-10', '2020-08-15')
        )
    ],
    responses=[SentinelHubRequest.output_response("default", MimeType.PNG)],
    bbox=aoi_bbox,
    size=aoi_size,
    config=config,
)

作成したリクエストをもとに、関心領域の画像データを取得・表示します。

# 画像データの取得
true_color_imgs = request_true_color.get_data()

# 正しく画像データを取得できているかを確認(正しく取得できていればnumpy配列の形式かつ長さ1のリストで返される)
print(f"Returned data is of type = {type(true_color_imgs)} and length {len(true_color_imgs)}.")
print(f"Single element in the list is of type {type(true_color_imgs[-1])} and has shape {true_color_imgs[-1].shape}")

# 画像の表示
image = true_color_imgs[0]
print(f"Image type: {image.dtype}")
# plot function
# factor 1/255 to scale between 0-1
# factor 3.5 to increase brightness
plot_image(image, factor=3.5 / 255, clip_range=(0, 1))

取得した画像がこちら

5. 画像データの加工

具体的な分析例に移る前に、この章では画像データの前処理方法や便利な加工方法について紹介します。というのも、これらの前処理や加工方法は、波長ごとに異なる衛星画像データを複雑に組み合わせて処理を行う衛星データ分析で大いに活用できる技術だからです。
まず前準備としてsentinelhubライブラリを用いて全バンドの画像データを取得し、ローカルに保存しておきます。以下のコードは全バンド取得〜データ保存までの一連の処理を実行するものになります。ユーザー認証と読み込み範囲の指定は設定済とします。

# evalscriptを作成
evalscript_all_bands = """
    //VERSION=3
    function setup() {
        return {
            input: [{
                bands: ["B01","B02","B03","B04","B05","B06","B07","B08","B8A","B09","B11","B12","AOT"],
                units: "DN"
            }],
            output: {
                bands: 13,
                sampleType: "INT16"
            }
        };
    }

    function evaluatePixel(sample) {
        return [sample.B01, sample.B02, sample.B03, sample.B04,
                sample.B05, sample.B06, sample.B07, sample.B08,
                sample.B8A, sample.B09, sample.B11, sample.B12,
                sample.AOT];
    }
"""

# 保存先フォルダの設定
FOLDER_NAME = "sample_all_bands"
# リクエストを作成
request_all_bands = SentinelHubRequest(
    data_folder=FOLDER_NAME,
    evalscript=evalscript_all_bands,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A.define_from(
                name="s2", 
                service_url="https://sh.dataspace.copernicus.eu"
            ),
            time_interval=('2020-08-10', '2020-08-15'), 
            other_args={"dataFilter": {"mosaickingOrder": "leastCC"}},
        )
    ],
    responses=[SentinelHubRequest.output_response("default", MimeType.TIFF)],
    bbox=aoi_bbox,
    size=aoi_size,
    config=config,
)

# データ取得を実行
all_bands_img = request_all_bands.get_data(save_data=True)

5.1. バンドの合成

バンドの合成という処理は、特にリモートセンシングや画像処理の分野でよく使われます。この手法は簡単に言うと、画像の色や情報を構成するバンドという層を組み合わせて、新しい画像を作り出す処理になります。
例えば、カラー画像は通常、赤、緑、青の3つのバンドで構成されます。これらを組み合わせることで、普段目にするようなカラー画像を表示することができます。衛星画像では可視光以外のバンド（赤外線など）を含むことがあり、それらを合成することで植生の健康状態や水の分布といった特定の情報を視覚化することができます。

それでは、実際にバンドの合成を行って画像を表示してみましょう。以下が実行コードになります。

# ファイルの読み込み
FILE = glob.glob(f"./{FOLDER_NAME}/*/*.tiff")[0]
dataset = rio.open(FILE)

# 青色（バンド2）, 緑色（バンド3）, 赤色（バンド4）の読み込み
band_02 = src.read(2)
band_03 = src.read(3)
band_04 = src.read(4)
# 各バンドのデータをスタックしてRGB画像を作成
bgr_image = np.dstack((band_02, band_03, band_04))

# 0-255の範囲に正規化
def normalize(image):
    min_val = np.min(image)
    max_val = np.max(image)
    
    return ((image - min_val) / (max_val - min_val) * 255).astype(np.uint8)

bgr_image = normalize(bgr_image)

# カラー画像を表示
rgb_image = cv2.cvtColor(bgr_image, cv2.COLOR_BGR2RGB)

plt.figure(figsize=(10, 10))
plt.imshow(rgb_image)
plt.title("RGB Image (Band 02, 03, 04)")
plt.axis('off')
plt.show()

カラー画像として表示することはできましたが、全体的に暗くこのままでは見辛くなってしまっています。

5.2. 画像の明度調整

前節で作成したカラー画像をより見やすくする方法として、ガンマ補正という手法を紹介します。
そもそも作成した画像が暗く見にくくなっているのは、色味を表すピクセル値が人間の目の感覚に適した形で扱われておらず、線形的に表されているためです。ガンマ補正を行うことで、ピクセル値を非線形に調整して人間の視覚に近づいた形に変換することができます。以下、実行コードになります。

# 関数定義
def gamma_correction(image, gamma):
    '''
    ガンマ補正を行う
    '''
    # 整数型で2次元配列を作成[256,1]
    lookup_table = np.zeros((256, 1), dtype = 'uint8')
    for i in range(256):
        # γテーブルを作成
        lookup_table[i][0] = 255 * (float(i) / 255) ** (1.0/gamma)

    # lookup Tableを用いて配列を変換
    return cv2.LUT(image, lookup_table)

    
# CLAHE
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(3, 3))
# ガンマ補正を実施
gamma_image = gamma_correction(bgr_image, 2)
# CLAHEを使ってV（明度）を平坦化
gamma_image_hsv = cv2.cvtColor(gamma_image, cv2.COLOR_BGR2HSV)
gamma_image_hsv[:, :, 2] = clahe.apply(gamma_image_hsv[:, :, 2])


# ガンマ補正後の画像を表示
gamma_image_bgr = cv2.cvtColor(gamma_image_hsv, cv2.COLOR_HSV2BGR)
gamma_image_rgb = cv2.cvtColor(gamma_image_bgr, cv2.COLOR_BGR2RGB)
plt.figure(figsize=(10, 10))
plt.imshow(gamma_image_rgb)
plt.show()

先ほどよりも全体的に明るくなり、細部まで見えるようになりました。

5.3. その他の加工方法

ここからは、そのほかの便利な画像処理方法を3つほど紹介します。

5.3.1. 画像の切り抜き

まずは画像を読み込みます。使用するのは赤色のドライバーを撮った写真データです。

img1 = cv2.imread("./example_data/driver.jpeg")
# 画像の表示
plt.imshow(cv2.cvtColor(img1, cv2.COLOR_BGR2RGB))

表示した画像の左側と下側にある目盛りを参考にして切り抜きたい範囲をスライスで指定します。この例ではドライバーが写っている部分を抜き出したいので、heightは200~2800、widthは1000~3000を範囲として指定しています。

# スライスによる切り抜き
img1_crop = img1[200:2800, 1000:3000, :]
# 切り抜いた画像の表示
plt.imshow(cv2.cvtColor(img1_crop, cv2.COLOR_BGR2RGB))

うまくドライバーが写る部分だけを抜き出すことができました。
この画像は次の節でも用いるので保存しておきます。

# crop後の画像を保存
cv2.imwrite("./example_data/driver_crop.png", img1_crop)

5.3.2. 特定の色合い部分の抽出

続いては画像内から特定の色合い部分を抽出する方法です。ここではHSV変換を利用するため、まずはHSV色空間について簡単にご紹介します。
HSV色空間は色相 (Hue)、彩度 (Saturation, Chroma)、明度 (Value, Brightness) の三つの成分からなる色空間です。色相は色の種類（例えば、赤、青、黄色）のことで、通常0–360° の範囲で表されます。彩度は色の鮮やかさのことです。彩度が低下するにつれて灰色さが顕著になり、くすんだ色が現れます。明度は色の明るさのことです。（Wikipediaより）
以下はH(色相)の値と色の関係をグラフとして表示するコードその表示結果になります。左側の0に近い値では赤っぽく、そこから右に向かって値が大きくなると黄→緑→青→紫と変化していき、最終的に赤に戻ることがわかるかと思います。

# HSVの値を入れるarray配列を作成
h_map = np.zeros((60, 180, 3), dtype=np.uint8)

# H(色相)は0~179までの連続する整数
for i in range(180):
    h_map[:, i, 0] = i
# S(彩度)は全て255
h_map[:, :, 1] = 255
# V(明度)は全て255
h_map[:, :, 2] = 255

# HSVからBGRに変換
h_map = cv2.cvtColor(h_map, cv2.COLOR_HSV2BGR)

# H(色相)の値と色の関係をグラフに表示
fig = plt.figure(figsize=(3, 1), dpi=200)
plt.imshow(cv2.cvtColor(h_map, cv2.COLOR_BGR2RGB))

それでは、このHSVの特性を用いて特定の色合い部分の抽出を行います。試みるのは画像内のドライバーの持ち手部分です。まず切り抜き後画像をBGR形式からHSV型式に変換した上で、H、S、Vそれぞれ下限値と上限値を設定し、各ピクセルが下限値と上限値の間に入るか否かをcv2のinRange関数で判定させ二値化を行います。ドライバーの持ち手は赤色なので、先ほどのH(色相)と色の関係グラフを参考に下限値を0、上限値を15に設定してみます。

# 切り抜き後画像の読み込み
img2 = cv2.imread('./example_data/driver_crop.png')

# HSV形式に変換
img2_hsv = cv2.cvtColor(img2, cv2.COLOR_BGR2HSV)

# H, S, Vそれぞれの下限値と上限値を設定
lower_red = np.array([0, 100, 100])
upper_red = np.array([15, 255, 255])
# 各ピクセルが上限値と下限値の間に入るか否かで二値化を行う
img2_mask = cv2.inRange(img2_hsv, lower_red, upper_red)

続いてマスキング用データを作成し、二値化データ(img2_mask)を用いてBGR全てに対して二値化を実行します。

# 画像の高さ、幅、層数を変数に代入
img2_width, img2_height, img2_col = img2_hsv.shape
# マスキング用データを作成
img2_mask_img = np.zeros((img2_width, img2_height, 3), dtype=np.uint8)
# BGR全てで二値化を実行
for i in range(3):
    img2_mask_img[:, :, i] = img2_mask

# 二値化した画像を表示
plt.imshow(cv2.cvtColor(img2_mask_img, cv2.COLOR_BGR2RGB))
plt.show()

最後に、元の画像と二値化した画像をもとに、cv2のbitwise_and関数を用いてマスキングを施します。

# 処理前画像と二値化画像でマスキングを行う
img2_masked = cv2.bitwise_and(img2, img2_mask_img)

# マスキングした画像を表示
plt.imshow(cv2.cvtColor(img2_masked, cv2.COLOR_BGR2RGB))
plt.show()

概ね赤色の持ち手部分だけを抽出することができました。

5.3.3. 輪郭の検出

最後は輪郭を検出する方法です。今回紹介する方法では前節で作成したマスキング後の画像を用います。
まずはマスキング後の画像から輪郭の元になるものを抽出します。具体的には、グレースケール変換と二値化処理を行います。

# グレースケールに変換
img2_masked_gray = cv2.cvtColor(img2_masked, cv2.COLOR_BGR2GRAY)
# グレースケール画像を二値化
ret, img2_masked_gray = cv2.threshold(img2_masked_gray, 0, 255, cv2.THRESH_BINARY)

# 処理後画像を表示
plt.imshow(img2_masked_gray, cmap="gray")
plt.show()

続いて輪郭の抽出を行います。cv2のfindContours関数を用いて二値画像から輪郭を検出した後、面積の小さい輪郭データを除外しつつcv2のpolylines関数で輪郭線を描画します。

# 二値画像から輪郭を検出
contours, hierarchy = cv2.findContours(img2_masked_gray,       # グレースケール画像
                                       cv2.RETR_LIST,          # 全輪郭を検出
                                       cv2.CHAIN_APPROX_SIMPLE # 輪郭を簡略化
                                       )

# 元画像をコピー
img2_draw_line = np.copy(img2)

# 輪郭データの選別と輪郭線の描画
# 面積の小さい輪郭データが混じっているので、ループ処理を用いて選別する
for i in range(len(contours)):
    # 画面が大きくないものは無視する
    contour = contours[i]
    c_area = cv2.contourArea(contour)
    if c_area < 1000:
        continue
    # 長方形で囲む
    cv2.polylines(img2_draw_line, # 輪郭を描画する画像
                  contours[i],    # 描画する輪郭
                  True,           # 輪郭を閉じる指定
                  (255, 0, 0),    # 輪郭線の色
                  8               # 輪郭線の太さ
                 )

# 輪郭線の入った画像を表示
plt.imshow(cv2.cvtColor(img2_draw_line, cv2.COLOR_BGR2RGB))
plt.show()

補足：Canny法による輪郭の検出

輪郭（エッジ）を検出する別の方法としてCanny法もご紹介します（6章の分析例ではこの方法を使っています）。
Canny法では画像の高度変化が大きい部分を輪郭として検出することができます。この方法では二つの閾値を調整することが必要です。両者の役割は異なっており、一方はエッジがどうかを判断し、もう一方はエッジが続いているかどうかを判断する役割を持っています。詳細な考え方および計算方法はOpenCVの公式ドキュメントを参照してください。

以下はCanny法を用いた場合のコードです。

# 二つの閾値を設定
threshold1 = 150
threshold2 = 350
# Canny法で画像を処理
img2_edged = cv2.Canny(img2, threshold1, threshold2)
# 処理後の画像を表示
plt.imshow(img2_edged, cmap="gray")
plt.show()

検出した輪郭を膨張させて目立たせることもできます。膨張用の行列データを作成した上で、cv2のdilate関数を用いて輪郭を膨張させます。

# エッジ膨張用の行列を作成
kernel = np.ones((3, 3)) / 9.0
# エッジを膨張
img2_edge_dilated = cv2.dilate(img2_edged, kernel)
# 膨張させたエッジの画像を表示
plt.imshow(img2_edge_dilated, cmap="gray")
plt.show()

周囲にノイズが見られつつも、ドライバーの輪郭を検出できていることが確認できます。

6. 分析例

ここからは、実際の衛星データを用いて4章および5章で紹介した衛星データの取得方法や画像加工の方法を駆使しながら分析を行ってみます。
分析テーマは「地震活動による海岸線の変化」です。今回の分析例では2024年1月に石川県沖で発生した能登半島地震を題材とします。能登半島を撮影した異なる時期の衛星データを取得し、一連の画像処理を経て半島の海岸線が地震の発生前後でどのように変化しているかを可視化してみます。

6.1. データ準備

分析に入る前に、まずは利用する衛星データを取得します。下記のコードでは能登半島を対象地域に、また期間は地震発生の前と後である必要があるため、2023-12-01~31(1枚目)と2024-04-01~10-31(2枚目)に設定しています。
取得方法は4章でご紹介した通りですが、今回は一部のバンドだけでなく全バンドを、さらにtiff形式でローカルに落としている影響でコードの中身が若干異なるのでご注意ください。
なお、ライブラリのインポートとユーザ認証は実行済を前提としています。

# 解像度を設定
resolution = 30
# 関心領域の初期化・設定
aoi_coords_wgs84 = (136.644229,37.248371, 136.938113,37.439556)
aoi_bbox = BBox(bbox=aoi_coords_wgs84, crs=CRS.WGS84)
aoi_size = bbox_to_dimensions(aoi_bbox, resolution=resolution)

print(f'Image shape at {resolution} m resolution: {aoi_size} pixels')

# evalscriptの設定
evalscript_all_bands = """
    //VERSION=3
    function setup() {
        return {
            input: [{
                bands: ["B01","B02","B03","B04","B05","B06","B07","B08","B8A","B09","B11","B12","AOT"],
                units: "DN"
            }],
            output: {
                bands: 13,
                sampleType: "INT16"
            }
        };
    }

    function evaluatePixel(sample) {
        return [sample.B01,
                sample.B02,
                sample.B03,
                sample.B04,
                sample.B05,
                sample.B06,
                sample.B07,
                sample.B08,
                sample.B8A,
                sample.B09,
                sample.B11,
                sample.B12,
                sample.AOT];
    }
"""
# 
for folder_name, time_interval in zip(["noto_before", "noto_after"], 
                                      [["2013-12-01", "2023-12-31"], 
                                       ["2024-04-01", "2024-10-31"]]):
    # リクエストの設定
    request_all_bands = SentinelHubRequest(
        data_folder=folder_name, # 保存先のフォルダ名を設定
        evalscript=evalscript_all_bands,
        input_data=[
            SentinelHubRequest.input_data(
                data_collection=DataCollection.SENTINEL2_L2A.define_from(
                    name="s2", 
                    service_url="https://sh.dataspace.copernicus.eu")
                ,
                time_interval=(time_interval[0], time_interval[1])
                , 
                other_args={"dataFilter": {"mosaickingOrder": "leastCC"}} 
                ,
            )
        ],
        responses=[SentinelHubRequest.output_response("default", MimeType.TIFF)],
        bbox=aoi_bbox,
        size=aoi_size,
        config=config,
    )
    # データ取得を実行
    all_bands_img = request_all_bands.get_data(save_data=True)
    # 取得結果の確認
    print(
        "The output directory has been created and a tiff file with all 13 bands was saved into the following structure:\n"
    )

    for folder, _, filenames in os.walk(request_all_bands.data_folder):
        for filename in filenames:
            print(os.path.join(folder, filename))

6.2. データの格納

前節で取得した画像データを読み込みんでいきます。なお、この分析例ではデータの扱いやすさを考慮して、クラスにデータを保持する形を取ります。

# 衛星データを管理するクラスを作成
class Region:
    def __init__(self, path, bands=['B02','B03','B04','B08']):
        self.path = path
        self.bands = bands
        self.dataset = self._load_data(path, bands)

    def _load_data(self, path, bands):
        """
        データを読み込んで辞書形式で返す
        """
        # バンド名とリストのインデックスを対応させる辞書を作成
        band_dict = {"B01": 1, "B02": 2, "B03": 3, "B04": 4, 
                     "B05": 5, "B06": 6, "B07": 7, "B08": 8, 
                     "B8A": 9, "B09": 10, "B11": 11, "B12": 12, 
                     "AOT": 13,}
        return {b:rio.open(glob.glob(path)[0]).read(band_dict[b]) for b in bands}

    def _gamma_correction(self, image, gamma=3.0):
        """
        データに対してガンマ補正を行う
        """
        # 整数型で2次元配列を作成[256,1]
        lookup_table = np.zeros((256, 1), dtype = 'uint8')
        for i in range(256):
            # γテーブルを作成
            lookup_table[i][0] = 255 * (float(i) / 255) ** (1.0/gamma)
    
        # lookup Tableを用いて配列を変換
        return cv2.LUT(image, lookup_table)

    # 与えられたデータをもとにカラー画像を返す
    def get_rgb_img(self, dataset):
        # B02,03,04を合成したBGR形式の画像を用意する
        bgr = np.zeros((dataset['B02'].shape[0], dataset['B02'].shape[1], 3), dtype=np.uint8)

        for i, b in enumerate(self.bands[:3]):
            bgr[:, :, i] = dataset[b] // 256

        # HSV形式に変換する
        hsv = cv2.cvtColor(bgr, cv2.COLOR_BGR2HSV)
        # V（明度）をB08のデータで置き換える
        hsv[:, :, 2] = dataset['B08'] // 256
        # Vを置き換えたデータを再びBGR形式に置き換える
        rgb = cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB)

        # ガンマ補正を実施
        return self._gamma_correction(rgb)

このクラスは、画像の明度調整（ガンマ補正）を_gamma_correction関数で、バンドの合成をget_rgb_image関数で行う仕様になります。

では、実際にデータを格納します。

# 衛星データの取得
region = []
region.append(Region("./noto_before/*/*.tiff"))
region.append(Region("./noto_after/*/*.tiff"))

ちゃんと格納されているかどうか、確認も兼ねてカラー画像を表示してみます。

# 2023年と2024年の画像を表示する
for n in [0, 1]:
    plt.figure(figsize=(10,10))
    plt.imshow(region[n].get_rgb_img(region[n].dataset))
    plt.show()

　　　　　　　　　　　　　　　　　　　　2023年の画像

　　　　　　　　　　　　　　　　　　　　2024年の画像

ちゃんとデータが格納されていることが確認できました。

6.3. クラスタリングによる陸/水域の分類

領域データをクラスタリングを用いて陸域と水域に分類していきます。
前準備として、クラスタリングを行えるようにデータ構造を変換します。現状ではバンドごとに（高さ×幅）の形でデータが保持されていますが、これを（バンド数×領域のピクセル数）のNumpy配列に変換後、さらに転置して最終的に（領域のピクセル数×バンド数）の形にします。

for r in region:
    # 領域の辞書をNumpy配列に変換する
    tmp = np.array([r.dataset[b] for b in r.bands])
    print(tmp.shape)
    # 2次元にして転置する（行と列を入れ替える）
    r.region_reshaped = tmp.reshape(len(tmp), -1).T
    print(r.region_reshaped.shape)

データ構造を変換できたらクラスタリングを行います。今回の分析例ではk-means++法と言うアルゴリズムで実施します。

# ライブラリのインポート
from sklearn.cluster import KMeans

# 2023年と2024年のデータを結合
train = np.append(region[0].region_reshaped, region[1].region_reshaped, axis=0)
# 学習の実施
model = KMeans(n_clusters=2, init='k-means++', n_init=10)
model.fit(train)

# 領域の大きさを取得
height = region[0].dataset['B02'].shape[0]
width = region[0].dataset['B02'].shape[1]
print(height, width)

# クラスタリングを実施し、結果をregion_predictに代入する
for r in region:
    r.region_predict = model.predict(r.region_reshaped).reshape(height, width)
# クラスタリング結果を確認する
for n in [0, 1]:
    plt.figure(figsize=(10,10))
    plt.imshow(region[n].region_predict, cmap="tab20_r")
    plt.show()

　　　　　　　　　　　　　　　　2023年　クラスタリングした結果

　　　　　　　　　　　　　　　　2024年　クラスタリングした結果

クラスタリング後画像の陸域にはノイズが生じています。このままの状態だとこの後行う海岸線の検出が上手くできない可能性があるため、できる限りノイズの除去を行います。ノイズの除去にはcv2のGaussianBlur関数を用います。

# ノイズを除去し、regionオブジェクトのregion_blurred変数に格納
for r in region:
    tmp = r.region_predict.astype("uint8")
    r.region_blurred = cv2.GaussianBlur(tmp, (5, 5), 0)
# ノイズ除去後の画像を表示
for n in [0, 1]:
    plt.figure(figsize=(10,10))
    plt.imshow(region[n].region_blurred, cmap="tab20_r")
    plt.show()

　　　　　　　　　　　　　　　　　2023年　ノイズ除去した結果

　　　　　　　　　　　　　　　　　2024年　ノイズ除去した結果

6.4. 海岸線の検出

ノイズ除去後の画像を用いて、海岸線（エッジ）を検出します。エッジの検出には5章で紹介したCanny法を用いています。

# エッジ膨張用の行列データを作成
kernel = np.ones((1, 1), np.uint8)

for r in region:
    # Canny法でエッジを検出
    tmp = cv2.Canny(r.region_blurred, 2, 5).astype("float")
    # エッジを膨張させる
    tmp = cv2.dilate(tmp, kernel, iterations=1)
    # region_edge変数に格納　
    r.region_edge = tmp

# 検出したエッジの画像を表示
for n in [0, 1]:
    plt.figure(figsize=(15, 15))
    plt.imshow(region[n].region_edge, cmap="Set3_r")
    plt.show()

6.5. 検出した海岸線の比較

最後に検出した二つの海岸線を元の画像と重ね合わせ、二つの海岸線を比較します。

plt.figure(figsize=(15, 15))
# 元画像
plt.imshow(region[0].get_rgb_img(region[0].dataset))

# 検出した海岸線
# 2023年（ピンク）
edge_image_before = np.zeros((height, width, 4), dtype=np.uint8)
edge_image_before[:, :, 0] = region[0].region_edge * 40  # 赤チャンネル
edge_image_before[:, :, 1] = region[0].region_edge * 200  # 緑チャンネル
edge_image_before[:, :, 2] = region[0].region_edge * 40  # 青チャンネル
edge_image_before[:, :, 3] = (region[0].region_edge > 0).astype(np.uint8) * 255  # 透明度チャンネル
plt.imshow(edge_image_before, cmap="Set3_r")

# 2024年（シアン）
edge_image_after = np.zeros((height, width, 4), dtype=np.uint8)
edge_image_after[:, :, 0] = region[1].region_edge * 200  # 赤チャンネル
edge_image_after[:, :, 1] = region[1].region_edge * 40  # 緑チャンネル
edge_image_after[:, :, 2] = region[1].region_edge * 40  # 青チャンネル
edge_image_after[:, :, 3] = (region[1].region_edge > 0).astype(np.uint8) * 255  # 透明度チャンネル
plt.imshow(edge_image_after, cmap="Set1")

plt.show()

上記の画像だとわかりにくいので、一部を抜き出して詳しく見てみます。

# 一部領域を表示
plt.figure(figsize=(10,10))

plt.imshow(region[0].get_rgb_img(region[0].dataset)[200:500, 150:450])
plt.imshow(edge_image_before[200:500, 150:450], cmap="Set3_r")
plt.imshow(edge_image_after[200:500, 150:450], cmap="Set1")

plt.show()

上記の画像を見てみると、2024年（シアン）の海岸線が2023年（ピンク）よりも海側を走っている部分があります。地面が隆起したことで海岸線が変化したことが考えられます。ここまでの一連の分析によって、部分的に海岸線が変化し前進していることがわかりました。

7. おわりに

本記事では、衛星データの取得から基本的な分析手法、さらには地震による海岸線の変化を可視化する分析例までを一通り解説しました。衛星データの扱いに初めて触れる方でも、この記事が具体的な手順や活用のヒントを得る一助になれば幸いです。
もし実際に試された結果や新たな疑問が生じた際は、ぜひフィードバックをお寄せください。
また、今後は異なる衛星データの取得方法や、具体的な衛星データの分析事例についてもさらに掘り下げた記事を執筆していく予定です。

最後までお読みいただき、ありがとうございました。

8. 参考文献

空畑 SORABATAKE,「衛星データのキホン~分かること、種類、頻度、解像度、活用事例~」,https://sorabatake.jp/279/
一般社団法人リモート・センシング技術センター,「Sentinel-2A/2B/2C/2D」,https://www.restec.or.jp/satellite/sentinel-2-a-2-b
Copernicus「Copernicus Data Space Ecosystem|Europe's eyes on Earth」,https://dataspace.copernicus.eu/
SINERGISE「Documentation of sentinelhub Python package」,https://sentinelhub-py.readthedocs.io/en/latest/index.html
田中康平,田村賢哉,玉置慎吾　著,宮﨑浩之　監修,『Pythonで学ぶ衛星データ解析基礎』,技術評論社,2022年