気象レーダーによるエコー強度・エコー頂高度を描画していきたいと、思います。
#1. 全国合成レーダーについて
▶気象庁保有の全国20台の気象レーダーで観測した10分間隔のレーダーエコー強度・エコー頂高度が京都大学生存圏研究所より提供されている。
▶10分間隔の降水強度や雨雲の高さから降水システムを把握できる。
| エコー強度 | エコー頂高度 |
 |
 |
| エコー強度 | エコー頂高度 |
 0.0125°×0.0083°(1km格子)
|
 0.03125°×0.025°(2.5km格子)
|
radar_wget.py
import subprocess
import os
import glob
host = "http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synthetic/original"
#10分間隔のリスト
t_list = ["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20","21","22","23"]
#取得したい年月日を入力
yyyy = "××××"
mm = "××"
dd = "××"
for t in t_list:
time = [f"{t}0",f"{t}1",f"{t}2",f"{t}3",f"{t}4",f"{t}5"]
for i in range(0,6):
if os.path.exists(f"./{yyyy}{mm}/{dd}"):
if os.path.exists(f"./{yyyy}{mm}/{dd}/Z__C_RJTD_{yyyy}{mm}{dd}{time[i]}000_RDR_JMAGPV__grib2.tar"):
print(f"{yyyy}{mm}{dd}{time[i]}000 already exists")
else:
url = f"{host}/{yyyy}/{mm}/{dd}/Z__C_RJTD_{yyyy}{mm}{dd}{time[i]}000_RDR_JMAGPV__grib2.tar"
subprocess.run(f"wget {url} -P ./{yyyy}{mm}/{dd}/", shell=True)
else:
os.makedirs(f"./{yyyy}{mm}/{dd}")
url = f"{host}/{yyyy}/{mm}/{dd}/Z__C_RJTD_{yyyy}{mm}{dd}{time[i]}000_RDR_JMAGPV__grib2.tar"
subprocess.run("wget {}".format(url), shell=True)
##2.2 wgrib2によるデータ切り出し
取得したデータを展開すると、GRIB2形式として出力されるので、wgrib2コマンドを使ってbinデータを切り出します。
こちらの記事を参考にさせて頂きました。（https://qiita.com/kinosi/items/6f4f754f3f88f1fd74b4）
# wgrib2 (元のGRIB2データ) [オプション] (出力ファイル) ← オプションによってフォーマットを指定できる。
$ wgrib2 Z__C_RJTD_**************_RDR_JMAGPV__grib2.bin -order we:ns -no_header -bin hoge.bin
※pygribというgrib形式を読み込むモジュールでもやってみたが、今回のレーダーデータは対応していないと思われる。できたら方がいらっしゃいましたら教えてください笑
#3. 全国合成レーダー描画
あとは切り出したデータを地図上で表現するだけですが、ここで注意したいのが生データを降水強度やエコー頂高度に変換しなければいけないことです。以下が変換のルックアップテーブルで、生データがデータ代表値に対応します。
| エコー強度 | エコー頂高度 |
|
|
データ代表値⇒降水強度
・0⇒0 mm/h
・0.1⇒0.1 mm/h
・0.25⇒0.2 mm/h
：
・260⇒256 mm/h
みたいな感じで変換していきたいと思います。
まずエコー強度の描画プログラムがこちら。
echo_intensity.py
# -*- coding: utf-8 -*-¥
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
#↓ルックアップテーブル作成
def lookup_table(start1, end1, step1, start2, end2, step2):
keys = []
values = []
for i in np.arange(start1, end1, step1):
i = str(round(i, 2))
keys.append(i)
for j in np.arange(start2, end2, step2):
j = round(j, 2)
values.append(j)
d = dict(zip(keys, values))
return d
class Radar_echo:
def __init__(self):
f = open("Z__C_RJTD_20141216000000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_.bin", mode="rb") #切り出したデータを読み込んでから配列変換
echo = np.fromfile(f, dtype="float32",sep="").reshape(3360,2560)
#とりあえず0,0.1,260のデータ代表値を降水強度へ変換
echo = np.where(echo < 0, 0, echo)
echo = np.where(echo == 0.1, 0.1, echo)
self.echo = np.where(echo == 260, 256, echo)
def lut_to_echo(self):
d1 = lookup_table(0.25, 2.05, 0.1,0.2, 2.0, 0.1)
d2 = lookup_table(2.13, 5.13, 0.25, 2.0, 4.95, 0.25)
d3 = lookup_table(5.25, 10.25, 0.5, 5.0, 10.0, 0.5)
d4 = lookup_table(10.5, 180.5, 1.0, 10, 180, 1.0)
d5 = lookup_table(181, 257, 2.0, 180, 256, 2.0)
d1.update(d2)
d1.update(d3)
d1.update(d4)
d1.update(d5)
#ルックアップテーブルをもとに降水強度へ変換
for k, v in d1.items():
self.echo = np.where(self.echo == float(k), v, self.echo)
def draw(self):
#東西・南北のグリッドを作成
lon = np.round(np.linspace(118, 150, 2560), 4)
lat = np.round(np.linspace(20, 48, 3360)[::-1], 4)
X, Y = np.meshgrid(lon,lat)
fig = plt.figure(dpi=100)
levels = [0.5, 1, 3, 5, 8, 10]
cmap = plt.get_cmap('jet')
cmap.set_under("w")
cmap.set_over("firebrick")
#Basemapによる描画
m = Basemap(projection="cyl", resolution="i", llcrnrlat=20,urcrnrlat=50, llcrnrlon=120, urcrnrlon=150)
m.drawcoastlines(color='black')
m.drawmeridians(np.arange(0, 360, 5), labels=[True, False, False, True],linewidth=0.0)
m.drawparallels(np.arange(-90, 90, 5), labels=[True, False, True, False],linewidth=0.0)
im = plt.contourf(X,Y, self.echo, levels, cmap=cmap, extend='both')
cb = m.colorbar(im, "right", size="2.5%")
cb.set_label('mm/hr')
plt.title("Radar Echo Intensity 9:00JST 16DEC2014")
plt.show()
if __name__ == '__main__':
echo = Radar_echo()
echo.lut_to_echo()
echo.draw()
続いてエコー頂高度の描画プログラムがこちら。
echo_height.py
# -*- coding: utf-8 -*-¥
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
class Radar_etop:
def __init__(self):
f = open("Z__C_RJTD_20141216000000_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_.bin", mode="rb") #切り出したデータを読み込んでから配列変換
etop = np.fromfile(f, dtype="float32",sep="").reshape(1120,1024)
etop = np.where(etop < 0, 0, etop)
etop = np.where(etop == 1, 1, etop)
def lut_to_etop(self):
keys = []
values = []
for i in np.arange(3.0, 17.0, 2.0):
i = round(i, 1)
keys.append(i)
for j in np.arange(2.0, 16.0, 2.0):
j = round(j,1)
values.append(j)
d = dict(zip(keys, values))
for k, v in d.items():
self.etop = np.where(self.etop == k, v, self.etop)
def draw(self):
#東西・南北のグリッドを作成
lon = np.round(np.linspace(118, 150, 1024), 5)
lat = np.round(np.linspace(20, 48, 1120)[::-1], 3)
X, Y = np.meshgrid(lon,lat)
fig = plt.figure(dpi=100)
interval = list(np.arange(2, 10, 2))
cmap = plt.get_cmap('jet')
cmap.set_under("w")
cmap.set_over("firebrick")
#Basemapによる描画
m = Basemap(projection="cyl", resolution="i", llcrnrlat=20,urcrnrlat=50, llcrnrlon=120, urcrnrlon=150)
m.drawcoastlines(color='black')
m.drawmeridians(np.arange(0, 360, 5), labels=[True, False, False, True],linewidth=0.0)
m.drawparallels(np.arange(-90, 90, 5), labels=[True, False, True, False],linewidth=0.0)
im = plt.contourf(X,Y, self.etop, interval, cmap=cmap, extend='both')
cb = m.colorbar(im, "right", size="2.5%")
cb.set_label('km')
plt.title("Radar Echo Height 9:00JST 16DEC2014")
plt.show()
if __name__ == '__main__':
etop = Radar_etop()
etop.lut_to_etop()
etop.draw()
| エコー強度 | エコー頂高度 |
|
|
|
|
最初に示したレーダー画像とだいたい一致してそうだったので正しく読み出せたと思います。ありがとうございました。