衛星データを使って琵琶湖の温度を観察する #Python

衛星データを扱うお仕事をしていますが、最近「海面水温とかとれるの？」と聞かれました。ということで、試しにやってみました。

感想としては、rioxarrayがあれば、簡単にこんな処理もできるなぁということでした。ありがとうrioxarray！！！

使ったデータ・開発環境など

使ったデータ

GCOM-C(しきさい)
海面水温
JAXA G-PORTALより取得

開発環境

google colab (Python 3.7.14)
rioxarray 0.9.1
geopandas 0.10.2
plotly 5.5.0

JAXAのGPORTALから画像をダウンロード

ユーザ登録

JAXAのGPORTALはIDを作ると衛星のデータが無償でダウンロードできます。ユーザー登録はここから出来るようです。

ダウンロードしたいデータを選び、ダウンロードする

登録が終わったら、GUIからデータを選択しダウンロードします（ここのページから）。

範囲は大きめに指定します。そうすると、候補のデータをあげてくれます。

G-PORTALより

よさげなのを加工してダウンロードします。今回はgeotiff形式で扱いたいので、その指定を行います。

G-PORTALより

ダウンロードしたデータを読み込み、指定範囲のデータを作る

データがダウンロードできました。ダウンロードしたgeotiffファイルをrioxarrayのopen_rasterioで読み込みます。

import rioxarray as rxr
import geopandas as gpd

data = rxr.open_rasterio('/data/GC1SG1_202209200148J05610_L2SG_SSTDK_3001_0000000001010830_SST.tif')

読み込み終わったデータを確認し、可視化してみます。

データ確認

読み込んだデータはx軸に経度、y軸に緯度を持つデータということが分かります。

データ可視化

data.plot()

気温のデータがちょっと変

気温のデータに6万を超える値が入っていました。そこでググってみると、QGISを使ったGCOM-C（しきさい）プロダクトの画像化手順というページが出てきました。

そちらの資料に画像データを海面水温に換算するという項目がありました。

次の数式で換算するようです。

SST(℃) = DN × 0.0012 - 10

というわけで、先ほど読み込んだデータも水温に換算します。

data = data * 0.0012 -10

換算することで、次のように気温っぽいデータが作れました。

琵琶湖のデータをとってきて範囲指定したデータを作成する

気温のデータはできたので、最後に琵琶湖のデータを取得し、範囲を切り抜き、その湖面気温を作成します。

琵琶湖の位置情報に関しては、はんなりPythonのGISハンズオンのドリさんパートのノートブックを参照ください。国土交通省のページから取得できます。

そのデータを使ってと思いましたが、ちょっとカクカクし過ぎなので時間がかかりそうなので、そのデータから琵琶湖を囲める四角形を作って、その場所のデータを代替で取得しました。

lake = gpd.read_file('/content/W09-05-g_Lake.shp')
biwako = lake.loc[418]
biwako_rote = biwako.loc[418, 'geometry'].minimum_rotated_rectangle # ここで琵琶湖を囲む四角形を作った
print(biwako_rote.__geo_interface__)

出力

ざっくりと琵琶湖風の形を切り出すために、データを処理します。

biwako_rote_bounds = biwako_rote.bounds

これを使って、tifデータを切り抜くと、琵琶湖周辺のデータが作れます。

コードをまとめていく

琵琶湖の水温用のコードをまとめます。


data = rxr.open_rasterio('GC1SG1_202209200148J05610_L2SG_SSTDK_3001_0000000001010830_SST.tif')
data = data.where(data < 60000) # 60000以上の数値をnaとする
data = data * 0.0012 - 10
lake = gpd.read_file('/content/W09-05-g_Lake.shp')
biwako = lake.loc[418]
biwako_rote = biwako['geometry'].minimum_rotated_rectangle
biwako_rote_bounds = biwako_rote.bounds
biwako_data = data.sel(
    x=slice(biwako_rote_bounds[0], biwako_rote_bounds[2]),
    y=slice(biwako_rote_bounds[3], biwako_rote_bounds[1])
)
biwako_data.plot()

可視化すると琵琶湖の形っぽくないが、一応データはできたようです。

QGISで位置があっているか確認する

うまく処理できたか分からないので、一先ずtifファイルで保存し、QGISでいい感じにデータが表示されるか確認します。

biwako_data.to_raster('biwako.tif') #保存

QGISでオープンストリートマップにかぶせてみます。

このデータが取得できていないのは何なのか？琵琶湖クリップをしていないデータを見てみると・・・

6万以上の数値には意味があった！！！（ドキュメントしっかり読め問題）

やはり、ちょっと抜けがあります。これは雲がかかっている可能性が考えられますね。と思ってJAXAの資料を見返していると、格納されている60000以上の数値は次のような意味があるようです。

最後にしきさいのデータを扱う際に注意すべき点として資料よりピックアップしておきます。

QGISを使ったGCOM-C(しきさい)プロダクトの画像化手順（STEP1-2）

アピール

衛星データを使ってみたいとか、会社のDXを進めたいとかございましたら、京都のデータテクノロジー企業　長目までおといあわせください。