漁業での衛星データ利用事例

この章で学習すること

  • 「しきさい」「ひまわり」衛星で観測した海面水温画像の可視化
  • rasterioEarthPyを使ったGeoTiffファイルの操作
  • matplotlibrasterioCartopy を組み合わせた衛星データの可視化
  • GeoPandas を使ったShapeファイル、GeoJSONファイルの操作

背景

持続可能な開発目標(SDGs)では、目標14に「海の豊かさを守ろう」が定められています。海洋と海洋資源を保全し、持続可能な形で利用しようというもので、例えば以下のようなターゲットが含まれています(数字はターゲット番号)。

  • 海洋ごみや富栄養化など海洋汚染の防止・削減(14.1)
  • 海洋及び沿岸の生態系の管理と保護(14.2)
  • 乱獲や違法・無報告・無規制(IUU)漁業および破壊的な漁業慣行の撤廃(14.4)
  • 科学的情報に基づく沿岸・海洋エリアの保全(14.6)
  • 漁業や養殖などを通じた海洋資源の持続的な利用による経済的利益の増加(14.7)

地球の7割以上を占める海洋において、効率的に観測して状況や変化を把握するためには地球観測衛星が力を発揮します。既に、富栄養化や重油流出事故などの海洋汚染は衛星でも観測されていますし、サンゴ礁など沿岸の生態系把握にも衛星データが利用されています。

この章では、衛星データを用いて漁場を推定する方法をご紹介します。漁場を推定することで、乱獲やIUU漁業を効果的に取り締まることが可能となりますし(14.4)、時間と燃油を節約できるので経済的利益の増加に繋がります(14.7)。漁業者が積極的に資源管理に取り組むことになり、海洋資源の持続的な利用に繋がることも期待されます。ここでは、衛星リモートセンシングによる海面水温のデータを用いて、Pythonの様々なライブラリを使った地理空間データのハンドリング・可視化について学んでいきます。

海面水温(SST: Sea Surface Temperature)と漁場の関係

漁業の現場では、以前から衛星リモートセンシングで観測した海面水温のデータが利用されてきました。これは好漁場と水温が密接に関係していると言われているからです。

変温動物である魚には種類によって適正水温があり、水温の分布パターンはターゲットの魚種が集まる海域の目安の一つとなります。また、異なる密度の水の塊がぶつかりあう場所である「潮目」は、魚の餌となる植物性プランクトンが集まりやすい場所となっており、経験的にも好漁場であることがよく知られています。

下の図は衛星で観測した海面水温分布図に、漁場を◯と△で重ねた図です。温度が急激に変化する潮目の部分に漁場が集中していることがわかります。

JAFIC_潮目.png

図:2012 年 10 月 26 日の東北沖衛星水温と漁場

水産海洋分野の衛星リモート センシングとICT - 東京水産振興会 図4より引用)

リモートセンシングによる海面水温の観測

赤外線の観測

海面水温は熱赤外線やマイクロ波で観測したデータを元に算出をしています。これらの観測は、太陽光の反射を計測する可視・近赤外線観測と異なり、物体自身が放射する電磁波を計測しており、これを利用して物体の温度を推定することができます。

大気の影響の補正

海面から放射される赤外線は大気の影響を受けて減衰し衛星センサに到達するため、観測する温度は実際の海面水温よりも低くなります。この影響を補正するために、「しきさい」や「ひまわり」では、波長の異なる2つ以上の赤外線チャンネルで得られた観測値の差から大気による減衰分を推定し、海面の水温を推定しています。

海面水温が観測できるセンサ

海面水温が観測できるセンサは、熱赤外線センサを搭載した「しきさい」や、「ひまわり」の他、マイクロ波センサを搭載した「しずく」等があります。マイクロ波センサは熱赤外線センサと比べて空間分解能は劣りますが、雲の下を観測できるという利点があります。

海面水温が観測できる衛星・センサは、宙畑のこちらの記事でも紹介されているので、ぜひご覧ください。

コラム1:リモートセンシングによる海色観測

人工衛星による海洋観測は、海面水温だけではありません。例えば、今回取り上げた「しきさい」では、植物プランクトン現存量の指標になるクロロフィル-a濃度の分布を表したプロダクトも公開されています。プランクトンの異常な増殖により水の色が変化する現象を「赤潮」と呼ぶように、プランクトンの量は海色に大きく影響を与えています。衛星ではこの海色を、可視光をベースに観測しています。

海の中では、植物プランクトン→動物プランクトン→イワシなどの小魚→大きい魚という食物連鎖があり、植物プランクトンが多い場所では漁の対象となる魚が集まっていることが予想されます。好漁場を推定する研究では、海面水温などの物理環境が利用されるケースが先行していますが、上記のクロロフィル-a濃度のような生物環境の条件を考慮したモデルの研究も進んでいます。

スクリーンショット 2020-12-02 165824.png

しきさいの海面水温とクロロフィル-aプロダクトは、Tellus OSでも閲覧できます。

しきさいで観測したSSTの解析

しきさい(GCOM-C)で観測した海面水温プロダクトを題材に、 Pythonを使った衛星データの可視化について学習していきます。

しきさいについて

しきさい/GCOM-C(Global Change Observation Mission-Climate)は、2017年にJAXAが打ち上げた地球観測衛星です。可視光から熱赤外まで19の観測チャンネルを持ち、地球の気候形成に影響を及ぼしている様々な物理量の観測を行っています。空間分解能は最高で250mです。

では、熱赤外(TIR)チャンネルを用いて観測したデータから作成された、海面水温のデータを使用していきます。

準備

In [ ]:
# google driveのマウント
# Tellusの開発環境やローカル環境を使用する想定であれば、このセルは削除可です。
from google.colab import drive
drive.mount('/content/drive')
Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).

必要なモジュールのインストール

In [ ]:
!apt-get install libproj-dev proj-data proj-bin
!apt-get install libgeos-dev
!pip install cython
!pip install cartopy
!apt-get -qq install python-cartopy python3-cartopy
!pip uninstall -y shapely    # cartopy and shapely aren't friends (early 2020)
!pip install shapely --no-binary shapely
In [ ]:
!pip install geopandas rasterio earthpy

データのダウンロード

今回はG-Portalからデータをダウンロードします。

大まかな手順は以下の通りです。

  1. ログインし、トップページから「衛星からの検索」をクリックする
  2. 「1. 絞り込み」の「衛星、センサから選ぶ」タブから、GCOM-C/SGLI -> LEVEL2 -> 海洋圏 -> L2-海面水温を選択する
  3. 「2. 期間指定」から希望の期間を指定する
  4. 「3. 範囲指定」の「矩形指定」で、地図上で希望の範囲をドラック指定する
  5. フォーマットを「GeoTiff」に指定し、加工要求を行う
  6. 生産終了後、通知メールに記載のリンクからデータをダウンロードする (ログイン後のトップページからもダウンロード可能です)

※ ダウンロードの方法について、詳細はこちらをご覧ください。

使用するデータの情報

  • 処理レベル: L2
  • プロダクトの種類:海面水温
  • 観測日:2018年11月04日
  • 空間分解能:250m

※ 空間分解能が小さいと、環境によってはメモリが足りずに処理が落ちてしまうことがあります。その場合は1kmのものをご使用ください。

使用する主なライブラリ

「基礎編」でご紹介したGDALライブラリは、地理空間データを扱う上で基礎となる重要なライブラリですが、GISの概念に慣れていないと理解しにくいパラメータがたくさん登場するため、衛星データを触り始めたばかりの頃は複雑に感じることもあるかと思います。

この教材では、地理空間データの処理がより便利になるPythonのライブラリをいくつかご紹介します。

  • rasterio

    rasterioは、オンラインマッププロバイダーのMapboxが中心となって開発しているOSSのライブラリで、衛星データでは最も一般的なGeoTiffフォーマットのデータを扱う際に、シンプルな書き方ができるのが特徴のライブラリです。

  • EarthPy

    米コロラド大学ボルダー校のEarthlabが中心となって開発しているライブラリで、特定のデータ操作についてはrasterioと比べてより直感的に行なえます。

  • Cartopy

    Pythonでのデータ可視化に欠かせないmatplotlibと連携して、地理空間データの描画ができるライブラリです。様々な地図投影法が用意されている他、国境線・海岸線・緯経線など独自のベクターデータも利用できます。

  • GeoPandas

    有名なデータ処理ライブラリpandasを、地理情報が取り扱いできるように拡張したライブラリです。pandasと同じ強力なデータ処理機能を持つだけでなく、他の地理空間データとの連携が可能です。

In [ ]:
 # 必要ライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt

# rasterio
import rasterio
from rasterio.plot import show

# earthpy
import earthpy.mask as em

# cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# GeoPandas
import geopandas as gpd

データの読み込みと可視化

それでは早速SSTデータを可視化していきましょう。

今回の対象はGeoTiffです。rasterioを使って簡単に可視化ができます。

In [ ]:
file_path = "/content/drive/MyDrive/GCOM-C/GC1SG1_201811040128M05010_L2SG_SSTDQ_1002_0000000000192925_SST.tif"

# rasterio.open(読み込むファイルパスを指定)
sst_object = rasterio.open(file_path)
fig, ax = plt.subplots(1,1, figsize=(10,10))
show(sst_object, ax=ax)
Out[ ]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fac2f9fee10>

このように、数行のコードで簡単に画像の表示までできました!

上の図では、縦軸が緯度、横軸が経度になっています。
日本の陸の形はわかりますが、これだけではどのようなデータになっているか分からないので、もう少し見やすくしていきましょう。

まずはSSTデータのフォーマットを確認します。

GCOM-Cプロジェクトのホームページには、SSTプロダクトのデータ仕様が記載されています。

Attribute情報->SSTを確認すると、画素値のうちいくつかには以下のように意味が与えられていることがわかります。

(一部抜粋)

key value
Error_DN 65535
Land_DN 65534
Cloud_error_DN 65533
Retrieval_error_DN 65532
Maximum_valid_DN 65531
Minimum_valid_DN 0

また、データとして有効な値の範囲は0~65531であることもわかります。

陸や雲、エラーの領域は解析で不必要な領域になるので、マスクしてしまいましょう。 マスクとは、画像の指定した領域を隠し、可視化や解析の対象外とする処理です。

画像のマスクには、earthpyを使います。

In [ ]:
# 元データの座標変換情報を保持
sst_trans = sst_object.transform

# マスクする値を設定
mask_values = [65535, 65534, 65533, 65532]

# 画像をnumpy.ndarrayとして取得
# rasterio_object.read(バンド番号)で読み込むバンド番号を指定
sst_dn = sst_object.read(1) # dn: digital number

# マスク処理
# em.mask_pixels(マスク対象の配列、マスクに使用する値を持つ配列、vals=[マスクする値])
sst_masked = em.mask_pixels(sst_dn, sst_dn, vals=mask_values)

画素値から温度[℃]への変換は、同じくAttribute情報->SSTから、以下の式で温度(SST)に変換できることがわかります。

$ SST[degree]=DN*Slope+Offset $

Slope, Offsetの値は、同じくAttribute情報->SSTから、以下の値になっています。
(一部抜粋)

key value
Slope 0.0012
Offset -10
In [ ]:
# DN値をSSTに変換
slope = 0.0012
offset = -10.0

sst_temp = sst_masked * slope + offset

色付けの仕方も、デフォルトではなく温度が感覚的にわかりやすいものに置き換えます。今回は、matplotlib既定のjetというカラーマップを使います。

In [ ]:
fig, ax = plt.subplots(1,1, figsize=(10,10))
ax.set_xlabel('longtitude[degree]') # X軸名を追加
ax.set_ylabel('latitude[degree]')   # Y軸名を追加
# show()でnp.ndarrayを指定する場合は、"transform="でrasterio_object.transformを指定する
ax = show(sst_temp, transform=sst_trans, cmap='jet')

画像の範囲外と陸・雲の領域が白くなっており、マスク処理に成功したことがわかります。

画素値がどのように分布しているか確認するために、ヒストグラムを作成します。
ここではrasterio.plotshow_histという関数を使用します。

In [ ]:
from rasterio.plot import show_hist

fig, ax = plt.subplots(1,1, figsize=(10,5))

show_hist(
    sst_temp, bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", label="SST", ax=ax)
ax.set_xlabel("Temperature[℃]")
Out[ ]:
Text(0.5, 0, 'Temperature[℃]')

この画像の範囲では、海面水温は0~30℃に分布していることがわかります。

それでは、こちらの温度範囲に合わせてカラー表示をしてみましょう。

In [ ]:
temp_min = 5.0 # deg.C
temp_max = 28.0 # deg.C

# 表示範囲の最大値、最小値の設定
options = {
    'vmin': temp_min,
    'vmax': temp_max,
    'cmap': 'jet',
}

fig, ax = plt.subplots(1,1, figsize=(10,10))
ax.set_xlabel('longtitude[degree]')
ax.set_ylabel('latitude[degree]')
show(sst_temp, transform=sst_trans, ax=ax, **options)
Out[ ]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff11648b3c8>

温度変化のある部分がよりはっきりと分かります。

ちなみに、earthpyでは次のように更に簡潔に可視化できます.

In [ ]:
import earthpy.plot as ep
ep.plot_bands(sst_temp, **options)
Out[ ]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff1176b17f0>

シンプルに記述でき、デフォルトでカラーバーも表示してくれます。

コラム2:データの品質について


今回使用している「しきさい」のデータは、G-Portalで加工要求を行ったGeoTiffファイルです。オリジナルのデータ(HDF5フォーマット)には、海面水温の値以外にも「QAフラグ」というデータが含まれています。

このQAフラグには、上記でマスクした陸・雲の領域だけでなく、アルゴリズムで雲の有無の判断が困難だった領域や、雲の可能性がある領域などが示されており、既存より厳しく雲のマスクをかけたい時に有効活用できます。

詳細については、JAXA/GCOM-Cプロジェクトのよくある質問 | QAフラグとはをご覧ください

潮目の確認

拡大表示と表示する温度によって、コントラストの違いから潮目の見え方が違うことを確認します。

今回は、黒潮と親潮がぶつかり好漁場となる、岩手宮城・三陸沿岸にフォーカスしてみましょう。

TellusOSを利用して、関心領域のGeoJSONファイルを作成します。

まず、TellusOSにアクセスし、地図上クリップする領域を描画します。

Inkedtellus_aoi_LI.jpg

画面右上に出現するウィンドウのダウンロードボタンをクリックし、 GeoJSONフォーマットを選択してダウンロードします。

選択した領域を修正したい場合は、「編集」ボタンから調整ができます。

画像保存.png

ダウンロードしたGeoJSONファイルを使って、SSTの表示範囲を決定します。

In [ ]:
import cartopy.crs as ccrs
import geopandas as gpd
import cartopy.feature as cfeature

# geojsonをインポート
fname_aoi =  "/content/drive/MyDrive/GCOM-C/rectangle-1.geojson"
df_aoi = gpd.read_file(fname_aoi, encoding='utf-8')

# geometoryから範囲を取得
bounds = df_aoi['geometry'].bounds

fig = plt.figure(figsize=(21,7))

# パターン1
temp_min = 5.0 # deg.C
temp_max = 28.0 # deg.C
ax1 = plt.subplot(1,3,1, projection=ccrs.PlateCarree()) # 座標系をPlateCarreeに指定
ax1.set_extent((bounds['minx'], bounds['maxx'], bounds['miny'], bounds['maxy']))
ax1.set_title('SST: ' + str(temp_min) + '-' + str(temp_max) + '℃')
show(sst_temp, transform=sst_trans, ax=ax1, vmin=temp_min, vmax=temp_max, cmap='jet')
ax1.add_feature(cfeature.COASTLINE) # 海岸線を追加

# パターン2
temp_min = 15.0 # deg.C
temp_max = 20.0 # deg.C
ax2 = plt.subplot(1,3,2, projection=ccrs.PlateCarree())
ax2.set_extent([bounds['minx'], bounds['maxx'], bounds['miny'], bounds['maxy']])
ax2.set_title('SST: ' + str(temp_min) + '-' + str(temp_max) + '℃')
show(sst_temp, transform=sst_trans, ax=ax2, vmin=temp_min, vmax=temp_max, cmap='jet')
ax2.add_feature(cfeature.COASTLINE)

# パターン3
temp_min = 18.0 # deg.C
temp_max = 23.0 # deg.C
ax3 = plt.subplot(1,3,3, projection=ccrs.PlateCarree())
ax3.set_extent([bounds['minx'], bounds['maxx'], bounds['miny'], bounds['maxy']])
ax3.set_title('SST: ' + str(temp_min) + '-' + str(temp_max) + '℃')
show(sst_temp, transform=sst_trans, ax=ax3, vmin=temp_min, vmax=temp_max, cmap='jet')
ax3.add_feature(cfeature.COASTLINE)
Out[ ]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7ff12158fb38>

異なるデータを重ねて表示させるには、投影する座標を揃えなければいけません。ここではcartopyを使って、matplotlib上の座標を正距円筒図法(Plate Carree)という、緯線と経線が直角・等間隔に交差する座標系を指定しています。同時にcartopyで用意している海岸線も合わせて表示しています。

左の図は全体表示させたときと同じ温度範囲ですが、拡大すると同じ色合いになっており、局所的な変化を確認するには適当でないことがわかります。

中央の図は、ヒストグラムをみて温度変化が大きそうな範囲を可視化しました。暖かい部分と冷たい部分の領域がはっきりとわかります。

右の図は、高い温度の範囲の変化に焦点を当ててみました。画面の中央では赤く塗りつぶされている領域内でも、温度の変化があることがわかります。

海流のプロット

潮目は、水温や塩分の違いなどにより密度が異なる水の塊がぶつかるところに発生し、海流と深く関係しています。

海面水温の画像と海流のデータを重ねて見てみましょう。

海流のデータは、海上保安庁の海洋速報&海流推測図のページから、シェープファイルが取得できます。今回は津軽暖流のデータを使います。

2019年11月4日の海況図は以下になります。

Shapefileの読み込みにはGeoPandasを使います。

In [ ]:
# 津軽暖流のシェープファイルをインポート
file_tsugaru =  "/content/drive/MyDrive/GCOM-C/qboc2018Nov/tg/qboc2018208tg_Line.shp"
df_tsugaru = gpd.read_file(file_tsugaru)

# SSTを表示
plt.figure(figsize=(10, 10))

# 地図投影法をCartopyで指定
ax = plt.axes(projection=ccrs.PlateCarree())
temp_min = 10.0 # deg.C
temp_max = 20.0 # deg.C
show(sst_temp, transform=sst_trans, ax=ax, vmin=temp_min, vmax=temp_max, cmap='jet')
ax.add_feature(cfeature.COASTLINE)
df_tsugaru.plot(ax=ax, color='purple', linewidth=3.0, label="Tsugaru Warm Current")
ax.set_extent([500.66345214843744, 503.99230957031244, 40.55137419871514, 42.65820178455669])
ax.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=18)
ax.set_title('SST with Tsugaru Warm Current')
Out[ ]:
Text(0.5, 1.0, 'SST with Tsugaru Warm Current')

上の図で紫色の線が津軽暖流です。海流図と比べて見ると津軽海峡→太平洋の方向の流れであることがわかります。このように、海面水温と海流のシェープファイルを重ねることで、海流の流れにそって暖かい水が流れている様子が確認できました。

ひまわりSSTの解析


ひまわり8/9号について

ひまわり8/9号は気象庁がもつ気象衛星で、静止軌道上からで日本を含むアジア・太平洋地域を広く観測しています。10分に一度と高頻度に観測できるのが特徴で、時々刻々と変化する雲の動き等を観察するのに適しています。その分、地球から遠く離れた軌道上にいるため、「しきさい」に比べて空間分解能が劣ります。観測波長帯は可視~近赤外~熱赤外と「しきさい」と同じような帯域を持っています。

ひまわりSSTについて

Tellusマーケットで公開予定のひまわりSSTデータ(©RESTEC, 気象庁データからRESTECが作成)は、10分に一度観測されるデータを3時間分合成し、雲で観測できない領域を埋めています。

使用するデータの情報

  • プロダクトの種類:海面水温(3時間合成)
  • 観測日:2019年11月12日 00:00-24:00(JST)間の8枚の画像を使用
  • 空間分解能:約2km

ファイルの命名規則は、

  • SST_3H_{YYYYMMDD}_{HHMM}_{hhmm}_{X}.tif
  • YYYYMMDD: データ取得年月日
  • HHMM: データ合成開始時刻
  • hhmm: データ合成終了時刻
  • X: プロダクトバージョン

になります。10分毎に観測したデータを3H分合成して雲の欠損を埋めたデータです。

本教材用にサンプルデータをこちらのリンクからダウンロードいただけますので、ご利用ください。当該データの利用で発生した事故等で利用者が受けた損害について、教材提供者は一切の責任を負いませんので予めご了承ください。なお、利用に当たって下記の点に留意してください。

  • 学習のための内部利用に限定し、二次配布は禁止

読み込みと可視化

ファイルフォーマットはGeoTiffですので、基本はしきさいと一緒です。 rasterioを使ってファイルの読み込み・可視化を行います。

In [ ]:
# モジュールのインポート
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
import earthpy.mask as em
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import geopandas as gpd
from osgeo import gdal
In [ ]:
file_path = "/content/drive/MyDrive/Himawari/SST_3H_20191113_0000_0250_A.tif"

with rasterio.open(file_path) as himawari_sst:
    himawari_trans = himawari_sst.transform # 座標変換情報
    himawari_sst_array = himawari_sst.read(1) # 画像情報

# マスク処理
mask_values = [-9999, -8888] # (雲、陸)
himawari_sst_masked = em.mask_pixels(himawari_sst_array, himawari_sst_array, vals=mask_values)
In [ ]:
# 画像を表示
fig, ax = plt.subplots(1,1, figsize=(10,10))
ax.set_xlabel('longtitude[degree]')
ax.set_ylabel('latitude[degree]')
show(himawari_sst_masked, transform=himawari_trans, ax=ax, cmap='jet')
Out[ ]:
Text(0, 0.5, 'latitude[degree]')

画像の下に見えるのがオースラリア大陸です。緯度と経度の範囲からも、しきさいと比べて広い範囲を撮影していることがわかると思います。 そのため、画像のピクセル数が多く処理に時間がかかるので、予め解析に使う領域だけ取り出して処理したいと思います。

切り出しの方法は、基礎編と同じくGDALライブラリを使います。

In [ ]:
def clip_tif(src, bounds):
    """GeoTiffファイルを範囲指定して読み込む

    Args:
        src(str): ソースファイルのパス
        bounds(list): 切り取り範囲 [ulx, uly, lrx, lry]

    Returns:
        clipped_array(numpy.array): 切り取り後の画像
        clipped_array(rasterio.Affine): 切り取り後の座標変換情報

    """
    trans_options = gdal.TranslateOptions(
        format = 'vrt',
        projWin = bounds,
        resampleAlg = 'nearest'
    )
    ds = gdal.Translate('', src, options=trans_options)

    # 画像データをnumpy.array型で保持
    clipped_array = ds.ReadAsArray()
    # 座標変換情報をrasterio.Affine型で保持 
    clipped_trans = rasterio.Affine.from_gdal(*ds.GetGeoTransform())

    return clipped_array, clipped_trans
In [ ]:
# 日本近辺
minx = 128.0
maxx = 150.0
miny = 29.0
maxy = 46.0
bounds = [minx, maxy, maxx, miny] # 切り取り範囲を指定

file_path = "/content/drive/MyDrive/Himawari/SST_3H_20191113_0000_0250_A.tif"
# 上のセルで作成したclip_tif関数を使用
himawari_sst_clipped, himawari_trans_clipped = clip_tif(file_path, bounds)
In [ ]:
# マスク処理
mask_values = [-9999, -8888] # (雲、陸)
himawari_sst_clipped_masked = em.mask_pixels(himawari_sst_clipped, himawari_sst_clipped, vals=mask_values)

# 可視化
fig, ax = plt.subplots(1,1, figsize=(10,10))
show(himawari_sst_clipped_masked, transform=himawari_trans_clipped, ax=ax, cmap='jet')
ax.set_xlabel('longtitude[degree]')
ax.set_ylabel('latitude[degree]')
Out[ ]:
Text(0, 0.5, 'latitude[degree]')

日本近くの領域に拡大できました。

時系列変化の可視化

ひまわりは高頻度に観測をしているため、時間変化を観測することができることができます。

ファイル命名規則は以下のとおりです。

SST_3H_YYYYMMDD_hhmm_hhmm_A.tif

  • YYYYMMDD:観測日(YYYY:年、MM:月、DD:日)
  • hhmm_hhmm:観測開始~終了時刻(hh:時間、mm:日時)
In [ ]:
import glob

# 特定の日付のファイルリストを取得
file_list = glob.glob('/content/drive/MyDrive/Himawari/SST_3H_201911*.tif')
In [ ]:
# ファイル一覧の表示
for file in file_list:
    print(file)
/content/drive/MyDrive/Himawari/SST_3H_20191112_1500_1750_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191112_1800_2050_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191112_2100_2350_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191113_0000_0250_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191113_0300_0550_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191113_0600_0850_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191113_0900_1150_A.tif
/content/drive/MyDrive/Himawari/SST_3H_20191113_1200_1450_A.tif

それでは、3時間合成のデータを一日分(8枚分)連続的に表示させて、海面水温の変化を可視化していきます。

In [ ]:
mask_values = [-9999, -8888] # (雲、陸)
temp_min = 5.0 # deg.C
temp_max = 28.0 # deg.C
options = {
    'vmin': temp_min,
    'vmax': temp_max,
    'cmap': 'jet',
}

fig = plt.figure(figsize=(7, 5))
images = []
for file_path in file_list: 
    # 画像ファイルを読み込み
    time = file_path[-24:-6]
    himawari_sst_clipped, himawari_trans_clipped = clip_tif(file_path, bounds)
    himawari_sst_clipped_masked = em.mask_pixels(himawari_sst_clipped, himawari_sst_clipped, vals=mask_values)
    img = plt.imshow(himawari_sst_clipped_masked, animated=True, **options)
    plt.axis('off')
    images.append([img])
plt.close()
In [ ]:
from matplotlib import animation, rc

# アニメーションの作成
am = animation.ArtistAnimation(fig, images, interval=1500, repeat_delay=2000)
rc('animation', html='jshtml')
am
Out[ ]:

時刻が進むにつれて温度分布が変化しており、また雲が動いている様子もわかります。

さらに拡大して時刻変化をみてみましょう。

In [ ]:
# 岩手・宮城三陸沖のGeoJSONをインポート
# しきさいで使用したものと同じ
fname_aoi =  "/content/drive/MyDrive/GCOM-C/rectangle-1.geojson"
df_aoi = gpd.read_file(fname_aoi, encoding='utf-8')

# geometoryから範囲を取得
bounds = sum(df_aoi['geometry'].bounds.values.tolist(), [])
bounds[1], bounds[3] = bounds[3], bounds[1] # gdal向けに順番入れ替え
In [ ]:
mask_values = [-9999, -8888] # (雲、陸)
temp_min = 10.0 # deg.C
temp_max = 22.0 # deg.C
options = {
    'vmin': temp_min,
    'vmax': temp_max,
    'cmap': 'jet',
}

images = []
fig = plt.figure(figsize=(7, 5))
for file_path in file_list: 
    # 画像ファイルを読み込み
    time = file_path[-24:-6]
    himawari_sst_clipped, himawari_trans_clipped = clip_tif(file_path, bounds)
    himawari_sst_clipped_masked = em.mask_pixels(himawari_sst_clipped, himawari_sst_clipped, vals=mask_values)
    img = plt.imshow(himawari_sst_clipped_masked, **options, animated=True)
    plt.axis('off')
    images.append([img])
plt.close()
In [ ]:
# アニメーションの作成
am = animation.ArtistAnimation(fig, images, interval=1500, repeat_delay=2000)
rc('animation', html='jshtml')
am
Out[ ]:

雲がかかっている時間帯もありますが、潮目が時間で変化していることが確認できます。

まとめ


この章では、「しきさい」と「ひまわり」の海面水温データについて、主に可視化を通してPythonの様々なライブラリの使い方を学習しました。この教材で紹介した以外のデータを表示してみたり、温度の範囲・拡大する場所を変更することで、同じ種類のデータでも違った見え方になってくると思います。ぜひ挑戦してみてください。

コラム3:漁場はどこにある?


冒頭で潮目は好漁場になりやすいことを説明しましたが、実際にどの場所で漁が行われたかを調べるのは難しいです。

Global Fishing WatchというWebサービスでは、船舶に搭載している自動船舶識別装置(Automatic Identification System, AIS)から取得した位置情報や衛星で観測した夜間光のデータ等から船の動きを見ることができますが、漁を行った場所を示すものではありません。

直接漁場を探すのは難しいですが、各漁港の水揚量を見てみると、漁がうまくいったかの大まかな傾向は見ることができます。

例えば、今回取り上げた岩手県の漁港の水揚量は、いわて大漁ナビ の「市況日報」で過去分のデータが閲覧できます。潮目が発生した時期と魚種別の水揚量の関係などを調べてみると、面白い結果が得られるかもしれませんね。