気象庁の1kmメッシュ解析雨量GPV：GRIB2形式をpython：hvplotでレンダリングする

---

---

解析雨量について

気象庁が配信している解析雨量は、レーダー雨量と地上雨量を使って作成した面的な雨量グリッドデータです。

参考：気象庁｜解析雨量

今回は、気象業務支援センターのHP：技術資料（実際に配信されたデータによるサンプル提供）で公開されているサンプルデータを使用します。 直リンクはこちら

ファイル名はSRF_GPV_Ggis1km_Prr60lv_ANAL*となっており、 SRF（Short-Range Forecasting of Precipitation）の1km間隔のGPV（Grid Point Value）で前60分間雨量のanalysisデータを示しています。データの間隔は30分間です。

データ型の変換

元のデータはgrib2形式というバイナリで拡張はbinです。 今回は後の処理で扱いやすいようにNetCDF型に変換します。

grib2形式をデコードする方法は以下のとおりです。

• バイナリを直接読む

  • 気象業務支援センターで購入したDVDに付属するCのプログラム
  • pythonだと 気象庁解析雨量(GRIB2形式)をPythonで扱う - Qiita等
  • 他にも色々あります。

• ソフトウェアを使って型を変換する：有料、無料と大量にありますのでググって下さい。

今回はpythonを使いますが、grib2用のメジャーなモジュールがcfgribやpygribは、型が対応しておらず、使用することができませんでした。

そのため、今回はNOAA（National Oceanic and Atmospheric Administration）が公開しているgrib2用のユーティリティ：wgrib2を使います。

wgrib2のwindowsでのインストール方法（2022/7/10時点）

公式サイトClimate Prediction Center - wgrib2: grib2 utilityからダウンロードしてインストールします。

ページ最下部のSource Code and Compling Hints のSource codeより、「Windows10」-「v3.0.2」とリンクをたどります。

その中にある、dllファイルとwgrib2.exeをダウンロードし、例えば、C直下にwgrid2というとフォルダを作成して保存します。

これでwgrib2の設定はO.K.です。

wgrib2による解析雨量バイナリデータのNetCDF型への変換

サンプルデータをwgrib2でNetCDF型へ変換するスクリプトは次のとおりです。

@echo off
set PATH=C:\wgrid2;%PATH%

wgrib2 input.bin -netcdf output.nc

日ごとにフォルダで分類されたデータを一括で処理する場合は次のとおりです。

@echo off
set PATH=C:\wgrid2;%PATH%

for /D %%p in (*) do (
    echo %%p
    for %%f in ("%%p/*.bin") do (
    wgrib2 %%p/%%~nxf -netcdf %%~nf.nc
    )
)

参考：Windows環境でGRIB2形式データをCSV形式に変換する（ある地点のデータを抽出する）｜お天気データサイエンス｜日本気象株式会社

hvplotによるレンダリング

作成したNetCDFファイルを図化します。

hvplotの解説は以下の記事に詳しく書いています。

computational-sediment-hyd.hatenablog.jp

クリップと結合

解析雨量データの日本全域のため、使用する範囲のデータをクリップします。

今回は兵庫県三木市を対象としました。行政界データは、国土数値情報 | 行政区域データ からダウンロードしました。

また、NetCDFファイルは時間ごとの個別ファイルのため、1ファイルに結合します。

使用するライブラリ

以下のとおりです。

import xarray as xr
import rioxarray
import pandas as pd
import geopandas as gpd
import glob
import numpy as np
from shapely.geometry import mapping
 
rioxarray.__version__
# '0.0.21'
gpd.__version__
# '0.6.1'

少し珍しいものとしてはrioxarrayがあります。Raster I/Oという名前のとおり、xarrayをラスタ（ジオメトリ情報付きも可）のように取り扱うことができるモジュールです。

クリップ

抽出する範囲のポリゴンデータをgeojson形式として、xarray.rio.clipの引数に与えること、矩形領域でクリップしたxarray.datasetが作成されます。なお、矩形領域内のポリゴン範囲外はnanとなります。

参考：python - Extracting data within geometry (shape) - Geographic Information Systems Stack Exchange

gdf = gpd.read_file("N03-19_28_190101.geojson")
poly = gdf[gdf['N03_004']=='三木市'].geometry.values[0]

f = 'input.nc'
ds = xr.open_dataset(f)
d = ds.rio.write_crs('epsg:6668', inplace=True)
dsp = ds.rio.clip([mapping(poly)], gdf.crs)
dsp

雨量の単位はmm/hrです。

値はレベルに応じたデータ代表値になっています。対応は以下を参照下さい。

• 気象庁解析雨量(GRIB2形式)をPythonで扱う - Qiita
• http://sumisumi.cocolog-nifty.com/sumisumi/files/Rader-AME.pdf#page=5

データはfloat32型になっていますが、小数点以下はレベル1の0.4のみであるため、容量を節約したい場合は10倍にしてint型にできます。

結合

すべてのNetCDFファイルをクリップして結合します。

fs = glob.glob('*.nc')

dss = []
for f in fs:
    ds = xr.open_dataset(f)
    d = ds.rio.write_crs('epsg:6668', inplace=True)
    dsp = ds.rio.clip([mapping(poly)], gdf.crs)
    dss.append(dsp)

dsAll = xr.concat(dss, dim='time')

レンダリング

使用するライブラリ

hvplot、holoviews、geoviewsをインポートします。

import hvplot.xarray
import hvplot.pandas
import geoviews as gv
import holoviews as hv
import cartopy.crs as ccrs

gv.extension('bokeh', logo=False)
hv.extension('bokeh', logo=False)

座標系を設定

hvplotの混沌としているところでcrsの設定が独自です。バージョンに記述が変わっているので公式を確認したほうが良いです。

s = dsAll.attrs
s['crs'] = '+init=epsg:' + str(6668)
dsAll.attrs = s

タイムゾーンの変更：UTCからJST

気象系のデータはUTCで作成されているため、JSTに変更します。

computational-sediment-hyd.hatenablog.jp

time = dsAll['time'].values

s = pd.Series(time)
s = s.dt.tz_localize('UTC') 
s = s.dt.tz_convert('Asia/Tokyo')
s = s.dt.tz_localize(None)

dsAll['time'] = s.values

なお、ここではdsAll['time']のCoordinatesやAttributesは不要のため付けていないです。

レンダリング

作成したデータdsAllをhvplotを使ってレンダリングします。 なお、背景図には地理院タイルマップを使用しています。

gall = dsAll['var0_1_200_surface'].hvplot(groupby='time', clim=(0,25)
                                           , width=500, height=500
                                           , alpha=0.6, cmap='jet'
                                           , geo=True, rasterize=True, project=True, dynamic=False)

# 背景図
geomap = gv.WMTS('https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{Z}/{X}/{Y}.jpg')
# 行政界ポリゴン
gpoly = gv.Polygons(poly, crs=ccrs.PlateCarree()).options(fill_color=None, line_color='red')

go = geomap*gall*gpoly
go

フルサイズ版はこちら

ソースコード

• GitHub:RenderingGridPrecipitationData/rendering.ipynb at main · computational-sediment-hyd/RenderingGridPrecipitationData · GitHub

• nbviewer:Jupyter Notebook Viewer

GitHub

github.com

その他

• gif作成はScreenToGifを使用

computational-sediment-hyd.hatenablog.jp

• グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

• はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp