2015年12月2日水曜日

Earthの風データをFOSS4Gで可視化

FOSS4G Advent Calendar 2015の12/2の記事です。
2013年、2014年は参加してなかったので2012年以来となります。


ご存じの方も多いと思うのですが、気象データの美しい可視化サイトで有名なEarthというサイトがあります。


こちらはCameron Beccarioさんという方がNOAAやNASAなどの研究機関が公開している気象データからD3などのツールを使って可視化を行っており、大気や海洋の様子をわかりやすいアニメーションで見ることができます。

もちろんこのまま表示をいじっていろいろ見るだけでも十分楽しいのですが、研究者的には元データってどんなデータなんだろう?とか、どこからダウンロードして手元でいじれるのかなどと考えてしまうわけです。

そこでこの記事では、EARTHで使われている気象データのうち、風データのみに着目してデータ配布サイトから入手、フォーマット変換をしてからQGISで可視化するまでを紹介したいと思います。


データの入手

一応EARTHの説明ページを見るとデータの種類ごとに、元データの配布機関へのリンクは張っています。風データの場合はNational Centers for Environmental Prediction (NCEP) により提供されている気象予測モデルGlobal Forecast System (GFS)の気象データを使っているようです。


ぱっと見ただけではそれらしいデータへのリンクがたくさんあってよくわかりませんorz.
そこで、githubに公開されているEARTHのソースコードを確認してみます。


ページの下の方に行くと "getting weather data" というそれらしい説明部分がありました。NOAA Operational Model Archive and Distribution System(NOMADS)というデータ配布システムからダウンロードできるみたいです。
※そういえばRにも rNOMADS というこのサイトのデータへアクセスできるライブラリがありますね。

でNOMADSのページから "GFS 0.25 Degree"の "grib filter"というリンクをたどり。
ダウンロードしたい日時(GMT)を選びます。ここでは2015/12/01 00時のデータを使いました。
(↑だんだんWebページのスクリーンショットを貼るのが面倒になってきましたので文字ベースで)


GRIB Filterのページでは以下の指定を行います。

  • ドロップダウンリストで "gfs.t00z.pgrb2.0p25.f000" (fのあとはその日時から何時間後の予測データなのかを示す。f003だと3時間後)
  • levelとして "10 m above ground" にチェック
  • variable(変数)として "UGRD" と "VGRD" を選択してダウンロードを開始

※とりあえず地上の風(1000 hPa)を見たいので10 mの高さで、風の東西成分(UGRD)、南北成分(VGRD)を指定しています。

ダウンロードされたデータは、 "gfs.t00z.pgrb2.0p25.f000" という名前のGRIB2形式のデータです。GRIB2形式はGDALでも対応しているのでQGISで開くことができます。

Band1のU(東西成分)が赤色、Band2のV(南北成分)が緑色で表示されているみたいですが、このままではよくわかんないですね…。


データの変換

この2バンドのラスターGRIB2を以下の2種類のデータに変換します。

  1. 風向・風速を示すベクターデータ
  2. 風速を示すラスターデータ


風向・風速を示すベクターデータへの変換

WindowsにQGIS 2.12.1(スタンドアロン版)がインストールされているとして話を進めていきます。

先ほどのGRIB2データがDドライブの直下に置いてあるとして、MSYS (QGISと一緒にインストールされるUNIXライクなコンソール)を起動して以下のコマンド入力していきます。#で始まる行はコマンドじゃなく説明です。

  cd d

#gdalinfo コマンドでGRIB2データの情報を確認

  gdalinfo gfs.t00z.pgrb2.0p25.f000

#GRIB2形式のband1(UGRD)とband2(VGRD)をそれぞれCSVに変換

  gdal2xyz.py -band 1 -csv gfs.t00z.pgrb2.0p25.f000 wind_x.csv
  gdal2xyz.py -band 2 -csv gfs.t00z.pgrb2.0p25.f000 wind_y.csv

#改行コードをUNIX形式に変換

  dos2unix wind_x.csv
  dos2unix wind_y.csv

#1つのCSVファイルに結合

  echo "lon,lat,u,v" > wind.csv
  paste -d"," wind_x.csv wind_y.csv | awk -F, '{print $1","$2","$3","$6}' >> wind.csv

これで以下のように、longitude、latitude、U成分、V成分のCSVファイルが出来上がります。
lon,lat,u,v
0.000,90.000,-2.79,-1.51
0.250,90.000,-2.8,-1.5
0.500,90.000,-2.8,-1.48

あとはQGISでCSVファイルを読み込み、Shapefile形式で別名で保存を行って下さい。
CSVファイル読み込みウィンドウ
X fieldにlon、Y fieldにlatを指定

風速を示すラスターデータへの変換

#再びMSYS上で以下のようにコマンドを入力して、GRIB2形式のファイルを風速のラスターデータ(GeoTIFF)に変換します。
  gdal_calc.py -A gfs.t00z.pgrb2.0p25.f000 --A_band=1 -B gfs.t00z.pgrb2.0p25.f000 --B_band=2 --calc="sqrt(A*A+B*B)" --outfile=scalar.tif


変換したデータを読み込んで色を設定すると以下の様になります。
風速が強いところが赤、弱いところが青



QGISでの操作

風のベクトルプロットをするにあたって追加でプラグインが必要になります。
QGISのプラグインメニューより "Vector field renderer" プラグインをインストールして下さい。

その後、風速のラスターデータの上に、風向・風速のベクターデータを読み込んで重ねて下さい。ポイントの数が多いので、読み込んだ直後は以下の図のように真っ黒に見えます。日本周辺などに拡大してみましょう。

データ全体表示
よくわからない…。

日本の北海道あたりに拡大するとポイントデータが格子状に並んでいるのがわかります。


このそれぞれのポイントデータが風の東西成分と南北成分を持っていいるのでその値を元にそれぞれのポイントにおけるベクトルを表示していきます。

ベクトルデータのプロパティをみると、スタイルタブで "ベクターフィールド" という項目が選べるようになっています。
そこで以下のように "X attribute" に "u"、 "Y attribute" に"v" を指定して、"Base size" を "0" に指定してOKをクリックします。

”Vector field” の設定


すると最終的に以下のように風向・風速の分布を表示することができます。


まあ実際解析する際はこのようなスナップショットではなく、一定期間ごとに統計処理を行って解析するのですが、そこまでやる場合はプログラミングでの処理が必要になってきますね。そこら辺は追々Rを使った処理についてまたブログに書きたいと思います。

ちょっと学生への手順書としてもかねての記事なので説明が冗長な部分もありますが…。
こんなところで終わります。