気象レーダーによる降水量メッシュデータを、pythonでダウンロード、描画していきます。（Ubuntu-20.04 on WSL2）
気象庁全国合成レーダーGPVとは
気象庁は全国に20台の気象レーダを保有しており、10分間隔で雨の観測が行われています。1つ1つのレーダは半径約数100kmの範囲の観測に限られますが、これらの多数のレーダからの情報をひとまとめにしたメッシュ状のマップデータが提供されています。これを全国合成レーダGPVといいます。GPVは「Grid Point Value (格子点値)」の略で、メッシュ状のデータという意味です。
・エコー強度（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