はじめに

本記事では衛星データを用いて、船舶の数をカウントできるか調査します。
結論から言うと、以下の無料データでカウントを行うにはいくつか課題があり、今回は精緻なカウントではなく2024年1月の能登地震前後での船舶数トレンドを確認することにしました。

利用環境：GoogleColaboratory
プラットフォーム：Google Earth Engine
データセット：COPERNICUS/S1_GRD
データタイプ：SAR

Pythonコードと実行結果

依存パッケージのインストール（Colab）

!pip install folium
!pip install geemap

解説

Google Colab 環境で必要なライブラリをその場でインストールしています。
folium: 地図描画
geemap: Earth Engine と Python/folium の橋渡し

ライブラリのインポート／EE初期化／Driveマウント

import ee
import numpy as np
import matplotlib.pyplot as plt
import geemap
import json
from datetime import datetime
from datetime import timedelta
import folium
import cv2
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.ticker as ticker
import math

try:
    ee.Initialize()
except Exception as e:
    ee.Authenticate()
    ee.Initialize(project='ご自身のプロジェクト名')

from google.colab import drive
drive.mount('/content/drive')

解説

Earth Engine (ee) や数値計算・可視化のライブラリを読み込みます。
ee.Initialize() で EE を初期化。失敗した場合は ee.Authenticate() → ee.Initialize(project='ご自身のプロジェクト名') で再初期化。
Colab 上の Google Drive を /content/drive にマウントして、入出力先として利用します。

関心領域（AOI）の定義

# 関心領域を設定
A =  {
  "type": "FeatureCollection",
  "features": [
    {
      "properties": {
        "note": "",
        "distance": "N/A",
        "drawtype": "rectangle",
        "area": "N/A"
      },
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [137.4068223,37.4608783],
            ...
            [138.0797349,37.1834869],
            [137.4068223,37.4608783]
          ]
        ]
      }
    }
  ]
}

#Google Driveのワークディレクトリを指定
output_dir = "/content/drive/MyDrive/Colab Notebooks/SAR/ToyamaWan-Senpaku/"

#今後使用する任意のファイル名をセットする． 例えば，地域の名前など．
object_name = 'ToyamaWan'

with open(output_dir + str(object_name) +'_2.geojson', 'w') as f:
    json.dump(A, f)

json_file = open(output_dir + str(object_name) +'_2.geojson')
json_object = json.load(json_file)

#jsonから関心域の緯度・経度情報のみを抽出する．
AREA = json_object["features"][0]["geometry"]['coordinates'][0]

"""関心域の確認
"""
m = folium.Map([(AREA[0][1]+AREA[len(AREA)-2][1])/2,(AREA[0][0]+AREA[len(AREA)-3][0])/2], zoom_start=8)
folium.GeoJson(output_dir + str(object_name) +'_2.geojson').add_to(m)
m

# region=ee.Geometry.Polygon(AREA)
region=ee.Geometry.Polygon(A["features"][0]["geometry"]['coordinates'][0])
region['coordinates'][0]

解説

関心領域（AOI）を GeoJSON 形式でスクリプト内に定義し、output_dir 配下に保存。
保存した GeoJSON を読み直し、座標配列 AREA を取り出します。
folium で AOI を地図上に可視化して確認。
ee.Geometry.Polygon(...) から Earth Engine の region を生成します。

結果
上記mの表示は以下であり、富山湾全体を選択していることがわかります。

create_daily_images（日別平均画像の生成）

# 日ごとの画像を生成
def create_daily_images(start_date, end_date, collection, region):
    results = []
    current_date = datetime.strptime(start_date, "%Y-%m-%d")
    end_date = datetime.strptime(end_date, "%Y-%m-%d")
    while current_date <= end_date:
        next_date = current_date + timedelta(days=1)
        daily_images = collection.filterDate(
            current_date.strftime("%Y-%m-%d"), next_date.strftime("%Y-%m-%d")
        ).map(lambda img: img.clip(region)).mean().focalMean(radius=1,kernelType='square')  # 日ごとに平均を取る

        results.append({
            'date': current_date.strftime("%Y-%m-%d"),
            'image': daily_images
        })
        current_date = next_date
    return results

解説

指定期間 [start_date, end_date] を 1 日刻みで走査し、毎日の画像を生成する関数
collection.filterDate(当日, 翌日) → clip(region) → mean() → focalMean(radius=1, kernelType='square') で日平均・平滑化画像を作成
※周期的に1日に2回以上撮影されることはないはずなのでmean()は無くても変わらないと思われる。focaMeanは周り1ピクセル分の平均化によりノイズを軽減させている。
返り値は {date: 'YYYY-MM-DD', image: ee.Image} のリスト

split_into_tiles（領域のタイル分割）

# regionをのtileに分割
def split_into_tiles(region, num_x, num_y):
    """region を num_x × num_y のタイルに分割（位置情報付き）"""
    bounds = region.bounds().coordinates().get(0).getInfo()  # 領域の境界座標を取得
    x_min, y_min = bounds[0]
    x_max, y_max = bounds[2]

    x_step = (x_max - x_min) / num_x
    y_step = (y_max - y_min) / num_y

    tiles_grid = []  # 各タイルの位置情報を保持するリスト（行 × 列）

    for row in range(num_y):  # 行方向
        row_tiles = []
        for col in range(num_x):  # 列方向
            x1 = x_min + col * x_step
            x2 = x_min + (col + 1) * x_step
            # 南半球用
            # y1 = y_min + row * y_step
            # y2 = y_min + (row + 1) * y_step
            # 北半球用
            y1 = y_max - row * y_step
            y2 = y_max - (row + 1) * y_step

            tile = ee.Geometry.Rectangle([x1, y1, x2, y2])
            row_tiles.append((tile, row, col))  # タイルとその位置情報をセットで保存

        tiles_grid.append(row_tiles)  # 行ごとに格納

    return tiles_grid  # 位置情報付きのタイルリストを返す

解説

region の外接矩形を取得し、num_x × num_y のタイルに分割する関数（北半球の扱いに対応）
※南半球の場合は追加する方向が逆になるのでご注意
返り値は [(tileGeom, row, col), ...] を行単位の二次元配列にしたグリッド

ship_count（船舶抽出〜可視化〜集計）

def ship_count(from_date, to_date):
  Sentinel1 = (ee.ImageCollection('COPERNICUS/S1_GRD')
             .filterBounds(region)
             .filterDate(from_date, to_date)
            #  .select(['VH'])
             .select(['VV'])
             .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING')))

  # 日ごとの画像を生成
  daily_images = create_daily_images(from_date, to_date, Sentinel1, region)

  # 人工物（船舶）を抽出
  threshold = 0.5  # VV強度の閾値（人工物が検出される値）
  # 分割数（横方向）
  num_x = 8
  # 分割数（縦方向）
  num_y = 8
  # 解像度
  scale = 10

  for daily in daily_images:
      date = daily['date']
      image = daily['image']

      # 画像が空かどうかを確認
      if image.bandNames().size().getInfo() == 0:
          print(f"{date}: 空の画像をスキップしました。")
          continue

      # regionをタイルに分割
      tiles_grid = split_into_tiles(region, num_x=num_x, num_y=num_y)

      # regionマスク
      region_mask = ee.Image.constant(1).clip(region)

      # バイナリ画像（閾値より大きい部分）にも同じマスクを適用
      binary_image = image.gt(threshold)

      # 高解像度とバイナリ画像をタイルごとに保存する配列（行・列情報を考慮）
      high_res_grid = [[None for _ in range(num_x)] for _ in range(num_y)]
      binary_grid = [[None for _ in range(num_x)] for _ in range(num_y)]

      for row_tiles in tiles_grid:  # 各行ごとに処理
          for tile, row, col in row_tiles:  # 各タイルとその位置情報
              try:
                  img_array = geemap.ee_to_numpy(image, region=tile, scale=scale)
                  bin_array = geemap.ee_to_numpy(binary_image, region=tile, scale=scale)

                  img_array = np.squeeze(img_array)
                  bin_array = np.squeeze(bin_array)

                  high_res_grid[row][col] = img_array
                  binary_grid[row][col] = bin_array
              except Exception as e:
                  high_res_grid[row][col] = None
                  binary_grid[row][col] = None

      # 各行ごとにタイルを横方向に結合し、その結果を縦方向に結合して全体画像を作成
      if any(any(tile is not None for tile in row) for row in high_res_grid):
          stitched_high_res = np.vstack([
              np.hstack([tile for tile in row if tile is not None])
              for row in high_res_grid if any(tile is not None for tile in row)
          ])

          stitched_binary = np.vstack([
              np.hstack([tile for tile in row if tile is not None])
              for row in binary_grid if any(tile is not None for tile in row)
          ])

          # NumPy配列の前処理：NaNを0に変換し、データ型をuint8に変換
          binary_array = np.nan_to_num(stitched_binary).astype(np.uint8)

          # 連結成分解析の実施
          num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(binary_array)

          # 背景ラベルを除いた船舶数
          ship_count = num_labels - 1
          print(f"{date}: 検出された船舶数: {ship_count}")

          # カラーマップ設定（N/A を白にする）
          cmap_gray = plt.cm.gray.copy()
          cmap_gray.set_bad(color="white")  # N/A を白に

          cmap_binary = mcolors.ListedColormap(["black", "red"])  # 0: 黒, 1: 赤
          cmap_binary.set_bad(color="white")  # N/A を白に

          # 画像表示
          fig = plt.figure(figsize=(15, 15))
          gs = gridspec.GridSpec(3, 2, figure=fig)

          # 高解像度画像（全体）
          ax1 = fig.add_subplot(gs[0:2, 0])
          ax1.set_title("High-Res Image (Stitched)", fontsize=14)
          im1 = ax1.imshow(stitched_high_res, cmap=cmap_gray)
          ax1.axis("off")
          plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)

          # バイナリ画像（全体）
          ax2 = fig.add_subplot(gs[0:2, 1])
          ax2.set_title("Binary Image (Stitched)", fontsize=14)
          im2 = ax2.imshow(stitched_binary, cmap=cmap_binary)
          ax2.axis("off")
          plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04)

          # 検出された船舶の輪郭と重心を重ねる
          ax3 = fig.add_subplot(gs[2, :])
          ax3.set_title("Detected Components (Contours & Centroids)", fontsize=14)
          ax3.imshow(stitched_high_res, cmap=cmap_gray)
          ax3.axis("off")

          # バウンディングボックスや重心の描画
          cropped_images = []
          for i in range(1, num_labels):
              x, y, w, h, area = stats[i]
              cx, cy = centroids[i]

              # バウンディングボックスの拡張（周囲に余白を追加）
              margin = 40  # ピクセル単位の余白
              x1 = max(x - margin, 0)
              y1 = max(y - margin, 0)
              x2 = min(x + w + margin, stitched_high_res.shape[1])
              y2 = min(y + h + margin, stitched_high_res.shape[0])

              # 領域を切り出し
              cropped = stitched_high_res[y1:y2, x1:x2]

              # 切り出した領域を拡大（任意）
              # scale_factor = 2
              # resized = cv2.resize(cropped, None, fx=scale_factor, fy=scale_factor, interpolation=cv2.INTER_NEAREST)

              # リストに追加
              # cropped_images.append(resized)
              cropped_images.append(cropped)

          # 横方向に画像を並べて表示
          num_images = len(cropped_images)
          # 横方向の展示数
          cols = 4
          # 必要な行数を計算
          rows = math.ceil(num_images / cols)

          if num_images > 0:
              fig, axes = plt.subplots(rows, cols, figsize=(4*cols, 4*rows))
              axes = axes.flatten() if num_images > 1 else [axes]
              for idx, img in enumerate(cropped_images):
                  axes[idx].imshow(img, cmap=cmap_gray)
                  axes[idx].axis("off")
              for ax in axes[num_images:]:
                  ax.axis("off")
              plt.tight_layout()
              plt.show()
          else:
              print("表示する画像がありません。")

          ship_counts.append((date, ship_count))

解説

解析の中心となる関数です：
1. Sentinel-1 GRD コレクション（VV, 降り）を region と期間で絞り込み
  ※人工物はVVで反射が強い傾向があり、水面で弱くなる傾向があるためVVを使って反射が高いピクセルを船舶と判断しています。
2. create_daily_images で日別画像を生成
3. threshold = 0.5 で二値化（image.gt(threshold)）
  ※閾値の決定は正直恣意性が強いです。後述の結果を何度か確認し、一番適当と思われる値にしています。
  本来はgt0で良いと思われますが、0付近にノイズが多かったため0.5にしています。
4. split_into_tiles で領域を分割し、各タイルを geemap.ee_to_numpy(..., scale=10) で NumPy へ変換。
5. 行列でタイルを並べ直して表示（元画像・二値画像）
6. cv2.connectedComponentsWithStats により二値画像の閾値超過ピクセルを連結→船舶候補の個数をカウント
7. matplotlib で表示（元画像、二値画像、検出個所のハイライト）
各日付の検出数を ship_counts に (date, count) で蓄積し、連結して表示

大きな範囲を一括で処理しようとするとee_to_numpyのデータ容量制限に引っかかるため、一度タイルに分割してnumpy化した後に再度結合することによってエラーを回避しています。

解析の実行と時系列グラフ描画

ship_counts = []

# 2022年度を計算
ship_count(from_date= "2022-11-15", to_date= "2023-02-15")

# 2023年度を計算（能登半島地震発生）
ship_count(from_date= "2023-11-15", to_date= "2024-02-15")

# 2024年度を計算
ship_count(from_date= "2024-11-15", to_date= "2025-02-15")

date_list = [row[0] for row in ship_counts]
ship_list = [row[1] for row in ship_counts]

# 折れ線グラフを描画
plt.figure(figsize=(10, 5))
plt.plot(date_list, ship_list, marker='o')
plt.title('Shipping Count')
plt.xlabel('Date')
plt.ylabel('Ship Count')
plt.xticks(rotation=45)
plt.grid(True)
plt.tight_layout()
plt.show()

解説

地震発生年度前後1年ずつ、計3つの期間を対象に ship_count を実行：
- 2022-11-15〜2023-02-15
- 2023-11-15〜2024-02-15（能登半島地震期を含む）
- 2024-11-15〜2025-02-15
ship_counts から日付と個数を取り出し、折れ線グラフで推移を表示

結果
おそらく外れ値を含んでいるであろう2022-12-18の結果を貼り付けます。
①元画像(タイル連結後の高解像度画像)

雲のようなノイズが見て取れます。よく目を凝らせばその下に所々光っている箇所があります。
この日のRGBデータを確認してみましたが、SARとRGBで撮影時間が同じものがなかったので推測になりますが、RGBでは快晴のように見えたのでこれは雲ではなく波等のノイズだと思っています。
②二値画像

この日に限った話ではないですが、本来は閾値超過のピクセルが赤い点で表示されますが、範囲が広すぎるので潰れて見えなくなっています。対象や範囲を変えれば閾値超過ピクセルがよく見えるかと思います。
③ヒストグラム

どの強度の反射地がどれくらいあるかをプロットしています。このヒストグラムを見て分かる通り、閾値を0.5にするのが正しいのかは懸念が残ります。他の日にちのヒストグラムを見ていると大体、-0.5〜+0.5あたりの反射値が外れ値的にカウント数が抜きん出て多かったため、そのあたりがノイズなのだろうと判断しています。
④ハイライト

ハイライトはCV2によって連結された箇所を表示しています。つまり、本プログラムで船舶と判断された画像となります。ただし、Component1,2を見ると分かる通り、おそらく同じ対象をダブルカウントしていたりします。目で見ると分かりづらいですが、ピクセル的に見ると途切れているために1つの連携対象にならなかったのだと思います。連結範囲を広げるためにfocalMean(radius=1,kernelType='square')を行っていますが、やはりやりすぎると本来2つ以上に分けるべき対象も1つになってしまう為難しいところです。一応正方形で半径1を対象に平滑化しており、解像度的に1マス10mなので船の大きさ的には上記のパラメータが適当だと考えました。
⑤トレンドグラフ

2022/11/15-2025/02/15までを飛び石でみてみましたが、2024/1/1付近での大きな差異は見て取れず、移動平均をとってもそこまで有意な差は無いかと思います。

まとめ

悩まされた課題は主に以下でした。

ノイズの除去
閾値決定
ピクセル連結処理
上記解決のため試せそうな方法があればコメント等でご教示いただけると幸いです。
災害日時前後の船舶数トレンドを押さえればその付近の港の回復具合がわかるかもしれないと思って確認してみましたが、まだまだ精度に対しての課題があり、今回の結果だけを見て「災害時の港復興具合は船舶トレンドでは追えない」、もしくは「能登半島地震では地震前後で港を出ている船舶数に変化がないため出港に関しての被災程度は低い」と結論付けることはできないかと思いました。
そのため、今回の記事は「大きな範囲のSARデータを小分けにして処理、表示することができるようになるまで」として参照いただければ幸いです。