気象庁全国合成レーダーGPVの描画(python)

気象レーダーによる降水量メッシュデータを、pythonでダウンロード、描画していきます。(Ubuntu-20.04 on WSL2)

気象庁全国合成レーダーGPVとは

気象庁は全国に20台の気象レーダを保有しており、10分間隔で雨の観測が行われています。1つ1つのレーダは半径約数100kmの範囲の観測に限られますが、これらの多数のレーダからの情報をひとまとめにしたメッシュ状のマップデータが提供されています。これを全国合成レーダGPVといいます。GPVは「Grid Point Value (格子点値)」の略で、メッシュ状のデータという意味です。 

気象庁のレーダー配置図 令和3年7月現在 (気象庁HPより引用)

・エコー強度(1kmメッシュ)

・エコー頂高度(2.5kmメッシュ)

の2種類のデータが存在していますが、ここではエコー強度の方を描画していきましょう。

気象庁HP 気象レーダ解説ページ https://www.jma.go.jp/jma/kishou/know/radar/kaisetsu.html

全国合成レーダーGPV https://www.data.jma.go.jp/add/suishin/catalogue/format/ObdObs001_format.pdf

データのダウンロード・変換

京都大学・生存圏研究所(http://database.rish.kyoto-u.ac.jp/arch/jmadata/)よりデータを取得します。

(教育研究機関向けにデータが公開されています。「企業活動等のためにデータを頻繁に必要とされる方は,気象業務支援センターからデータを直接購入し,データ提供スキーム全体の維持発展にご協力ください.」とのこと。)

以下は、テータをダウンロード・読み取りやすい形式に変換するプログラムです。wgetとwgrib2を事前にインストールしておく必要があります。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import subprocess
import os
import pandas as pd

def download_time(time):
    ## 気象庁全国合成レーダーGPVのダウンロード

    ## 生存圏研究所ダウンロード元サイト
    http  = "http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synthetic/original"
    
    ## 保存先ディレクトリの指定
    Opath = "./data"
    day_dir = time.strftime("%Y/%m/%d")
    basename = "Z__C_RJTD_{}00_RDR_JMAGPV__grib2.tar".format(time.strftime("%Y%m%d%H%M"))
    fname    = "{}/{}/{}".format(Opath, day_dir, basename)

    ## すでにファイルが存在しなければ、ダウンロードを行う
    if os.path.exists(fname):
        print("{} Already exists.".format(basename))
    else:
        os.makedirs("{}/{}".format(Opath, day_dir), exist_ok=True)
        url      = "{}/{}/{}".format(http,  day_dir, basename)

        ## wgetコマンドでデータのダウンロード
        subprocess.run("wget {} -P {}/{}/".format(url, Opath, day_dir), shell=True)
        ## tarコマンドでダウンロードした圧縮ファイルの解凍
        subprocess.run("tar -xvf {} -C {}/{}/".format(fname, Opath, day_dir), shell=True)
    
    ## wgrib2コマンドでgrib2データをバイナリファイルに変換
    outfile  = "{}/{}/Z__C_RJTD_{}00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL.dat".format(Opath, day_dir, time.strftime("%Y%m%d%H%M"))
    if not os.path.exists(outfile):
        GgisFile = "{}/{}/Z__C_RJTD_{}00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin".format(Opath, day_dir, time.strftime("%Y%m%d%H%M"))
        subprocess.run("wgrib2 {} -d 1 -no_header -bin {}".format(GgisFile, outfile), shell=True)

    return

if __name__=="__main__":
    ## ダウンロードしたいデータの日時(10分間隔)を指定
    time = pd.Timestamp(2021,2,15,12,0)
    download_time(time)

上のコードの解説

wgetコマンドで生存圏研究所のURLからデータをダウンロードします。

ダウンロードしたデータは

Z__C_RJTD_20210215120000_RDR_JMAGPV__grib2.tar

のような圧縮ファイルになっているので、

tar -xvf Z__C_RJTD_20210215120000_RDR_JMAGPV__grib2.tar

で展開します。すると、

Z__C_RJTD_20210215120000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin

Z__C_RJTD_20210215120000_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_grib2.bin

という2つのファイルができます。上が1kmメッシュのエコー強度、下が2.5kmメッシュのエコー頂高度のデータに対応しています。ここでは上を使います。

これらはgrib2というデータ形式になっており、このままでは読み込むことができないので、wgrib2コマンドで通常のバイナリファイルへと変換します。

wgrib2 Z__C_RJTD_20210215120000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin -d 1 -no_header -bin 出力ファイル名

上のコードの実行

実行します。

$ ./download_RDR_JMAGPV.py
--2021-07-11 22:02:29--  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synthetic/original/2021/02/15/Z__C_RJTD_20210215120000_RDR_JMAGPV__grib2.tar
database.rish.kyoto-u.ac.jp (database.rish.kyoto-u.ac.jp) をDNSに問いあわせています... 133.3.8.137
database.rish.kyoto-u.ac.jp (database.rish.kyoto-u.ac.jp)|133.3.8.137|:80 に接続しています... 接続しました。
HTTP による接続要求を送信しました、応答を待っています... 200 OK
長さ: 225280 (220K) [application/x-tar]
`./data/2021/02/15/Z__C_RJTD_20210215120000_RDR_JMAGPV__grib2.tar' に保存中

Z__C_RJTD_20210215120000_RDR_JMAGPV__grib2.tar             100%[========================================================================================================================================>] 220.00K   424KB/s    時間 0.5s

2021-07-11 22:02:30 (424 KB/s) - `./data/2021/02/15/Z__C_RJTD_20210215120000_RDR_JMAGPV__grib2.tar' へ保存完了 [225280/225280]

Z__C_RJTD_20210215120000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin
Z__C_RJTD_20210215120000_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_grib2.bin
*** Update check_pdt_size: pdt_size unknown, pdt=50008 ***
1:0:d=2021021512:var discipline=0 center=34 local_table=1 parmcat=1 parm=201:surface:-10-0 min acc fcst:

この例では、2021年2月15日12時0分のデータが、./data/2021/02/15/にダウンロードされ、

Z__C_RJTD_20210215120000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.dat

というバイナリファイルを生成しています。

全国合成レーダーGPVの描画

前項で作ったバイナリファイルを読み込みます。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import pandas.tseries.offsets as offsets
import os, sys

def mkdata_time(time):
    utc = time - offsets.Hour(9)
    day_dir = utc.strftime("%Y/%m/%d")
    infile_path  = "{}/{}/Z__C_RJTD_{}00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL.dat".format(Opath, day_dir, utc.strftime("%Y%m%d%H%M"))

    if not os.path.exists(infile_path):
        import download_RDR_JMAGPV
        download_RDR_JMAGPV.download_time(utc)

    ## 入力ファイルの読み込み
    rain = np.fromfile(infile_path, dtype="float32", sep="").reshape((mlon0, mlat0), order="F")
    print("read ({}) OK".format(infile_path))
    rain[rain>10000]=-10

if __name__=="__main__":
    ## データの日時(10分間隔)を指定
    time = pd.Timestamp(2021,2,15,21,0)
    download_time(time)

Z__C_RJTD_20210215120000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.dat

をinfile_pathに指定して読み込んでいます。

格納されている値は降水量(mm/hr)です。

また、先ほどはスルーしましたが、ファイル名のタイムスタンプ「20210215120000」は2021年2月15日12時は世界標準時(UTC)のため、時差が+9時間の日本の地方標準時(JST)にすると2021年2月15日21時になります。

続いて、cartopyを使って地図上に描画します。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
np.set_printoptions(threshold=np.inf)
import pandas as pd
import pandas.tseries.offsets as offsets
import os, sys, glob
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shapereader

import matplotlib.ticker as mticker
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter

import download_RDR_JMAGPV

def plot_2Dgrid_map(LON, LAT, data, slon, elon, slat, elat, time):

    fig = plt.figure(figsize=(8, 6))
    plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.9)
    
    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
    ax.set_title(time.strftime("%Y-%m-%d %H:%M"), fontsize=15)
    ax.set_extent([slon, elon, slat, elat])
    ax.coastlines(resolution='10m', lw=0.5)
    ax.add_feature(cfeature.LAND)

    xticks = np.linspace(slon,elon,7)
    yticks = np.linspace(slat,elat,7)

    # 経度線・緯度線を描く
    gl = ax.gridlines(crs=ccrs.PlateCarree())
    gl.xlocator = mticker.FixedLocator(xticks)
    gl.ylocator = mticker.FixedLocator(yticks)

    #目盛を描く緯度経度の値を設定
    ax.set_xticks(xticks,crs=ccrs.PlateCarree())
    ax.set_yticks(yticks,crs=ccrs.PlateCarree())
    
    #目盛の表示形式を度数表記にする
    latfmt=LatitudeFormatter()
    lonfmt=LongitudeFormatter(zero_direction_label=True)
    ax.xaxis.set_major_formatter(lonfmt)
    ax.yaxis.set_major_formatter(latfmt)
    
    #目盛のサイズを指定
    ax.axes.tick_params(labelsize=10)

    #国境の描画
    shpfilename = shapereader.natural_earth(resolution='10m',
                                            category='cultural',
                                            name='admin_1_states_provinces')
    reader = shapereader.Reader(shpfilename)
    provinces = []
    for province in reader.records():
        if(province.attributes["admin"] == "Japan"):
            provinces.append(province)

    data[data==0]=np.nan
    clevs = np.array([0,0.5,1,2,3,4,5,7,10,15,20,30,50,100,9999])
    cont_level = clevs[:-1]
    cont_norm = clevs

    newcmp = cm.get_cmap('cool', 512)
    newcmp.set_under("silver")
    cnorm = colors.BoundaryNorm(cont_norm,512)
    
    cs = ax.contourf(LON, LAT, data, cmap=newcmp, levels=cont_level, norm=cnorm, extend="both")
    cb = plt.colorbar(cs, orientation="vertical", ticks=clevs[:-1])    
    cb.ax.tick_params(labelsize=15)

    plt.savefig("{}/RDR_{}.png".format(Ipath, time.strftime("%Y%m%d-%H%M")))
    plt.show()
    plt.close() 
    return

def mkdata_time(time):
    utc = time - offsets.Hour(9)
    day_dir = utc.strftime("%Y/%m/%d")
    infile_path  = "{}/{}/Z__C_RJTD_{}00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL.dat".format(Opath, day_dir, utc.strftime("%Y%m%d%H%M"))
    # return
    if not os.path.exists(infile_path):
        download_RDR_JMAGPV.download_time(utc)

    ## 入力ファイルの読み込み
    rain = np.fromfile(infile_path, dtype="float32", sep="").reshape((mlon0, mlat0), order="F")
    print("read ({}) OK".format(infile_path))
    rain[rain>10000]=-10

    ## 描画を行う
    ## plot_2Dgrid_map(LON0, LAT0, rain, 経度西端, 経度東端, 緯度南端, 緯度北端)
    plot_2Dgrid_map(LON0, LAT0, rain, 130, 145, 30, 45, time)
    return

##-------------------------------------------------------------
if __name__=="__main__":

    ## 気象庁全国レーダー合成図の読み込み・図示を行うプログラム

    ## 入力ファイルのディレクトリ
    Opath = "./data"
    Ipath = "./img/single"
    os.makedirs(Ipath, exist_ok=True)

    ## 入力ファイルの次元の指定(水平解像度約1km)
    ## 経度方向: 118°〜150° (1/80°間隔)
    ## 緯度方向: 20°〜48° (1/120°間隔)
    mlon0, mlat0 = 2560, 3360
    slon0, elon0, rlon0 = 118, 150, 1/80
    slat0, elat0, rlat0 = 20, 48, 1/120

    lon0 = np.arange(slon0, elon0, rlon0)
    lat0 = np.arange(slat0, elat0, rlat0)
    LON0, LAT0 = np.meshgrid(lon0, lat0)
    LON0, LAT0 = LON0.T, LAT0.T

    ## 描画する時間の指定(年,月,日,時,分)
    ## データは10分ごと(前10分の雨量が記録されている)

    ## 単一時間のデータを描画
    time = pd.Timestamp(2021,7,2,10,0)
    mkdata_time(time)

    ## 連続データの描画
    start_time = pd.Timestamp(2021,5,1,12,0)
    end_time = pd.Timestamp(2021,5,1,13,0)

    time_list = pd.date_range(start=start_time, end=end_time, freq="10min")
    print(time_list)
    for time in time_list:
        mkdata_time(time)

このプログラムを実行すると、以下のような図ができます。

参考にさせていただいたサイト

Quitta 気象×Python 〜全国合成レーダー〜https://qiita.com/OSAKO/items/ef042f80ec63dd288225