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

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

$ python3 get-and-plot-himawari.py

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

#------------------------------------------------------------------------------
#  千葉大学リモートセンシング研究センターから取得する衛星画像の設定    
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)