Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
11
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

気象×Python 〜全国合成レーダー〜

気象レーダーによるエコー強度・エコー頂高度を描画していきたいと、思います。

1. 全国合成レーダーについて

▶気象庁保有の全国20台の気象レーダーで観測した10分間隔のレーダーエコー強度・エコー頂高度が京都大学生存圏研究所より提供されている。
▶10分間隔の降水強度や雨雲の高さから降水システムを把握できる。

エコー強度 エコー頂高度
キャプチャ.PNG キャプチャ1.PNG

デジタル台風@国立情報学研究所より引用。

データセット詳細

エコー強度 エコー頂高度
1km.PNG0.0125°×0.0083°(1km格子) 2.5km.PNG0.03125°×0.025°(2.5km格子)

2. データ取得方法

2.1 京都大学生存圏研究所からtar fileのダウンロード

今回自分はwgetコマンドを用いていますが、データ取得はお好きなやり方でどーぞ(wgetコマンドはインストールする必要あり)。

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. 全国合成レーダー描画

あとは切り出したデータを地図上で表現するだけですが、ここで注意したいのが生データを降水強度やエコー頂高度に変換しなければいけないことです。以下が変換のルックアップテーブルで、生データがデータ代表値に対応します。

エコー強度 エコー頂高度
1kmlut.PNG 2.5kmlut.PNG

つまり、

データ代表値⇒降水強度
・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()
エコー強度 エコー頂高度
Figure_1.png Figure_2.png
キャプチャ.PNG キャプチャ1.PNG

最初に示したレーダー画像とだいたい一致してそうだったので正しく読み出せたと思います。ありがとうございました。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
11
Help us understand the problem. What are the problem?