いるかのボックス: Pythonで海面水温衛星データを可視化する

G-Portalという地球観測衛星データ提供システムがあり、登録するとさまざまな衛星データを取得できる。今回は気候変動観測衛星「しきさい」（GCOM-C）の海面水温プロダクト(SST) を取得して可視化してみる。

環境

WSL2(Ubuntu20.04)とJupyter。

$ lsb_release -dr
Description:    Ubuntu 20.04.4 LTS
Release:        20.04
$ python3 -V
Python 3.8.10
$ jupyter --version
jupyter core     : 4.7.1
jupyter-notebook : 6.2.0
qtconsole        : 5.0.2
ipython          : 7.20.0
ipykernel        : 5.4.3
jupyter client   : 6.1.12
jupyter lab      : 3.0.9
nbconvert        : 6.0.7
ipywidgets       : 7.6.3
nbformat         : 5.1.2
traitlets        : 5.0.5

G-Portalの登録

衛星データをダウンロードするために必要なユーザ登録をしておく。ユーザ登録画面から、アカウント名やメールアドレスなどを登録する。

衛星データのダウンロード

G-Portalトップページのメニューからログインにアクセスして、登録したユーザーアカウントとパスワードでログインする。［衛星からの検索］をクリックして、ダウンロードするデータを選択する画面に遷移する。

1.絞り込み

ダウンロードする衛星データの種類を指定する。

［1.絞り込み］タブの［衛星、センサから選ぶ］を選択する。そして、［GCOM-C/SGLI］＞［LEVEL2］＞［海洋圏］＞［L2-海面水温］のチェックボックスをチェックする。

2.期間指定

取得するデータを検索する期間を指定する。

［2.期間指定］タブで期間を指定する。

3.範囲指定

取得するデータの範囲を指定する。

［3.範囲指定］タブを開いて、右側に表示されている地図で指定するか、緯度・経度を数値で指定する。

検索するデータの指定ができたら［検索する］ボタンを押す。対象データの一覧が表示されるので、使用するデータをダウンロードする。ダウンロードしたファイル名はGC1SG1_L2SG_SSTDQ_3001.h5とする。

ダウンロードしたファイルはHDF形式で、可視化する際にはGeoTIFFが必要になる。データ一覧の［加工］ボタンを押すと、GeoTIFFへの加工要求ができる。加工には少し時間がかかるので、今回は後述するツールでGeoTIFFに変換する。

衛星データをGeoTIFFに変換する

ダウンロードしたHDF形式の衛星データをGeoTIFFに変換する。

まずはG-Portalから変換ツールをダウンロードする。トップページから［ツール・ドキュメント］画面へ遷移し、［GCOM-C］を選択。SGLI 地図投影・GeoTIFF出力ツール(Ver.1.1)をダウンロードする。ここではLinux版を使う。

ファイル名をsgli_geo_map_linuxとしてダウンロードし、先にダウンロードした衛星データを、WSL2上の同じディレクトリに配置する。

変換する前に、ダウンロードした衛星データ内のデータセット名を確認しておく。

! pip3 install h5py

import h5py

def print_dataset_name(name, obj):
    if isinstance(obj, h5py.Dataset):
        print(name)

with h5py.File('./GC1SG1_L2SG_SSTDQ_3001.h5', mode='r') as f:
    f.visititems(print_dataset_name)

上記コードを実行すると、以下のようにHDFファイル内のデータセット一覧が表示される。

Geometry_data/Latitude
Geometry_data/Longitude
Geometry_data/Obs_time
Geometry_data/Sensor_zenith
Geometry_data/Solar_zenith
Image_data/Cloud_probability
Image_data/Line_tai93
Image_data/QA_flag
Image_data/SST

上記データセットのうち、今回使用する海面水温プロダクト(SST)はImage_data/SST。変換時にはこのデータセットを指定する。

以下コマンドでHDF形式からGeoTIFFに変換する。はじめにダウンロードしたツールに実行権限を付与しておく。変換されたGeoTIFFのファイル名はGC1SG1_L2SG_SSTDQ_3001_SST.tif となる。

$ chmod u+x sgli_geo_map_linux
$ ./sgli_geo_map_linux GC1SG1_L2SG_SSTDQ_3001.h5 -d Image_data/SST

SSTデータの値を摂氏（℃）に変換

SSTの値は一般的に使われる温度の単位（摂氏）ではないので、摂氏に変換する必要がある。しきさい（GCOM-C）のページにSSTデータの仕様とその変換方法が記載されている。トップページから［データを使う］＞［標準プロダクト&アルゴリズム］を開き、Ver.3の海面水温を選択する。開いた画面のAttribute情報にあるSSTを開く。

この仕様によると、有効な値は0~65531で、65532～65535は陸地や雲、その他のエラーなどで水温が取得できなかった場合。SSTの値を摂氏（℃）に変換する式は以下の通り。DNがSSTの値で、Slopeは0.0012、Offsetは-10となっている。

SST[degree]=DN*Slope+Offset

海面水温データを可視化する

ダウンロードした海面水温衛星データを、Jupyterで可視化してみる。地理空間データ用のPythonライブラリであるrasterioとEarthPyを使う。

! pip3 install rasterio # ラスター形式の地理空間データを扱うためのPythonライブラリ
! pip3 install earthpy # ラスターおよびベクター形式空間データを扱うためのPythonライブラリ

import rasterio
from earthpy.mask import mask_pixels
from earthpy.plot  import plot_bands

file_path = './GC1SG1_202206100147B05610_L2SG_SSTDQ_3001_SST.tif'
sst = rasterio.open(file_path)

# バンド番号（1からはじまる）を指定してデータ読み込み
# https://rasterio.readthedocs.io/en/latest/topics/reading.html
sst_array = sst.read(1)
print(sst_array.shape)

# SSTの非有効値をマスクする
# https://earthpy.readthedocs.io/en/latest/api/earthpy.mask.html
mask_values = [65532, 65533, 65534, 65535] # マスクする値
sst_masked = mask_pixels(sst_array, sst_array, vals=mask_values)

# SSTの値を摂氏（℃）に変換
slope = 0.0012
offset = -10.0
sst_deg = sst_masked * slope + offset
sst_deg_min = sst_deg.min() # 最低水温
sst_deg_max = sst_deg.max() # 最高水温

# 海面水温を可視化
plot_bands(sst_deg, vmin=sst_deg_min, vmax=sst_deg_max, cmap='rainbow')

以下のような画像が表示される。陸地や雲などで有効値がない箇所は白色になっている。

いるかのボックス

2022年6月14日火曜日

Pythonで海面水温衛星データを可視化する

環境