講義のページ/Pythonのお勉強/RadarData

気象庁合成レーダーデータ †

概要
気象庁は合成レーダーデータを作成し、提供している。
データの読み取り
wgrib2 で netCDF 形式に変換してから読み込む。これを実行するには、Linux のような環境を用いることが一般的で、Windows で実現している例は少ない。本格的なデータを扱うためには、いずれは、Linux のような環境が必要になると考えられる。
しかし、当面は、Windows マシンに何かをインストールすることはせずに、実行することにする。そこで、Google Colaboratory で wgrib2 を実行する。

参考URL †

全般
- 増田先生によるプログラム
  - 解析雨量を Pythonで読み、作図する (試作版)
  - 全国合成レーダーデータを Pythonで読み、作図する (試作版)
    網羅的な情報が掲載されており、実行環境も含めてプログラムが示されている。まずは参照したいページ。地図は cartopy。
- pythonでGRIB2形式のバイナリファイルからレーダーエコー強度画像を作成する（画像作成編）
  地図に basemap を使った場合。
- 気象庁の1kmメッシュ解析雨量GPV：GRIB2形式をpython：hvplotでレンダリングする
  hvplot による解析雨量の結果で、レーダーデータそのものでない。情報も豊富。地図は geoview。
wgrib2
- Windows環境でGRIB2形式データをCSV形式に変換する（ある地点のデータを抽出する）
  wgrib2 は、単にデータ形式の変換だけでなく、データの抽出にも使える。
- Colaboratoryでwgrib2を使えるようにした
  wgrib2 を Google Colaboratory で使うとき、wgrib2 を Google ドライブに保存する方法。
- Colaboratoryを使ってブラウザ上だけで気象レーダー画像を作成する
  wgrib2 を Google Colaboratory で使う、というだけでなく、図にするまでの一連の流れ。
- ecCodesに入ったGRIB2 Template 5.200/7.200サポートを用いてpygribで気象庁全国合成レーダーGPVを読む
  grib2 という形式は、下位の仕様にバリエーションがあるために、grib2 用の python ライブラリ pygrib が使えるとは限らない。その辺の事情も含めて。
データの解説
- 気象×Python　〜全国合成レーダー〜
  エコー強度（降水強度と関連）、エコー頂高度の区別と、格納されているデータとの対応表。

Google Colaboratory でのプログラミング概要 †

Google Colaboratory で wgrib2 を動かすための準備

Google Colaboratory は、Debian 系のLinux として動く。そこで、実行可能なwgrib2 を、作成して保存しておくことが考えられる。ところが、Google Colaboratory では、接続が切れると保存したファイルは削除される。そこで、Googleドライブに保存し、それを実行できるようにする。まずは、通常のGoogleドライブとの混用が起こらないように、Google Colaboratory で使うツールを置く専用フォルダ Colab Env というフォルダを作る。 wgrib2 のファイルをその下で展開してコンパイル、設置する。

Google Colaboratory から、Google ドライブをマウントする際、許可を求められる。許可する。

この操作は、同じ環境（同じGoogle のID）を使っている間は、1度きりの実行で良い。

from google.colab import drive
drive.mount ('/content/drive') 
 
!mkdir -p /content/drive/MyDrive/ColabEnv/src
!mkdir -p /content/drive/MyDrive/ColabEnv/bin
!mkdir -p /content/drive/MyDrive/ColabEnv/DATA
%cd /content/drive/MyDrive/ColabEnv/src
!wget https://www.ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/wgrib2.tgz.v3.1.2
!tar xvzf wgrib2.tgz.v3.1.2
%cd grib2
%set_env CC gcc
%set_env FC gfortran
!make
!cp wgrib2/wgrib2 ../../bin/
!ls -l ../../bin/

気象庁レーダーデータのダウンロード

研究目的で、京大生存圏研究所からデータをありがたくダウンロードする。ダウンロードしたデータも、普通に使うと、毎回毎回、Google Colaboratory の接続が切れると削除されてしまう。そこで、通信負荷を減らすために、一度ダウンロードしたデータは、保存しておくようにする。また、上記の wgrib2 を用いて、ダウンロードしたデータを読み取りやすい netCDF 形式に変換する。これを関数を使って実行する。

次の操作は、Google Colaboratory の接続ごとに行う。一度行っておけば、Google Colaboratory との接続が切れるまでは、再び実行する必要はない。

-
!

from google.colab import drive
drive.mount ('/content/drive') 
 
import os
os.environ['PATH'] += ":/content/drive/MyDrive/ColabEnv/bin"
# print(os.environ['PATH'])
!chmod u+x /content/drive/MyDrive/ColabEnv/bin/wgrib2

def getradar(datehour, n10):
    datehour = str(datehour)
    n10 = str(n10)
    %cd /content/drive/MyDrive/ColabEnv/DATA
    !echo Z__C_RJTD_{datehour}{n10}000_RDR_JMAGPV__grib2.tar
    !if [ ! -e Z__C_RJTD_{datehour}{n10}000_RDR_JMAGPV__grib2.tar ]; then wget http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synthetic/original/{datehour[0:4]}/{datehour[4:6]}/{datehour[6:8]}/Z__C_RJTD_{datehour}{n10}000_RDR_JMAGPV__grib2.tar; tar xf Z__C_RJTD_{datehour}{n10}000_RDR_JMAGPV__grib2.tar; fi
    !for f in *.bin; do if [ ! -e "${f%.bin}.ncdf" ]; then wgrib2 "$f" -netcdf "${f%.bin}.ncdf"; fi; done

その後次のコマンドでファイルをダウンロードし、読み込みやすい netcdf 形式に変換できる。

getradar(2023062102,2)

この例では、2023年6月21日2時(UTC) の 10 分ごとのデータ( 0 ～ 5 )のうち、3番目を取得している。

取得したnetCDF 形式のファイルは、次のようにして確かめられる。

-
!
-
!

# エコー強度のデータ
!ls *1km*.ncdf
# エコー頂高度のデータ
!ls *5km*.ncdf

ファイ名とデータ名称の対応関係は、気象庁情報カタログ気象レーダー観測にある。

データ名称	ファイル名
2.5kmメッシュ全国合成レーダーエコー頂高度GPV	Z__C_RJTD_yyyyMMddhhmmss_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_grib2.bin
1kmメッシュ全国合成レーダーエコー強度GPV	Z__C_RJTD_yyyyMMddhhmmss_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin

レーダーデータの読み込み
netCDF 形式のファイルを読み込む。

準備
次のコードは、接続が切れるまで、一度だけ行えば良い。

-
|
!
 
-
!

# 必要なライブラリを導入する。
# ! apt install python3-netcdf4 : これでは何故か使えない。
! pip install netCDF4 cartopy
 
# データがある場所に移動する。
%cd /content/drive/MyDrive/ColabEnv/DATA

netCDF では変数名が重要である。ある特定の日時のデータを読み込み、netCDF のデータの変数名を抽出する。
これは今後、実行する必要がない。

-
|
!
-
!
 
-
!
 
-
|
!
-
!
 
-
!

## エコー頂高度
# netCDFファイルを開く
dataset = netCDF4.Dataset('Z__C_RJTD_20230621022000_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_grib2.ncdf')
# 変数名を表示する
variable_names = dataset.variables.keys()
print(variable_names)
# ファイルを閉じる
dataset.close()
 
## エコー強度
# netCDFファイルを開く
dataset = netCDF4.Dataset('Z__C_RJTD_20230706022000_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.ncdf')
# 変数名を表示する
variable_names = dataset.variables.keys()
print(variable_names)
# ファイルを閉じる
dataset.close()

データの名前：実行した結果、次のことがわかる。これは、ファイルからデータを引き出すために重要な情報になる。
エコー頂高度のファイルの中のデータ名
```
'latitude', 'longitude', 'time', 'var0_15_192_surface'
```
エコー強度のファイルの中のデータ名
```
'latitude', 'longitude', 'time', 'var0_1_201_surface'
```

図を描く

xarray と matplotlib を使った例
これを使うと、netCDF4 はいらない。

-
!
 
 
-
!
 
-
!
 
 
 
-
!
 
 
 
 
 
 
 
-
!

## エコー頂高度の図を描く
import xarray as xr
import matplotlib.pyplot as plt
 
# netCDFファイルを開く
dataset = xr.open_dataset('Z__C_RJTD_20230621022000_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_grib2.ncdf')
 
# データ変数、緯度と経度の値を取得する（例: 'var0_15_192_surface'はデータ変数の名前）
data = dataset['var0_15_192_surface']
lat  = dataset['latitude']
lon  = dataset['longitude']
 
# データをプロットする
plt.contourf(lon, lat, data[0,:,:], levels=20)  # 最初の時間スライスを描画（0は時間のインデックス）
plt.colorbar()           # カラーバーを追加
plt.xlabel('Longitude')  # 軸ラベル
plt.ylabel('Latitude')
plt.grid(True)           # グリッド
 
plt.show()               # プロットを表示
 
# ファイルを閉じる
dataset.close()

netCDF, matplotlib, cartopy を使った例

 
 
 
 
 
-
!
 
-
!
 
 
 
-
|
!
 
 
-
!
-
!
-
!
 
 
-
!

import netCDF4 as nc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
 
# netCDFファイルを開く
dataset = nc.Dataset('Z__C_RJTD_20230621022000_RDR_JMAGPV_Gll2p5km_Phhlv_ANAL_grib2.ncdf')
 
# データ変数を取得する（例: 'data_var'はデータ変数の名前）
data = dataset.variables['var0_15_192_surface']
lat  = dataset.variables['latitude'][:]
lon  = dataset.variables['longitude'][:]
 
# データをプロットする
# fig = plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.contourf(lon, lat, data[0, ...], levels=20, transform=ccrs.PlateCarree(), cmap='jet')
 
# 世界地図を描画する
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
# ax.add_feature(cfeature.BORDERS, linewidth=0.5)
 
# プロットを表示する
plt.title('NetCDF Data with World Map')
plt.show()
 
# ファイルを閉じる
dataset.close()

参照先

講義のページ/Pythonのお勉強/RadarData

メニュー †

Arduino †

電子工作 †

研究 †

気象庁合成レーダーデータ †

参考URL †

Google Colaboratory でのプログラミング概要 †