本記事は、株式会社Tellusにインターンとして参加した大阪成蹊大学　データサイエンス学部3年の伊藤宗真さんが、初めて衛星データの可視化に取り組んだ奮闘の過程をご紹介する記事です。衛星データにあまり触れたことがない方の参考になれば幸いです。

今回衛星データを使ってチャレンジしたのは、真珠養殖に適した新しい海域を探すための、海の状態の可視化です。

私は真珠の養殖・加工を生業とする家庭に生まれ、真珠とともに生活してきました。

ただ、近年、真珠の生産量が減ってきており、問題解決につながる様々な取り組みが行われています。

そこで、これまで衛星データやプログラミングとは無縁の生活を送ってきましたが、Tellusにインターンとして参加したことをきっかけに、真珠の生産量をあげるための新たな可能性を模索するために、衛星データの可視化に挑戦してみることにしました。

まず海面水温に注目し、GCOM-C衛星のデータを使って英虞湾周辺の水温の変化を調べることにしました。

問題背景：真珠養殖が抱える問題と現時点の解決策

三重県の英虞湾や愛媛県の宇和島湾は、真珠を育てる母貝であるアコヤ貝の養殖が古くから盛んな地域として知られています。しかし近年、真珠の生産量は減少傾向にあり、真珠の養殖業者や真珠を愛する人々にとって深刻な課題となっています。
https://www.jfa.maff.go.jp/j/saibai/attach/pdf/shinju-2307-1.pdf

こうした状況を引き起こしている要因としては、環境的・社会的なさまざまな問題が挙げられます。例えば、赤潮などの異常現象、稚貝の病気の発生、気候変動に伴う海水温の上昇や海洋環境の変化、高齢化による人手不足などがその代表です。

そのような課題に対処するための解決例をひとつ挙げると、現在、真珠養殖業者や研究機関によって、水質管理の徹底や病気に強いアコヤ貝の開発を目指したゲノム研究が進められています。

こうした中で、アコヤ貝の生産地を新たに増やすことによって、真珠の生産量全体を底上げできるのではないかと考えました。

課題設定：衛星データを利用した海洋状況の可視化

従来の養殖海域だけで対応するのではなく、新たに環境条件の整った海域を見つけ出すことが、将来的な持続可能な真珠養殖に向けた鍵となるはずです。

しかし、新たな生産地を見つけるには、まず「どの海域がアコヤ貝の養殖に適しているか」を科学的に判断する必要があります。水温や水質といった海況は日々変化しており、それらを現場レベルで網羅的に観測することは現実的には困難です。

そこで注目したのが、衛星リモートセンシングを用いた海況の可視化というアプローチです。衛星から得られるデータを使えば、広範囲かつ長期的な視点で海の状態を把握することが可能になると考えています。

本記事では、その第一歩として、海面水温（SST）に着目し、衛星データを用いて英虞湾周辺の水温変動を時系列で解析することに取り組みました。海面水温は真珠養殖において最も基本的かつ重要な指標の一つであり、まずこのデータを正確に把握することで、将来的な適地選定の基盤を築くことができます。

最終的に水温以外にも、塩分濃度、クロロフィルa濃度、溶存酸素量といった他の海洋環境要因も統合的に分析し、真珠養殖適地選定モデルの構築へとつなげていくことを目標としています。

実施内容：GCOM-Cの時系列海面水温データを抽出・可視化

先ほども述べたように本記事では、海面水温データに注目し、これを英虞湾周辺の時系列データとして抽出・可視化する処理を構築しました。

Pythonを用いたデータ抽出・補間処理のプログラムを作成し、結果をグラフで可視化することで、海況の変化を直感的に捉えることを試みました。

使用する衛星「GCOM-C」とは

今回使用する衛星は GCOM-C（しきさい） です。
GCOM-C は、JAXA（宇宙航空研究開発機構）が打ち上げた地球観測衛星で、正式名称は気候変動観測衛星「しきさい」（GCOM-C：Global Change Observation Mission- Climate）です。主に地球表面や大気、海洋の状態を高頻度かつ広範囲に観測することを目的としています。

この衛星には「SGLI（Second-generation Global Imager）」という光学センサーが搭載されており、可視・近赤外・熱赤外の複数の波長帯で地球を撮影することができます。そのデータからは、海面水温（SST）、植生の変化、大気中のエアロゾルや雲の情報など、多様な地球環境データを取得できます。
https://www.satnavi.jaxa.jp/ja/project/gcom-c/index.html

【関連記事】

データの取得について

では実際に作成したコードを説明していこうと思います。

今回使用したGCOM-CデータはG-portalから入手できます。GCOM-CではLEVEL1、 LEVEL2、 LEVEL3と処理レベル別にデータを選択できます。
　
簡単に処理レベルごとに説明します。GCOM-Cの観測データは、処理の進み具合に応じて複数のレベルに分類されており、それぞれのデータは異なる用途に適しています。

まず、Level-1Aは、重複パケットの削除やシーンの切り出し、放射補正用情報の付加、幾何情報の算出など、観測データに対して初歩的な整備を行ったものであり、主に再処理や補正の前段階で使用されます。

次に、Level-1Bは、Level-1Aをもとに放射補正や幾何補正を実施したデータで、観測値としての信頼性が高められた製品です。波長ごとの補正情報も付加されており、解析に直接利用可能なレベルの品質となっています。

Level-2では、Level-1Bなどの観測データから、海面水温（SST）やクロロフィルa濃度といった具体的な物理量が計算され、ピクセル単位での科学的分析が可能になります。これらのデータは、250m〜4kmグリッドで空間的に整備されており、全球スケールでの比較や傾向把握に適しています。

さらに、Level-3では、Level-2データを時間的・空間的に合成し、8日平均や月平均といった統計処理を加えた製品が提供されます。これにより、日々の観測のばらつきをならし、広域的・長期的な環境変化の把握が容易になります。

こちらの内容はJAXA公式サイト内のGCOM-Cデータに関するプロダクト説明ページに記載されています。こちらのサイトにはGCOM-Cによって提供される各データプロダクトの処理レベルや内容、特性について、詳しく説明されています。
URL：https://suzaku.eorc.jaxa.jp/GCOM_C/data/product_def_j.html

今回は沿岸部のデータを分析するため解像度の高いLEVEL2の250ｍ分解能のデータを選択しました。

次にデータファイルの取得方法を説明します。GCOM-CのデータはJAXAの公式ポータルサイト「G-Portal」（https://gportal.jaxa.jp/gpr/?lang=ja）から取得することができます。

無料のアカウントを作成し、「衛星からの検索」→「GCOM-C/SGLI」のプルダウンを開きます。前述の通り今回はlevel2のデータを使用するので、さらに「LEVEL2」→「海洋圏」のプルダウンを開きます。そうすると「L2-海面水温」がありますのでそちらにチェックを入れてください。

次に詳細な設定を行いたいのでチェックボックス右側にある歯車のアイコンをクリックしてください。下記写真のように品質を「Good」、解像度を「250m」に設定すると完了です。

後は取得したい期間（今回は雲が比較的少なく感じる2025年2月1日から28日に指定）を期間指定の欄から選択し、範囲指定の欄から短形指定を選択、地図上で取得したい地点（今回は英虞湾周辺）を覆うように◻︎をざっくり指定し検索を押すと、条件にあったファイルの一覧が表示されるので、1つずつダウンロードしていってください。

また、出力されるデータには2種類の軌道方向（昇交・降交）の情報が含まれているため、同じ日付であっても2枚のファイルが生成されます。今回は、指定地点よりも広い範囲でデータを抽出したため、1日あたり3〜4枚のファイルが出力されている場合がありますが、指定地点のファイル以外はエラーファイルとして扱われるため、特に問題はありません。

以下に取得したデータファイルを羅列しておきます。

GC1SG1_202502010145V05810_L2SG_SSTDQ_3001.h5
GC1SG1_202502011258V29602_L2SG_SSTNQ_3001.h5
GC1SG1_202502011303D29603_L2SG_SSTNQ_3001.h5
GC1SG1_202502020119E04910_L2SG_SSTDQ_3001.h5
GC1SG1_202502020123K04911_L2SG_SSTDQ_3001.h5
GC1SG1_202502021236J28703_L2SG_SSTNQ_3001.h5
GC1SG1_202502030052L04010_L2SG_SSTDQ_3001.h5
GC1SG1_202502030056Q04011_L2SG_SSTDQ_3001.h5
GC1SG1_202502031209Q27803_L2SG_SSTNQ_3001.h5
GC1SG1_202502040206Q06510_L2SG_SSTDQ_3001.h5
GC1SG1_202502041319Q30302_L2SG_SSTNQ_3001.h5
GC1SG1_202502041323U30303_L2SG_SSTNQ_3001.h5
GC1SG1_202502050139V05610_L2SG_SSTDQ_3001.h5
GC1SG1_202502191322H30402_L2SG_SSTNQ_3001.h5
GC1SG1_202502191326N30403_L2SG_SSTNQ_3001.h5
GC1SG1_202502200142P05710_L2SG_SSTDQ_3001.h5
GC1SG1_202502201255P29502_L2SG_SSTNQ_3001.h5
GC1SG1_202502201259T29503_L2SG_SSTNQ_3001.h5
GC1SG1_202502210115U04810_L2SG_SSTDQ_3001.h5
GC1SG1_202502210120C04811_L2SG_SSTDQ_3001.h5
GC1SG1_202502211233C28603_L2SG_SSTNQ_3001.h5
GC1SG1_202502220049D03910_L2SG_SSTDQ_3001.h5
GC1SG1_202502220053H03911_L2SG_SSTDQ_3001.h5
GC1SG1_202502221206H27703_L2SG_SSTNQ_3001.h5
GC1SG1_202502230203H06410_L2SG_SSTDQ_3001.h5
GC1SG1_202502231316H30202_L2SG_SSTNQ_3001.h5
GC1SG1_202502231320M30203_L2SG_SSTNQ_3001.h5
GC1SG1_202502240136N05510_L2SG_SSTDQ_3001.h5
GC1SG1_202502241249N29302_L2SG_SSTNQ_3001.h5
GC1SG1_202502241253S29303_L2SG_SSTNQ_3001.h5
GC1SG1_202502250109T04610_L2SG_SSTDQ_3001.h5
GC1SG1_202502250114B04611_L2SG_SSTDQ_3001.h5
GC1SG1_202502251227B28403_L2SG_SSTNQ_3001.h5
GC1SG1_202502260043C03710_L2SG_SSTDQ_3001.h5
GC1SG1_202502260047G03711_L2SG_SSTDQ_3001.h5
GC1SG1_202502260224B07110_L2SG_SSTDQ_3001.h5
GC1SG1_202502261200G27503_L2SG_SSTNQ_3001.h5
GC1SG1_202502270157G06210_L2SG_SSTDQ_3001.h5
GC1SG1_202502271310F30002_L2SG_SSTNQ_3001.h5
GC1SG1_202502271314K30003_L2SG_SSTNQ_3001.h5
GC1SG1_202502280130M05210_L2SG_SSTDQ_3001.h5
GC1SG1_202502281243L29002_L2SG_SSTNQ_3001.h5
GC1SG1_202502281247Q29003_L2SG_SSTNQ_3001.h5
GC1SG1_202502051252V29402_L2SG_SSTNQ_3001.h5
GC1SG1_202502051257D29403_L2SG_SSTNQ_3001.h5
GC1SG1_202502060113F04710_L2SG_SSTDQ_3001.h5
GC1SG1_202502060117K04711_L2SG_SSTDQ_3001.h5
GC1SG1_202502061230J28503_L2SG_SSTNQ_3001.h5
GC1SG1_202502070046L03810_L2SG_SSTDQ_3001.h5
GC1SG1_202502070050Q03811_L2SG_SSTDQ_3001.h5
GC1SG1_202502070227K07210_L2SG_SSTDQ_3001.h5
GC1SG1_202502071203Q27603_L2SG_SSTNQ_3001.h5
GC1SG1_202502080200Q06310_L2SG_SSTDQ_3001.h5
GC1SG1_202502090133V05410_L2SG_SSTDQ_3001.h5
GC1SG1_202502091246V29202_L2SG_SSTNQ_3001.h5
GC1SG1_202502091251D29203_L2SG_SSTNQ_3001.h5
GC1SG1_202502100107F04510_L2SG_SSTDQ_3001.h5
GC1SG1_202502100111K04511_L2SG_SSTDQ_3001.h5
GC1SG1_202502101224J28303_L2SG_SSTNQ_3001.h5
GC1SG1_202502110040L03610_L2SG_SSTDQ_3001.h5
GC1SG1_202502110044Q03611_L2SG_SSTDQ_3001.h5
GC1SG1_202502110221K07010_L2SG_SSTDQ_3001.h5
GC1SG1_202502111157Q27403_L2SG_SSTNQ_3001.h5
GC1SG1_202502120154Q06110_L2SG_SSTDQ_3001.h5
GC1SG1_202502121307Q29902_L2SG_SSTNQ_3001.h5
GC1SG1_202502121311U29903_L2SG_SSTNQ_3001.h5
GC1SG1_202502130127V05210_L2SG_SSTDQ_3001.h5
GC1SG1_202502131245D29003_L2SG_SSTNQ_3001.h5
GC1SG1_202502140101E04310_L2SG_SSTDQ_3001.h5
GC1SG1_202502140105K04311_L2SG_SSTDQ_3001.h5
GC1SG1_202502141218J28103_L2SG_SSTNQ_3001.h5
GC1SG1_202502150215K06810_L2SG_SSTDQ_3001.h5
GC1SG1_202502151151P27203_L2SG_SSTNQ_3001.h5
GC1SG1_202502151328J30602_L2SG_SSTNQ_3001.h5
GC1SG1_202502151332N30603_L2SG_SSTNQ_3001.h5
GC1SG1_202502160148Q05910_L2SG_SSTDQ_3001.h5
GC1SG1_202502161301P29702_L2SG_SSTNQ_3001.h5
GC1SG1_202502161305T29703_L2SG_SSTNQ_3001.h5
GC1SG1_202502170121V05010_L2SG_SSTDQ_3001.h5
GC1SG1_202502170126D05011_L2SG_SSTDQ_3001.h5
GC1SG1_202502171239C28803_L2SG_SSTNQ_3001.h5
GC1SG1_202502180055E04110_L2SG_SSTDQ_3001.h5
GC1SG1_202502180059J04111_L2SG_SSTDQ_3001.h5
GC1SG1_202502181212J27903_L2SG_SSTNQ_3001.h5
GC1SG1_202502190209J06610_L2SG_SSTDQ_3001.h5

作成したコードについて

このコードの全体の処理の流れは以下のようになります。

1．初期設定
2．補完処理
3．補完後データの検証
4．対象地点に近いピクセルの特定
5．データの異常値の特定
6．sstに保存されているデータを物理量に変換
7．1枚のHDFファイルから必要な情報の抽出・記録
8．時系列グラフのプロット
9．main処理

では説明していきます。

【1】初期設定：ライブラリのインポートとユーザ設定

# 使用するライブラリをインポート
import h5py
import pandas as pd
import matplotlib.pyplot as plt
import glob
import re
import numpy as np
import cv2
from collections import Counter
from scipy import stats
# ── ユーザー設定 ───────────────────────────────────
TARGET_LAT   = 34.2
TARGET_LON   = 136.7
DATA_FOLDER  = "your_folder"
SLOPE        = 0.0012
OFFSET       = -10
# ────────────────────────────────────────────

まず最初に初期設定として、今回の処理で必要となる各種ライブラリをインポートしています。h5pyはHDF5形式の衛星データを読み込むために、pandasやnumpyはデータ処理のために使用されます。また、matplotlib.pyplotは可視化、cv2（OpenCV）は画像の補間処理、scipy.statsやcollections.Counterは統計処理に利用されます。

その後、ユーザー設定として、対象とする緯度経度、読み込むデータのフォルダパス、およびSST（海面水温）のスケーリングに使用する係数を指定しています。これらの値は後続の処理において、データの読み込み・変換・分析の基準として用いられます。

【2】interpolate_with_cv2()：衛星データを滑らかに補間する

def interpolate_with_cv2(coarse:np.ndarray, interval:int)-> np.ndarray:
    array = np.array(coarse, dtype=np.float32)
    rows, cols = array.shape
    new_size = ((cols - 1) * interval + 1, (rows - 1) * interval + 1)
    resized = cv2.resize(array, new_size, interpolation=cv2.INTER_LINEAR)
    return resized

この関数は、衛星から得られた粗い格子状のデータ（緯度や経度など）を滑らかに補間するためのものです。今回使用した衛星データでは、緯度や経度が空間的に滑らかに変化するため、すべての画素に対して高解像度で保持する必要がなく、データ容量を削減し効率的に処理できるように設計されています。そのため、ある程度「間引きされた」低解像度の座標データが用いられています。

しかし、私たちが特定の緯度・経度を正確に扱いたい場合には、より細かく補間された座標情報が必要です。

この補間には、画像処理ライブラリである OpenCV の cv2.resize() 関数を使います。これは画像を拡大・縮小する際に使う関数ですが、数値の格子データにも応用できます。

線形補間（INTER_LINEAR）によって、もとのデータと連続性を保ちながら新しい細かい格子を生成します。たとえば、下記の図のようにもとの緯度データが3×3のグリッドであっても、補間によって7×7のような細かな緯度マップを得ることができます。

【3】verify_interpolation_bounds()：補間後のデータが正しいか検証する

def verify_interpolation_bounds(coarse_array, fine_array, interval, label=""):
    rows, cols = coarse_array.shape
    fine_rows, fine_cols = fine_array.shape
    # 想定される位置
    i_fine_max = (rows - 1) * interval
    j_fine_max = (cols - 1) * interval
    corners = [
        ("左上",    (0, 0),                   (0, 0)),
        ("右上",    (0, cols - 1),            (0, j_fine_max)),
        ("左下",    (rows - 1, 0),            (i_fine_max, 0)),
        ("右下",    (rows - 1, cols - 1),     (i_fine_max, j_fine_max))
    ]
    print(f"\n=== {label} 補間境界の検証 ===")
    for name, (i_coarse, j_coarse), (i_fine, j_fine) in corners:
        v_coarse = coarse_array[i_coarse, j_coarse]
        v_fine   = fine_array[i_fine, j_fine]
        ok = np.isclose(v_coarse, v_fine, atol=1e-4)
        #print(f"{name}: coarse = {v_coarse:.6f}, fine = {v_fine:.6f} → {'OK' if ok else 'NG'}")

補間処理を行う際には、「データが勝手に書き換わってしまっていないか」という検証も重要です。この関数は、補間後のデータの四隅が、補間前の値と一致しているかをチェックするために使われます。

具体的には、補間前の左上・右上・左下・右下の値を記録し、補間後の同じ位置に対応する値と比較します。誤差があまりに大きければ、補間処理に問題がある可能性があります。

今回は検証用のため、出力はコメントアウトされていますが、開発段階でのデバッグに役立つ関数です。

【4】interpolate_and_verify()：対象地点に最も近いピクセルを特定する

def interpolate_and_verify(lat_coarse, lon_coarse, lat_int, lon_int, sst_dn)-> tuple[int, int]:
    # 緯度・経度間隔の補完
    full_lat = interpolate_with_cv2(lat_coarse, lat_int)
    full_lon = interpolate_with_cv2(lon_coarse, lon_int)
    rows_full, cols_full = sst_dn.shape
    full_lat = full_lat[:rows_full, :cols_full]
    full_lon = full_lon[:rows_full, :cols_full]
    dists = (full_lat - TARGET_LAT)**2 + (full_lon - TARGET_LON)**2
    i_min, j_min = np.unravel_index(np.argmin(dists), dists.shape)
    print(f"  補間位置: ({i_min}, {j_min})")
    # # 検証
    # verify_interpolation_bounds(lat_coarse, full_lat, lat_int, label="緯度")
    # verify_interpolation_bounds(lon_coarse, full_lon, lon_int, label="経度")
    return i_min, j_min

この関数では、interpolate_with_cv2()で補間した緯度・経度のデータを使って、対象地点（今回は緯度34.2度、経度136.7度）に最も近いピクセル位置を探し出します。

まず、補間された緯度・経度のマップをそれぞれ取得し、それぞれのピクセルでの距離（ユークリッド距離の2乗）を計算します。全体の中から、その距離が最も小さいピクセル、つまり「観測地点に最も近い地点」を選び出します。

このようにして、衛星データの中から「自分が知りたい位置の情報」を特定できるようにしています。

【5】analyze_error()：データに含まれる異常値の種類を特定する

def analyze_error(block: np.ndarray) -> str:
    error_types = {
        65535: "error_DN",
        65533: "Cloud_mask_DN",
        65532: "Retrieval_error_DN",
        0: "Missing_data_DN"
    }
    flat = block[~np.isnan(block)].astype(int).flatten()
    mode_val, _ = stats.mode(flat, keepdims=False)
    label = error_types.get(int(mode_val), f"Valid_or_Unknown : {mode_val}")
    return label

衛星データには、天候条件やセンサーの不具合などによって、雲の影・欠損・エラー値などの異常データが含まれていることがあります。

この関数は、ある小さな領域（10×10ピクセル程度）のデータブロックを調べて、最も頻出している異常コードを分析するものです。たとえば「65535」はセンサーエラー、「65532」は取得エラー、「0」はデータなし、といったように、あらかじめ定義されたコードに応じてラベルをつけます。エラーについては下記URLのAttribute情報をご確認ください。

https://suzaku.eorc.jaxa.jp/GCOM_C/data/update/Algorithm_SST_ja.html

出力結果のデータが少なかった時にどういったエラーが出ているか簡単に可視化する、検証用の関数です。

【6】get_dn()：SST（海面水温）のDN値を実際の温度に変換する。

def get_dn(sst_dnblock: np.ndarray) -> tuple[float, str]:
    # SST変換
    if sst_dnblock.size == 0:
        print("block is empty")
        return None, "out of range"
    # dn値が65500以上の値は欠損値として扱う(参照URL: https://suzaku.eorc.jaxa.jp/GCOM_C/data/update/Algorithm_SST_ja.html)
    sst_dnblock_with_nan = np.where(sst_dnblock < 65500, sst_dnblock, np.nan)
    try:
        dn_val = int(np.nanmedian(sst_dnblock_with_nan))
        sst_val = dn_val * SLOPE + OFFSET  # DN値から実際の温度に変換(参考URL: https://suzaku.eorc.jaxa.jp/GCOM_C/data/update/Algorithm_SST_ja.html)
        print(f"  SST: DN={dn_val}, 実値={sst_val:.2f}°C")
        if not (-2 <= sst_val <= 40):
            print("  → SST out of realistic range, skip")
            return None, "error_SST"
        else:
            return sst_val, f"Valid_SST : {sst_val:.2f}°C"
    except ValueError as e:
        message = str(e)
        if "cannot convert float NaN to integer" in message:
            # print("  → NaN detected in SST block, analyzing error")
            return None, analyze_error(sst_dnblock)

この関数は、衛星が取得したDN値（Digital Number：センサーの出力）を物理量（摂氏温度）に変換する役割を持っています。

まず、65500以上のDN値は欠損値として扱うためにNaNに置き換えます。

次に、対象地点周辺の領域（【6】get_sst_data()関数で切り出したブロック）から得られたDN値の中からNaNを除外し、その中央値（代表値）を算出します。中央値を用いることで、外れ値やノイズの影響を抑え、安定した代表値を得ることができます。

この中央値に、JAXAが公開するプロダクト定義書に記載された係数（SLOPE = 0.0012、OFFSET = -10）を用いて「SST = DN × SLOPE + OFFSET」の式で物理量に変換を行い、実際の海面水温（°C）へと換算します。計算式と係数は下記URLに記載されています。
https://suzaku.eorc.jaxa.jp/GCOM_C/data/update/Algorithm_SST_ja.html

もし変換後の値が−2℃以下または40℃以上といった非現実的な範囲にある場合は、無効としてラベル付けされます。欠損が多くて中央値が取れないときは、analyze_error()に渡して異常の種類を調べます。

【7】get_sst_data()：1枚のHDF5ファイルから必要な情報を抽出・記録

def get_sst_data(fp: str):
    with h5py.File(fp, "r") as f:
        # h5pyファイルから必要なデータを取得
        m = re.search(r"GC1SG1_(\d{8})", fp)
        if not m:
            print("  → Date parse failed, skip")
            return None, None
        date = pd.to_datetime(m.group(1), format="%Y%m%d")
        # Geometryデータからそれぞれ取得
        lat_coarse = f["Geometry_data/Latitude"][:]  # 緯度の補間前データ
        lon_coarse = f["Geometry_data/Longitude"][:] # 経度の補間前データ
        lat_int = int(f["Geometry_data/Latitude"].attrs["Resampling_interval"][0])  # 緯度の補間間隔
        lon_int = int(f["Geometry_data/Longitude"].attrs["Resampling_interval"][0])  # 経度の補間間隔
        # obs_time_coarse = f["Geometry_data/Obs_time"][:]  # 観測時刻の補間前データ
        dt_from_fp = re.search(r"GC1SG1_(\d{12})", fp)
        obs_time_coarse = f["Geometry_data/Obs_time"][:]
        obs_dt = pd.to_datetime(dt_from_fp.group(1), format="%Y%m%d%H%M")
        sst_dn = f["Image_data/SST"][:]  # SSTのDN値
        cloud_prob = f["Image_data/Cloud_probability"][:]  # 雲確率
        i_min, j_min = interpolate_and_verify(lat_coarse, lon_coarse, lat_int, lon_int, sst_dn)
        kernel = 10
        sst_dnblock = sst_dn[i_min - kernel : i_min + kernel, j_min - kernel : j_min + kernel]
        sst_val, label = get_dn(sst_dnblock)
        # 軌道の取得
        orbit = "Descending" if 6 <= obs_dt.hour <= 18 else "Ascending"
        record = {
            "datetime": obs_dt,
            "SST": sst_val,
            "CloudProb": float(cloud_prob[i_min, j_min]),
            "Orbit": orbit,
            "label": label
        }
        return record

この関数は、1枚の衛星データファイル（.h5）から、指定地点のSSTや被雲率、観測時刻などの情報を抽出して記録する役割を担っています。
この関数で行われている手順としては以下の通りです：

1.ファイル名から観測日時の抽出
HDF5ファイルのファイル名には、観測日時が含まれています。この情報を正規表現で抽出し、pd.to_datetime() を使って日付・時刻形式に変換します。

2.緯度・経度の元データを読み取り、補間間隔を取得
衛星は、地球の表面を大まかに格子状に観測しています。この格子情報である緯度・経度の元データ（低解像度）を読み取ります。さらに、データに付属する属性情報（attribute）から、補間の倍率（何倍に細かく補完すべきか）も取得します。

3.SSTのDN値と雲確率データを読み取り
ここでは、SST（海面水温）の観測値（DN値）と、雲の被り具合（Cloud Probability）を取得します。

4.対象地点のピクセルを特定し、その周辺のデータブロックの抽出
対象地点（例：34.2°, 136.7°）に最も近いピクセルを補間によって探し、データの解像度をあげるためにその周囲を含む20×20ピクセルのブロックを抽出します。
【3】interpolate_and_verify() 関数を通じて、粗い座標データを補間し、目的地点に最も近い格子点の位置（行・列）を探します。その周囲（上下左右10ピクセルずつ）を切り出して、局所的な平均的な状況を把握します。このブロックが【5】get_dn()関数を用いた温度変換処理に使用されます。

5.軌道（昇り/降り）も判別してラベル付け
観測時刻をもとに、衛星が**昇交軌道（夕〜夜／南→北）か降交軌道（朝〜昼／北→南）のどちらを通過しているかを判定しています。
GCOM-Cは太陽同期準回帰軌道を採用しており、毎回ほぼ同じ地方時に同じ地域を通過します。
その特性を利用し、観測時刻が6〜18時（朝〜昼）であれば、衛星は地球を北から南に向かって通過している（降交軌道）**と判定します。
それ以外の時間帯（18〜翌6時）は、**南から北に向かって通過している（昇交軌道）**とみなしています。

下図はアセンディングとディセンディングをわかりやすく表した図になっております。上記で説明したように同じ地点でも軌道によって撮影向きとその地点を通過する時間は違います。

6.最終的に、辞書形式の「1日分の観測記録」として返す
抽出した観測情報を辞書形式（Pythonのデータ構造）にまとめて返します。この情報が、のちのグラフ描画や統計分析に利用されます。

【8】plot_sst_data()：海面水温の変化を時系列グラフで表示する

def plot_sst_data(records):
    df = pd.DataFrame([d for d in records if "Valid_SST" in d["label"]])
    if df.empty:
        print("有効なデータが取得できませんでした。")
        quit()
    df = df.sort_values("datetime")
    start = pd.to_datetime("2025-02-01")
    end = pd.to_datetime("2025-02-28")
    df_month = df[(df["datetime"] >= start) & (df["datetime"] <= end)]
    if df_month.empty:
        print("指定期間内に有効データがありません。")
        quit()
    plt.figure(figsize=(12, 6))
    for label, group in df_month.groupby("Orbit"):
        plt.plot(group["datetime"], group["SST"], marker="o",linestyle="--" if label == "Descending" else "-", label=f"{label} Pass")
    plt.xlabel("DateTime")
    plt.ylabel("Sea Surface Temperature (°C)")
    plt.title(f"GCOM-C SST at ({TARGET_LAT}, {TARGET_LON})")
    plt.minorticks_on()
    # datelabels = [date.strftime("%Y-%m-%d") for date in pd.date_range(start, end, freq="D")]
    # plt.set_xticklabels(datelabels, rotation=45)
    plt.xticks(rotation=45)
    plt.grid(which="major", linestyle="-", linewidth=0.8)
    plt.grid(which="minor", linestyle=":", linewidth=0.5)
    plt.legend()
    plt.tight_layout()
    plt.show()

この関数は、複数のファイルから取得した観測記録をもとに、2月1日から2月28日までの海面水温の変化を時系列でプロットします。データを「降交軌道（朝〜昼）」と「昇交軌道（夕方〜夜）」に分け、色や線種を変えて視覚的に違いがわかるようにしています。これにより、昼と夜で水温がどのように異なるかを一目で理解できるようにしています。

【9】main()：すべての処理を制御する中心の関数

最後の main() 関数は、プログラム全体の処理の流れを統括する役割を担っています。
まず、指定されたフォルダ内からすべての .h5 ファイルを検索・取得し、それぞれのファイルに対して 【6】get_sst_data() 関数を実行して、観測結果を記録として保存します。その後、取得したデータをまとめて 【7】plot_sst_data() 関数に渡し、海面水温の推移をグラフとして可視化します。

以上が、今回作成したコードにおける一連の処理の流れとなります。
なお、本コードを実行し、生成された可視化結果を以下に添付しております。

2025年2月の時系列データの結果から、対象海域の海面水温はおおむね15℃前後で推移していることが分かりました。

この水温は、アコヤ貝の活動がやや鈍るとされる15℃付近にあたります。養殖においては、アコヤ貝が活発に摂餌・成長するためには15〜25℃程度が理想とされているため、2月の水温環境は理想値から下回ってしまっていることも踏まえ決して良いものとは言い切れません。これは、想定していたよりも水温変動の幅が大きく、あくまで海面水温ですが水温管理の見直しが必要となる可能性を示唆しています。

一方で、データからは日中（降交軌道）と夜間（昇交軌道）で水温に大きな差がないことも確認でき、短期間での急激な変動は発生していないことから、ある程度安定した環境ではあるとも考えられます。

躓いた点と解決方法

今回のコード作成は初めての本格的なコード作成だったこともあり、自分の中で“調べて組めば「ある程度動くだろう」という楽観的な見込みからスタートしました。

しかし、いざ手を動かし始めると、想定以上に技術的・認識的な壁にぶつかる連続で、ほとんどが「何がわからないのかもわからない」状態でした。

衛星データの補間処理

最初に大きくつまずいたのは、衛星データを扱う際の「補間処理」が前提として必要であることにすら気づいていなかったことです。

GCOM-Cのプロダクトは格子状の観測点に沿ってデータを持っているため、観測点と観測点の間の値を推定する必要がありますが、当初、座標を指定すれば自動的にその場所の海面水温が得られるものだと誤解していました。結果として、指定した座標とは異なる位置のピクセルが「最近傍」として出力され、地図上で明らかにずれた値が出力されていました。

この段階ではまだその異常に気づいておらず、何枚も“誤ったマップ”を生成してしまう状態が続きました。その後、プロダクト仕様書（下記URL：p120）を改めて確認し、「補間処理が必要だ」という非常に基本的な事実にやっと気づきました。
https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/GCOM-C_Products_Users_Guide_entrylevel_jp.pdf

また、Pythonで解析することを想定していましたが、　プロダクト仕様書（p120）には MATLAB の interp2 を使った例しか掲載されておらず、Pythonで代替できる方法を模索する必要が出てきました。

Pythonでは scipy.interpolate.interp2d、numpy.interp、cv2.resize など複数の方法があるものの、それぞれ引数の扱いや補間方法が微妙に異なり、どれが「衛星データの補間」に適しているのか判断が困難でした。

OpenCVの cv2.resize で試してみると補間は滑らかに見えましたが、「ただピクセルを拡大してるだけではないか？本当に意味のある補間か？」と疑念が湧き、何度も出力画像と座標の関係を検証し直す必要がありました。

最終的には cv2.resize を使って線形補間を行う方法で落ち着きましたが、補間の仕組み自体を自分で理解し、納得できるまでに数日かかりました。

今回試した方法を表形式にまとめてみました。

データの欠損の多さ

さらに追い打ちをかけたのが、出力された時系列データにエラーDN値（655○○）が大量に含まれていた問題です。

最初に時系列グラフを描いたとき、数ヶ月分のデータを処理しても、グラフ上には点が数個しか出てこなかったです。以下にグラフも載せておきます。

その原因を調べるにあたって、エラーが出ている箇所をすべてプリントデバッグで一つずつ確認し、どのタイミングでどのような値が出力されているのかを丁寧に追っていきました。その結果、特定のタイミングで655○○という値が繰り返し出力されていることに気づき、GCOM-Cのプロダクト仕様を改めて確認したところ、この値は雲の影響などで観測が正常に行えなかった場合に記録されるエラーDN値であることが判明しました。

下部にプロダクト説明書のスクリーンショットと実際にエラーが出たファイルのsstをカラーマップとして出力したものです。このカラーマップの中では白色の部分がエラーDN値です。

こういった結果から、対象地点の1ピクセルだけに依存するのではなく、周囲の複数ピクセルを取得し、その中央値を取る方法を考案しました。作成したコードの【5】get_dn()の部分です。この方法により、たとえあるピクセルにエラー値が含まれていても、他のピクセルから有効な値を補間できる可能性が高まり、結果的に有効なデータ数を増やすことができました。

それでも、補間による“滑らかすぎる値”が本当に正しいのか、中央値による“代表値”が意味のあるものなのか、常にどこかに「これはごまかしているだけでは？」という不安がつきまとっていました。それでも、試行錯誤しながら出力されたグラフに明らかな傾向や周期が見えてきたとき、「ああ、やっとここまできた」と実感することができました。

全体として、コードを書く行為そのものよりも、「正しいとは何か」を理解しようとする作業の方が何倍も難しく、そして根気が必要な工程でした。この経験を通じて、データ処理の技術的理解と同時に、正確な解釈や可視化の意義を考え抜くことの大切さを強く学んだように思います。

まとめ

今回は、真珠養殖の持続的な発展に向けて、新たな生産地の可能性を科学的に評価する方法として、衛星リモートセンシングによる海況の可視化に着目しました。

特に、第一歩としてGCOM-Cの海面水温（SST）データを用い、英虞湾周辺における時系列解析を行いました。

取り組みを通して、データ補間の必要性や、エラーDN値の処理の難しさなど、実際のデータ活用における課題にも直面しました。

しかし、Python上での線形補間処理や、周辺ピクセルの中央値を利用する工夫により、ある程度の有効データを抽出することができ、海面水温データの有用性を確かめる基盤を築くことができました。

今回の結果は、まだ限られた指標にとどまっていますが、海面水温という観点から真珠養殖環境の一端を把握する第一歩となりました。ただし、養殖の適否を判断するには、これだけでは不十分であり、今後は塩分濃度やクロロフィルa濃度、溶存酸素量といった他の重要な環境指標も併せて検討する必要があります。

特にクロロフィルa濃度は、アコヤ貝の餌となる植物プランクトンの量と関係しており、貝の成長や真珠の品質に直接影響する要素です。これら複数のデータを組み合わせることで、より現実的かつ信頼性の高い養殖適地の評価が可能となります。

また、衛星データを活用することで、実際に現地に足を運ばなくても、広範囲にわたる海域を対象に環境条件を把握・比較することができ、養殖候補地の選定において大きな可能性を秘めています。

そのため今後は、別の観測期間や季節、さらには別の地理的地域においても同様の解析を実施し、環境の時空間的な変動を踏まえた包括的な評価を行っていく必要があります。あわせて、他の物理量データとの統合的解析を通じて、最終的には真珠養殖のための環境モデル構築を目指していきたいと考えています。

環境変化に強い真珠養殖の実現に向けて、衛星データの活用は今後ますます重要になると考えられます。

コラム(8日平均)

　今回用いたLevel2のデータでは、雲や天候の影響により観測が困難となり、取得できるデータ数が少ないという課題がありました。
そこで、観測値を8日間ごとに平均処理した「8日平均SSTデータ」を活用し、日単位データに比べて雲や欠測の影響を受けにくく、より広範囲の水温分布を安定して得られるという特性を活かして、時系列グラフを作成してみました。以下に簡単にLevel2のデータと8日平均のデータの違いを表でまとめてみました。

https://suzaku.eorc.jaxa.jp/GCOM_C/data/product_def_j.html

8日平均データの場合、1か月あたりに得られるデータ数は3〜4点程度に限られます。そのため、ある程度の傾向を把握するために、3か月分のデータを用意しました。

もちろん、3か月分であっても依然としてデータ数は多くはなく、統計的な解析や明確な傾向の抽出には限界があります。

それでも、8日平均データは観測の安定性が高く、データ欠損が少ないという利点があるため、全体的な変動傾向や長期的な推移を初期的に把握するための素材としては十分に有効であると考えられます。

では、GCOM-C衛星が提供する8日平均の海面水温（SST）データを用いて、特定地点における時系列変化をグラフとして可視化するPythonコードを紹介します。

import h5py
import numpy as np
import pandas as pd
import glob
import os
import matplotlib.pyplot as plt
import re


# ── ユーザー設定 ─────────────────────
TARGET_LAT = 34.2
TARGET_LON = 136.7
DATA_FOLDER = "your_folder"  # 適宜パス変更
SLOPE = 0.0012
OFFSET = -10
# ────────────────────────────────────

まずは必要なライブラリのインポートと、観測対象地点やデータフォルダの指定を行います。ここでは、先ほどと同じように三重県沖の志摩半島周辺を想定し、緯度34.2°、 経度136.7°を指定しています。

def find_nearest_pixel(lat_array, lon_array, target_lat, target_lon):
    """指定座標に最も近いピクセルのインデックスを返す"""
    dist = (lat_array - target_lat)**2 + (lon_array - target_lon)**2
    i, j = np.unravel_index(np.argmin(dist), dist.shape)
    return i, j

次に、観測地点に最も近いピクセルを特定するための関数を定義します。緯度・経度の2次元配列から、指定座標に最も近いピクセルのインデックスを算出します。

def extract_sst_from_file(file_path, target_lat, target_lon):
    """指定ファイルから対象地点のSSTを取得し、物理量に変換"""
    with h5py.File(file_path, 'r') as f:
        sst_dn = f['Image_data']['SST_AVE'][:]


    rows, cols = sst_dn.shape
    lat = np.linspace(90, -90, rows)
    lon = np.linspace(0, 360, cols)
    lon2d, lat2d = np.meshgrid(lon, lat)


    i, j = find_nearest_pixel(lat2d, lon2d, target_lat, target_lon)
    dn = sst_dn[i, j]


    if dn in [65535, -9999]:
        return np.nan
    return dn * SLOPE + OFFSET

HDF5ファイルからSSTの8日平均値（SST_AVE）を抽出し、対象地点に最も近い値を取得・変換します。65535 や -9999 は欠損値としてNaNに変換しています。

ile_list = sorted(glob.glob(os.path.join(DATA_FOLDER, "*.h5")))
dates_asc, sst_asc = [], []
dates_desc, sst_desc = [], []


for file_path in file_list:
    basename = os.path.basename(file_path)


    # 日付抽出
    match = re.search(r"(\d{8})", basename)
    if not match:
        continue
    date_str = match.group(1)
    date = pd.to_datetime(date_str, format="%Y%m%d")


    # 軌道タイプ判定（8桁日付の直後の文字）
    idx = basename.find(date_str)
    if idx == -1 or len(basename) <= idx + 8:
        continue
    orbit_type = basename[idx + 8]  # 'A' または 'D'


    # SST抽出
    sst = extract_sst_from_file(file_path, TARGET_LAT, TARGET_LON)


    if orbit_type == "A":
        dates_asc.append(date)
        sst_asc.append(sst)
    elif orbit_type == "D":
        dates_desc.append(date)
        sst_desc.append(sst)

データフォルダ内のHDF5ファイルをすべて取得し、各ファイルの名前に含まれる8桁の日付と、その直後に記載された軌道情報（”A”：昇交軌道、”D”：降交軌道）をもとに、観測日時と軌道種別を判別し、水温データを昇交・降交軌道ごとに分類しています。

# ── DataFrame作成 ─────────────────────
df_asc = pd.DataFrame({"Date": dates_asc, "SST": sst_asc}).sort_values("Date").set_index("Date")
df_desc = pd.DataFrame({"Date": dates_desc, "SST": sst_desc}).sort_values("Date").set_index("Date")


# NaNを除去（dropna）して線が途切れないようにする
df_asc_clean = df_asc.dropna()
df_desc_clean = df_desc.dropna()

取得したデータを pandas の DataFrame にまとめ、日付をインデックスに設定します。また、グラフをプロットする際にNaNが存在するとグラフが途切れてしまうため、NaNを除外することでグラフ上で線が途切れないようにします。

plt.figure(figsize=(12, 5))
plt.plot(df_asc_clean.index, df_asc_clean["SST"], marker='o', linestyle='-', label="Ascending")
plt.plot(df_desc_clean.index, df_desc_clean["SST"], marker='s', linestyle='-', label="Descending", color='orange')
plt.title(f"SST at ({TARGET_LAT}, {TARGET_LON}) - Ascending vs Descending")
plt.xlabel("Date")
plt.ylabel("SST [°C]")
plt.grid(True)
plt.legend()
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

matplotlibを使用して昇交軌道と降交軌道のSSTを色分けしてプロットします。X軸には観測日、Y軸には海面水温（℃）を配置し、軌道ごとの違いが視覚的に確認できます。以上で8日平均時系列データ作成のコードは終わりです。以下にgooglecolabでこのコードを実行したときのグラフを添付しておきます。

本コードを実行し、得られたグラフを確認することで、2025年2月から4月にかけての約3か月間の8日平均SSTデータをもとに、昇交軌道（夜間の観測）と降交軌道（昼間の観測）における水温変化を比較することができました。

8日平均という手法の特性により、日ごとのばらつきや観測欠損の影響が抑えられ、なめらかな水温変化の傾向が視覚的に確認できる点も利点です。

一方で、1か月あたりのデータ点数は3〜4件程度にとどまり、時系列解析や統計的検定を行うにはややデータ量が不足している側面もあります。しかし、初期的な傾向把握や、観測地点における環境変化のモニタリングを目的とする場合には、十分有用な情報が得られると考えられます。

今後は、さらに長期間のデータを組み合わせたり、複数地点を比較することで、季節変動や地域差の分析にもつなげられる可能性があります。
Tellusの有料サービス「Tellus Satellite Data Master with QGIS」を使うとQGIS上で、海面水温やクロロフィルa濃度の時系列解析を行うことができるツールを利用することができます。

【関連記事】

最後に

今回の活動を踏まえて、気候変動や環境の変化が進む中で、衛星データは真珠養殖という伝統産業に新たな可能性をもたらすと思いました。

科学的な視点を取り入れることで、これまで見えていなかった適地や兆候を捉え、持続可能な未来を築いていく第一歩になると信じています。

今後は水温以外にも、塩分濃度、クロロフィルa濃度、溶存酸素量といった他の海洋環境要因も統合的に分析し、真珠養殖適地選定モデルの構築を行っていきたいです。

コードはGithubで公開しております。