気象×Python 〜気象衛星ひまわり8号を用いたRGB合成〜

気象衛星ひまわりグリッドデータの描画です。とりあえず自分の備忘録として、ゆるく書きます。
#1.　静止気象衛星ひまわり８号の概要

▶2014年10月7日に打ち上げられ、2015年7月7日に運用を開始した。世界最先端の観測能力を有する可視赤外放射計（AHI：Advanced Himawari Imager）を搭載し、可視域〜赤外域の16バンドを記録することができる。（気象庁より）
▶可視域の赤・青・緑それぞれに感度を持つバンドの画像を合成することで、強調して表示したい現象を分かりやすく表示することができるようになった。https://www.jma-net.go.jp/sat/himawari/satobs.html#rgb-composite

#2.　使用データ
▶千葉大学環境リモートセンシング研究センター「ひまわり8/9号」の全球スキャン (FD) の gridded data (緯度経度直行座標系による精密幾何補正済データ)→詳細はこちらのページ
▶ちなみに今回選んだ事例は東日本に記録的な大雨をもたらした2019/10/12の台風19号。
#2.　RGB画像作成コード
こちらの資料を参考に実装しました。
あと、生データを輝度温度，反射率へ変換するためのルックアップテーブルは必要です。
ftp://hmwr829gr.cr.chiba-u.ac.jp/gridded/FD/support/
###（１）Airmass RGB
####画像の特徴と画像作成レシピ
気団の解析、ジェット気流の位置、上層トラフや上層渦、大気沈降域などの推定に用いられる。

色	赤	緑	青
利用するバンドの種類	バンド8 - バンド10	バンド12 - バンド13	バンド8 (反転)
階調の設定	-25～0K	-40～5K	243～208K
ガンマ値	1.0	1.0	1.0

airmass.py

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap


#最小値0、最大値1にする数量正規化メソッド
def normalization(x):
    x_min = x.min()
    x_max = x.max()
    x_norm = (x - x_min) / (x_max - x_min)
    return x_norm   #正規化した配列を返す

#階調反転メソッド
def reverse(color, min, max):
  for x in color[:]:
    y = ((0-1)/(max-min))*x + (max/(max-min))
    color = np.where(color == x, y, color)
  return color

class Himawari:
  def __init__(self):
    #必要なバンドの読み込み
    band08 = np.fromfile("201910120600.tir.06.fld.geoss", dtype='>u2').reshape(6000, 6000)
    band10 = np.fromfile("201910120600.tir.08.fld.geoss", dtype='>u2').reshape(6000, 6000)
    band12 = np.fromfile("201910120600.tir.10.fld.geoss", dtype='>u2').reshape(6000, 6000)
    band13 = np.fromfile("201910120600.tir.01.fld.geoss", dtype='>u2').reshape(6000, 6000)

    #欠損値置換
    band08 = np.where(band08 > 2000, 2000, band08)
    band10 = np.where(band10 > 4050, 4050, band10)
    band12 = np.where(band12 > 4050, 4050, band12)
    band13 = np.where(band13 > 4050, 4050, band13)

    #ルックアップテーブル参照
    #日本領域にトレミング(50→30, l25→150)
    band08 = np.loadtxt('count2tbb/tir.06',usecols=(1,))[band08][500:1500, 2000:3250]
    band10 = np.loadtxt('count2tbb/tir.08',usecols=(1,))[band10][500:1500, 2000:3250]
    band12 = np.loadtxt('count2tbb/tir.10',usecols=(1,))[band12][500:1500, 2000:3250]
    band13 = np.loadtxt('count2tbb/tir.01', usecols=(1,))[band13][500:1500, 2000:3250]

    #差分をとって階調域の限定
    self.R = (band08 - band10).clip(-25, 0)
    self.G = (band12 - band13).clip(-40, 5)
    self.B = band08.clip(208, 243)

  def rgb_synthesis(self):
    R = normalization(self.R)
    G = normalization(self.G)
    B = reverse(self.B, 208, 243)
    
    rgb = np.dstack((R, G, B))
    self.rgb = np.power(rgb, 1.0)

  def draw(self):
    fig = plt.figure(dpi=100)
    fname1 = 'gadm/gadm36_JPN_1'
    m = Basemap(projection="cyl", resolution="i", llcrnrlat=30, urcrnrlat=50, llcrnrlon=125, urcrnrlon=150, area_thresh=100, fix_aspect=True)
    m.imshow(self.rgb, origin="upper")
    m.drawcoastlines(color='k')
    m.drawmeridians(np.arange(0, 360, 10), labels=[True, False, False, True], linewidth=0)
    m.drawparallels(np.arange(-90, 90, 10), labels=[ True,False, True, False],linewidth=0)
    m.readshapefile(fname1, 'prefectural_bound1', color='k', linewidth=.8) #県境
    plt.title("Airmass RGB 20191012 15:00JST")
    plt.savefig('airmass_201910121500.jpeg', dpi=100)


if __name__ == '__main__':
  h = Himawari()
  h.rgb_synthesis()
  h.draw()

それっぽいのはできた、、

ちなみに気象衛星センターによると、暖気側の温かい気団は緑系色、寒気側の冷たい気団は青紫系色、トラフ後面の大気沈降域は赤系色で表されるとのこと。
###（2）Day convective storm RGB
####画像の特徴と画像作成レシピ
突風や竜巻等のシビア現象を伴う対流雲の判別に用いられる。

色	赤	緑	青
利用するバンドの種類	バンド8 - バンド10	バンド7 - バンド13	バンド5 - バンド3
階調の設定	-35～5K	-5～60K	-75～25%
ガンマ値	1.0	1.0	1.0
コードは省略しますが結果だけ。シビア現象を伴う積乱雲は黄色表示となるみたい。

Day convective storm & Radar	地上天気図

天気図はtenki.jpより引用してます。

データ容量は大きいので注意です。
以上！！！