![[PukiWiki]](image/ROBO192.png "[PukiWiki]")
#author("2023-11-01T15:53:51+09:00","external:moriat","moriat") #author("2023-11-01T16:27:46+09:00","external:moriat","moriat") #topicpath ** 気象庁合成レーダーデータ [#e7c3e7c1] - 概要~ 気象庁は合成レーダーデータを作成し、提供している。 -- [[気象庁情報カタログ 気象レーダー観測>https://www.data.jma.go.jp/add/suishin/cgi-bin/catalogue/make_product_page.cgi?id=Radar]] -- [[気象業務支援センター>http://www.jmbsc.or.jp/jp/online/file/f-online30100.html]] -- [[京都大学生存圏研究所>http://database.rish.kyoto-u.ac.jp/arch/jmadata/synthetic-original.html]] - データの読み取り~ wgrib2 で netCDF 形式に変換してから読み込む。これを実行するには、Linux のような環境を用いることが一般的で、Windows で実現している例は少ない。本格的なデータを扱うためには、いずれは、Linux のような環境が必要になると考えられる。~ しかし、当面は、Windows マシンに何かをインストールすることはせずに、実行することにする。そこで、Google Colaboratory で wgrib2 を実行する。 *** 参考URL [#i59c78b2] - 全般 -- 増田先生によるプログラム --- [[解析雨量を Pythonで読み、作図する (試作版)>http://macroscope.world.coocan.jp/ja/edu/compex/kaiseki_uryo_python.html]] --- [[全国合成レーダーデータを Pythonで読み、作図する (試作版)>http://macroscope.world.coocan.jp/ja/edu/compex/zenkoku_radar_python.html]]~ 網羅的な情報が掲載されており、実行環境も含めてプログラムが示されている。まずは参照したいページ。地図は cartopy。 -- [[pythonでGRIB2形式のバイナリファイルからレーダーエコー強度画像を作成する（画像作成編）>https://qiita.com/kinosi/items/6f4f754f3f88f1fd74b4]]~ 地図に basemap を使った場合。 -- [[気象庁の1kmメッシュ解析雨量GPV：GRIB2形式をpython：hvplotでレンダリングする>https://computational-sediment-hyd.hatenablog.jp/entry/2022/07/10/230000#:~:text=wgrib2%E3%81%AE%20windows%20%E3%81%A7%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB%E6%96%B9%E6%B3%95%EF%BC%882022%2F7%2F10%E6%99%82%E7%82%B9%EF%BC%89%20%E5%85%AC%E5%BC%8F%E3%82%B5%E3%82%A4%E3%83%88%20Climate%20Prediction%20Center%20-,Code%20and%20Compling%20Hints%20%E3%81%AE%20Source%20code%20%E3%82%88%E3%82%8A%E3%80%81%E3%80%8CWindows10%E3%80%8D-%E3%80%8Cv3.0.2%E3%80%8D%E3%81%A8%E3%83%AA%E3%83%B3%E3%82%AF%E3%82%92%E3%81%9F%E3%81%A9%E3%82%8A%E3%81%BE%E3%81%99%E3%80%82]]~ hvplot による解析雨量の結果で、レーダーデータそのものでない。情報も豊富。地図は geoview。 - wgrib2 -- [[Windows環境でGRIB2形式データをCSV形式に変換する（ある地点のデータを抽出する）>https://ods.n-kishou.co.jp/tech/blog/detail/2869]]~ wgrib2 は、単にデータ形式の変換だけでなく、データの抽出にも使える。 -- [[Colaboratoryでwgrib2を使えるようにした>https://qiita.com/kogetaenoki/items/219c1d34cf1846452f79]]~ wgrib2 を Google Colaboratory で使うとき、wgrib2 を Google ドライブに保存する方法。 -- [[Colaboratoryを使ってブラウザ上だけで気象レーダー画像を作成する>https://qiita.com/kogetaenoki/items/2b50f627db3769280dce]]~ wgrib2 を Google Colaboratory で使う、というだけでなく、図にするまでの一連の流れ。 -- [[ecCodesに入ったGRIB2 Template 5.200/7.200サポートを用いてpygribで気象庁全国合成レーダーGPVを読む>https://zenn.dev/noritada/articles/grib2-template5_200-support-on-eccodes-pygrib]]~ grib2 という形式は、下位の仕様にバリエーションがあるために、grib2 用の python ライブラリ pygrib が使えるとは限らない。その辺の事情も含めて。 - データの解説 -- [[気象×Python 〜全国合成レーダー〜>https://qiita.com/OSAKO/items/ef042f80ec63dd288225]]~ エコー強度（降水強度と関連）、エコー頂高度の区別と、格納されているデータとの対応表。 ** Google Colaboratory でのプログラミング概要 [#p3d1cd0f] - Google Colaboratory で wgrib2 を動かすための準備~ ~ Google Colaboratory は、Debian 系のLinux として動く。そこで、実行可能なwgrib2 を、作成して保存しておくことが考えられる。ところが、Google Colaboratory では、接続が切れると保存したファイルは削除される。そこで、Googleドライブに保存し、それを実行できるようにする。 まずは、通常のGoogleドライブとの混用が起こらないように、Google Colaboratory で使うツールを置く専用フォルダ Colab Env というフォルダを作る。 wgrib2 のファイルをその下で展開してコンパイル、設置する。~ ~ Google Colaboratory から、Google ドライブをマウントする際、許可を求められる。許可する。~ ~ この操作は、同じ環境（同じGoogle のID）を使っている間は、1度きりの実行で良い。~ ~ #code(python){{ 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 との接続が切れるまでは、再び実行する必要はない。~ ~ #code(python){{ 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 }} #code(python){{ 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 形式に変換できる。 #code(python){{ getradar(2023062102,2) }} この例では、2023年6月21日2時(UTC) の 10 分ごとのデータ( 0 〜 5 )のうち、3番目を取得している。~ ~ 取得したnetCDF 形式のファイルは、次のようにして確かめられる。 #code(python){{ # エコー強度のデータ !ls *1km*.ncdf # エコー頂高度のデータ !ls *5km*.ncdf }} ファイ名とデータ名称の対応関係は、[[気象庁情報カタログ 気象レーダー観測>https://www.data.jma.go.jp/add/suishin/cgi-bin/catalogue/make_product_page.cgi?id=Radar]]にある。 |データ名称|ファイル名|h |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| |250mメッシュ全国合成レーダー降水強度 | Z__C_RJTD_yyyyMMddhhmmss_RDR_GPV_Ggis0p25km_Pri60lv_Aper5min_ANAL_grib2.bin.gz| ~ - レーダーデータの読み込み~ netCDF 形式のファイルを読み込む。~ -- 準備~ 次のコードは、接続が切れるまで、一度だけ行えば良い。 #code(python){{ # 必要なライブラリを導入する。 # ! apt install python3-netcdf4 : これでは何故か使えない。 ! pip install netCDF4 cartopy # データがある場所に移動する。 %cd /content/drive/MyDrive/ColabEnv/DATA }} netCDF では変数名が重要である。ある特定の日時のデータを読み込み、netCDF のデータの変数名を抽出する。~ これは今後、実行する必要がない。 #code(python){{ ## エコー頂高度 # 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 はいらない。 #code(python){{ ## エコー頂高度の図を描く 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 を使った例 #code(python){{ 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() }} [[参照先>https://colab.research.google.com/drive/1jjgtFY4aT5wbDosBXCDWu9TRQSW0VC3u?usp=sharing]]
