【Python】気象衛星ひまわり画像データの取得と描画【Cartopy】

【Python】気象衛星ひまわり画像データの取得と描画【Cartopy】

前回の記事では,PythonでGISデータの描画が簡単にできるCartopyのインストールと,Cartopyを使って日本付近の地図を描く方法を紹介しました。今回は,気象衛星ひまわりが撮影した雲画像(赤外線画像)を地図上に描いて,日本上空の雲の様子を見てみましょう。

なお,今回使用する気象衛星ひまわりが撮影した画像データは,千葉大学環境リモートセンシング研究センターが,非営利目的での使用や出版物等の掲載に限り無償で提供しているものです。ご利用にあたっては,注意事項をよく読んで活用するようにしましょう。

動作検証済み環境

macOS Monterey 12.6.2, Python 3.9.15

気象衛星ひまわりが撮影した画像データの取得と描画

気象衛星ひまわりが撮影した画像データは,千葉大学リモートセンシング研究センター (CEReS) のサイトから,FTPで取得できます。FTP でファイルにアクセスしデータを取得するには Python 標準ライブラリの ftp が使用できます。

取得した画像データ(バイナリ)は bzip2 という形式で圧縮されています。したがって,解凍し,バイナリデータを読み込む必要があります。これには,Python 標準ライブラリ bz2 が使用できます。

気象衛星は様々な波長領域に感度を持つセンサーをもっています。今回は,バンド13と呼ばれる中心波長 10.4 μm の赤外線領域の画像データを使用します。夜でも見ることができるのでテレビ等でよく見る雲画像は,赤外線画像であることが多いと思います。詳しい特性については,気象衛星センターの説明がわかりやすいです。

画像データの緯度・経度,分解能等の情報は,千葉大学リモートセンシング研究センターのリリースノートで公開されています。下の表にデータ名の名前とバンドの対応,解像度の表を示します。これは,千葉大学リモートセンシング研究センターのリリースノート (2022/12/27閲覧)から転載したものです。

CEReS griddedひまわり8/9号バンドPixel x Line空間解像度
EXT 01Band03 (0.64 μm)24000 x 240000.005 degree (500 m 相当)
VIS 01Band01 (0.47 μm)12000 x 120000.01 degree (1 km 相当)
VIS 02Band02 (0.51 μm)12000 x 120000.01 degree (1 km 相当)
VIS 03Band04 (0.86 μm)12000 x 120000.01 degree (1 km 相当)
SIR 01Band05 (1.6 μm)6000 x 60000.02 degree (2 km 相当)
SIR 02Band06 (2.3 μm)6000 x 60000.02 degree (2 km 相当)
TIR 01Band13 (10.4 μm)6000 x 60000.02 degree (2 km 相当)
TIR 02Band14 (11.2 μm)6000 x 60000.02 degree (2 km 相当)
TIR 03Band15 (12.4 μm)6000 x 60000.02 degree (2 km 相当)
TIR 04Band16 (13.3 μm)6000 x 60000.02 degree (2 km 相当)
TIR 05Band07 (3.9 μm)6000 x 60000.02 degree (2 km 相当)
TIR 06Band08 (6.2 μm)6000 x 60000.02 degree (2 km 相当)
TIR 07Band09 (6.9 μm)6000 x 60000.02 degree (2 km 相当)
TIR 08Band10 (7.3 μm)6000 x 60000.02 degree (2 km 相当)
TIR 09Band11 (8.8 μm)6000 x 60000.02 degree (2 km 相当)
TIR 10Band12 (9.6 μm)6000 x 60000.02 degree (2 km 相当)
表1:千葉大学リモートセンシング研究センターが提供するグリッドデータとひまわり8/9号バンドの対応関係と解像度

緯度経度範囲東経 85度 – 西経155度 (85E – 155W (205E)), 北緯60度 – 南緯60度 (60N – 60S)
格納データ仕様ヘッダー無し2 byte 符号無し整数 (unsigned short), big endian data order のバイナリデータ
データ書き出し順西 -> 東 (左 -> 右),北 -> 南(上 -> 下)に書き出し
格納データ気象庁より配信されたひまわりスタンダード(HS) データ内のカウント値そのもの 
幾何補正で値が入っていないところには例外値として65535 が入る点に注意が必要
表2:画像データの仕様

画像データそのものの値は,デジタルナンバー値 (DN値) とよばれます。DN値を物理的に,より解釈しやすい情報にするための対応表が用意されています。python標準ライブラリのurllibを使って,httpでアクセスし,この対応表をダウンロードしてみます。そして,これを使って,DN値を黒体放射輝度温度になおしてCartopyを使ってデータを描画してみましょう。

実装方法

以上述べたことをPythonでコーディングしてみたいと思います。以下に実装例を示します。このコードを好きな名前 (ここでは,get-and-plot-himawari.py とします) をつけて,好きなディレクトリ (ここでは,~/Desktop/LabCode/python/himawari とします) に保存します。

import bz2
import os
import ftplib
import tarfile
import urllib.request

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np


#------------------------------------------------------------------------------
#  千葉大学リモートセンシング研究センターから取得する衛星画像の設定    
date = "20220917" # 8桁の数字で入力する
time = "0000" # 協定世界時刻であることに注意。つまり,日本時間はこれよりも9時間進んでいる
band = "tir" # バンド名 vis, tir 
ch = "01" # チャンネル番号

# 最大,最小の緯度と経度,縦横の画素数
lon_min = 85
lon_max = 205
lat_min = -60
lat_max = 60
n = 6000

# ひまわり画像データの範囲
extent_hmwr = (lon_min, lon_max, lat_min, lat_max)

#------------------------------------------------------------------------------
# ひまわりの画像データを読み込む
# ひまわりデータは千葉大学リモートセンシング研究センターがFTPで提供しているものを利用する。
server = "hmwr829gr.cr.chiba-u.ac.jp"
tail = "fld.geoss.bz2"
fname_hmwr = f"{date}{time}.{band}.{ch}.{tail}"
url = f"{server}/{date[0:6]}/{band.upper()}/{fname_hmwr}"

if not os.path.isfile(fname_hmwr):
    print('*** 画像ファイルをダウンロードします: ', fname_hmwr,)
    ftp = ftplib.FTP(server)
    ftp.login('', '') # anonymousユーザーでログイン
    ftp.cwd(f'gridded/FD/V20190123/{date[0:6]}/{band.upper()}')

    with open(fname_hmwr ,'wb') as fhandle:
        ftp.retrbinary(f'RETR {fname_hmwr}', fhandle.write)
    print('*** 画像ファイルのダウンロードが完了しました: ', fname_hmwr,)

    ftp.quit()


#------------------------------------------------------------------------------
# センサのカウント値を放射輝度へ変換する必要がある。
# そのための変換テーブルをダウンロードする
fname_cnt2tbb = 'count2tbb_v102.tgz'
url = f'http://www.cr.chiba-u.jp/databases/GEO/H8_9/FD/{fname_cnt2tbb}'

# ファイルのダウンロード
if not os.path.isfile(fname_cnt2tbb):
    print("*** 変換テーブルをダウンロードしています:", fname_cnt2tbb)
    urllib.request.urlretrieve(url, fname_cnt2tbb)

# ファイルの解凍
if not os.path.isdir(f'count2tbb_v102'):
    print("*** 変換テーブルを解凍しています")
    with tarfile.open(fname_cnt2tbb, "r:gz") as tarin:
        tarin.extractall()

# カウント値から黒体放射輝度温度(TBB)への変換テーブルの読み込み
print("*** 変換テーブルを読み込みます")
_, cnt2tbb = np.loadtxt(f'count2tbb_v102/{band}.{ch}', unpack=True)

#  標準ライブラリbz2を使用して,ひまわり画像を開く。
with bz2.open(fname_hmwr) as bz2fin:
    buf_cnt = bz2fin.read()
    data_cnt = np.frombuffer(buf_cnt, dtype='>u2').reshape(n, n)
    
    data_tbb = cnt2tbb[data_cnt] # 画像データを変換テーブルで変換する


#------------------------------------------------------------------------------
# プロット
fig = plt.figure(figsize=(294/25.4, 210/25.4))
out_fig_path = f'{date}{time}-himawari_{band}{ch}.png'

fig.text(0.10, 0.97, f"{date} {time}UTC", fontsize=14)
fig.text(0.10, 0.94, f"BAND:{band}, CH:{ch}", fontsize=14)
#  図法の設定 (適宜コメントアウトを外して使ってください)
# mapcrs = ccrs.PlateCarree()
mapcrs = ccrs.NorthPolarStereo(central_longitude=140.0, true_scale_latitude=60.0) # 気象庁「日々の天気図」に近いもの
# mapcrs = ccrs.Robinson(central_longitude=180.0) # ロビンソン図法
# mapcrs = ccrs.Mercator() # メルカトル図法
# mapcrs = ccrs.Mollweide() # モルワイデ図法

ax = fig.add_axes([0.10, 0.10, 0.75, 0.75], projection=mapcrs)


# 海岸線や国境を描画する
ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color="green")
ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', color="green")

# 地図を描く緯度経度の範囲の設定
ax.set_extent([117.0, 156.0, 18.0, 61.0], ccrs.PlateCarree())

# ひまわり衛星画像の描画
ax.imshow(data_tbb, origin='upper', extent=extent_hmwr, interpolation='none',cmap="gray_r", transform=ccrs.PlateCarree())

ax.gridlines(xlocs=np.arange(100, 180, 10), 
             ylocs=np.arange(10, 90, 10),
             linestyle='dashed', color='green', linewidth=0.5)

#------------------------------------------------------------------------------
# 図の保存
plt.savefig(out_fig_path, transparent=False)

プログラムを実行する

ターミナルを開き,

$ cd Desktop/LabCode/python/himawari

と入力し,先程プログラムを保存したディレクトリを移動します。あとは以下のコマンドでget-and-plot-himawari.pyを実行するだけです。( $マークは「プロンプト」です。入力する必要はありません。)

$ python3 get-and-plot-himawari.py

実行結果

~/Desktop/LabCode/python/himawari に,まず,202209170000.tir.01.fld.geoss.bz2 がダウンロードされて,次いでcount2tbb_v102.tgz がダウンロードされます。その後,解凍されてcount2tbbというフォルダができます。最後に,202209170000-himawari_tir01.png というファイルができて,以下のようになっていれば成功です!

生成される図

上のプログラムで描画される図は,2022年に史上最強クラスで日本に接近した台風14号を捉えた画像です。2022年9月17日,日本時間9時の段階では,九州の南の海上にあって,台風の眼がくっきり見えます。

気象庁ホームページ「日々の天気図」より令和4年9月17日の天気図と記事を抜粋。

参考のために,気象庁ホームページより「日々の天気図」から2022年9月17日,日本時間9時の地上天気図を転載しておきます。台風の中心位置に眼があることや,北緯50度付近の雲が低気圧に対応していることがわかります。

コードの解説

上に書いたソースコードの解説をしていきます。

import bz2
import ftplib
import os
import tarfile
import urllib.request

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np

まず,必要なモジュールをインポートしておきます。上段はPython標準ライブラリです。下段はpipなどで入れたライブラリです。

#------------------------------------------------------------------------------
#  千葉大学リモートセンシング研究センターから取得する衛星画像の設定    
date = "20220917" # 8桁の数字で入力する
time = "0000" # 協定世界時刻であることに注意。つまり,日本時間はこれよりも9時間進んでいる
band = "tir" # バンド名 vis, tir 
ch = "01" # チャンネル番号

# 最大,最小の緯度と経度,縦横の画素数
lon_min = 85
lon_max = 205
lat_min = -60
lat_max = 60
n = 6000

# ひまわり画像データの範囲
extent_hmwr = (lon_min, lon_max, lat_min, lat_max)

千葉大学リモートセンシング研究センターからデータを取得しますが,データの日付や時刻,バンド名,チャンネル番号などの設定をします。

日付と時刻は協定世界時であることに注意してください。つまり,日本時間で9時なら,time は 0000 としなければなりません。

緯度と経度の範囲,縦横の画素数はデータの仕様に合わせて設定します。これも,千葉大学リモートセンシング研究センターのリリースノートに従って設定します。

#------------------------------------------------------------------------------
# ひまわりの画像データを読み込む
# ひまわりデータは千葉大学リモートセンシング研究センターがFTPで提供しているものを利用する。
server = "hmwr829gr.cr.chiba-u.ac.jp"
tail = "fld.geoss.bz2"
fname_hmwr = f"{date}{time}.{band}.{ch}.{tail}"
url = f"{server}/{date[0:6]}/{band.upper()}/{fname_hmwr}"

if not os.path.isfile(fname_hmwr):
    print('*** 画像ファイルをダウンロードします: ', fname_hmwr,)
    ftp = ftplib.FTP(server)
    ftp.login('', '') # anonymousユーザーでログイン
    ftp.cwd(f'gridded/FD/V20190123/{date[0:6]}/{band.upper()}')

    with open(fname_hmwr ,'wb') as fhandle:
        ftp.retrbinary(f'RETR {fname_hmwr}', fhandle.write)
    print('*** 画像ファイルのダウンロードが完了しました: ', fname_hmwr,)

    ftp.quit()

まず,ひまわり画像データをFTPで取得します。Python標準ライブラリであるftplibを使用します。無駄なアクセスを避けるために,ファイルがある場合にはこの部分はスキップします。

#------------------------------------------------------------------------------
# センサのカウント値を放射輝度へ変換する必要がある。
# そのための変換テーブルをダウンロードする
fname_cnt2tbb = 'count2tbb_v102.tgz'
url = f'http://www.cr.chiba-u.jp/databases/GEO/H8_9/FD/{fname_cnt2tbb}'

# ファイルのダウンロード
if not os.path.isfile(fname_cnt2tbb):
    print("*** 変換テーブルをダウンロードしています:", fname_cnt2tbb)
    urllib.request.urlretrieve(url, fname_cnt2tbb)

# ファイルの解凍
if not os.path.isdir(f'count2tbb_v102'):
    print("*** 変換テーブルを解凍しています")
    with tarfile.open(fname_cnt2tbb, "r:gz") as tarin:
        tarin.extractall()

# カウント値から黒体放射輝度温度(TBB)への変換テーブルの読み込み
print("*** 変換テーブルを読み込みます")
_, cnt2tbb = np.loadtxt(f'count2tbb_v102/{band}.{ch}', unpack=True)

この部分は,センサのカウント値から放射輝度温度に変換するテーブルをダウンロードする部分です。httpでアクセスできるので,urllibをつかって取得します。この部分も,無駄なアクセスを避けるために,すでにある場合は実行しません。ファイルは tar で圧縮されているので,tarfile を使用して現在のディレクトリ上に解凍します。解凍したら,バンドとチャンネルに対応する変換テーブルを読み込みます。

#  標準ライブラリbz2を使用して,ひまわり画像を開く。
with bz2.open(fname_hmwr) as bz2fin:
    buf_cnt = bz2fin.read()
    data_cnt = np.frombuffer(buf_cnt, dtype='>u2').reshape(n, n)
    
    data_tbb = cnt2tbb[data_cnt] # 画像データを変換テーブルで変換する

この部分は,bzip2で圧縮されているひまわり画像データを開いて,data_cnt という二次元のndarray に格納している部分です。最終部分では,変換テーブルで輝度温度に変換しています。

#------------------------------------------------------------------------------
# プロット
fig = plt.figure(figsize=(294/25.4, 210/25.4))
out_fig_path = f'{date}{time}-himawari_{band}{ch}.png'

fig.text(0.10, 0.97, f"{date} {time}UTC", fontsize=14)
fig.text(0.10, 0.94, f"BAND:{band}, CH:{ch}", fontsize=14)
#  図法の設定 (適宜コメントアウトを外して使ってください)
# mapcrs = ccrs.PlateCarree()
mapcrs = ccrs.NorthPolarStereo(central_longitude=140.0, true_scale_latitude=60.0) # 気象庁「日々の天気図」に近いもの
# mapcrs = ccrs.Robinson(central_longitude=180.0) # ロビンソン図法
# mapcrs = ccrs.Mercator() # メルカトル図法
# mapcrs = ccrs.Mollweide() # モルワイデ図法

ax = fig.add_axes([0.10, 0.10, 0.75, 0.75], projection=mapcrs)


# 海岸線や国境を描画する
ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color="green")
ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', color="green")

# 地図を描く緯度経度の範囲の設定
ax.set_extent([117.0, 156.0, 18.0, 61.0], ccrs.PlateCarree())

# ひまわり衛星画像の描画
ax.imshow(data_tbb, origin='upper', extent=extent_hmwr, interpolation='none',cmap="gray_r", transform=ccrs.PlateCarree())

ax.gridlines(xlocs=np.arange(100, 180, 10), 
             ylocs=np.arange(10, 90, 10),
             linestyle='dashed', color='green', linewidth=0.5)

#------------------------------------------------------------------------------
# 図の保存
plt.savefig(out_fig_path, transparent=False)

この部分はCartopyで読み込んだひまわり画像データを描画している部分です。ファイル名や図の上部に日付やバンド,チャンネル番号をつけるようにしました。

気象庁ホームページの「日々の天気図」掲載の天気図と比較がしやすいように,ポーラーステレオ図法でプロットしています(ちなみに,気象庁の資料によると,「北緯 60度、東経 140 度を基準としたポーラーステレオ図法で投影している」とのことですが,微妙にずれます… )。

最後に

以上のように,Pythonをつかって気象衛星が撮影した日本付近の画像をプロットすることができました。今回は,赤外線画像を使って雲の様子を眺めてみただけですが,バンドの組み合わせ等で有用な情報が得られます。また,画像処理のテクニックを用いることで,台風の位置をトラッキングしたりすることも可能でしょう。

次回は,ひまわり画像データと数値予報データの重ね書きの方法を紹介したいと思います。

最後に,このような価値の高い画像を (非営利または,出版物での利用に限るとはいえ) 無料で得られるようにしてくださった,気象衛星の開発,維持,そして画像データ公開に携わっている方々に感謝を申し上げます。

参考にしたサイト・資料

この記事を書くにあたって,以下のウェブサイトや資料を参考にしました。このような有益な情報や資料を無料で公開しているウェブサイトの管理者や資料を作成された方に感謝を申し上げます。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です