こんにちは、梶原です。
これは TECHSCORE Advent Calendar 2015 の13日目の記事です。
衝撃のビジュアル表現
あれは忘れもしません。大阪の環状線を止めるほどの台風が来た日の事です。
テレビを観ていると、NHKが天気・災害情報を放送していましたが、とてもカッコよく、わかりやすいビジュアル表現をしていました。画面上には「DIGITAL EARTH」と記載されており、NHKの新システムだったようです。
(参考) 「NHKデジタルアース」が凄い
特に風の表現に衝撃を受けたので、色々調べてみると、世界風速は、とても良くできていて、フリーバージョンはソースが公開されています。
そこで世界風速のソースや、気象データ処理や地図表現についても、調べたのでまとめてみました。
ソースはこちらに置いておきます。
地図表現
ビジュアル表現に、地図表現はつきもの！
という訳で、こちらは色々な手法があると思いますが、ブラウザで表示したかったので、D3.jsを利用する方法が一番簡単に実現できそうでした。海岸線を非常に簡単にSVG描画できます。海岸線情報自体は、Natural Earthを利用させていただく事にしました。
上記サイトで手に入るShapeデータを、GeoJSONに変換して、さらにより効率的に処理する為TopoJSONに変換して描画するのが一般的のようです。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
# 50m cultural の Admin0 - Countriesを利用しました。 $ wget http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries.zip $ unzip ne_50m_admin_0_countries.zip # gdalとtopojsonはインストール済みの想定です。 # GeoJSON化。日本周辺のみ抽出。 $ ogr2ogr -f GeoJSON -where "adm0_a3 IN ('JPN', 'KOR', 'PRK', 'TWN', 'CHN', 'RUS')" map_geo.json ne_50m_admin_0_countries.shp # TopoJSON化。オプションの使い方が今ひとつ分かってません。。 $ topojson --id-property su_a3 -p NAME=name -p name -o map_topo.json map_geo.json |
出来上がったmap_topo.jsonを利用するhtml/javascriptを作成し、httpサーバに置くと
ね、簡単でしょ(汗)！ソースはgithubにおいてみました。
こだわりポイントとしては、
- 緯度経度線があるとかっこいい
- 海岸線はくっきり
- 地図図法の"それっぽさ"は大事なので、正射図法を利用
あたりでしょうか。なおモンゴルが欠けているのは秘密です。余力のある方は、Natural Earthの湖データを利用して琵琶湖を描いてみてください。
地図を描画できる事が一番重要ですが、座標変換ロジックが関数として再利用できる事も重要です。
なぜなら次の気象データは、「緯度A、経度Bの風向きは〜」という形で得られる為、緯度A, 経度B地点を、画面上の座標に変換できる必要があります。上記の場合、下記のような変数に入れて後で利用します。
|
1 2 3 4 5 6 7 |
... var projection = d3.geo.orthographic() .scale(width * 2.0) .translate([width / 2, height / 2]) .clipAngle(90) .rotate([225, -38]); ... |
気象データ処理
最近気象データの多くは、「GRIB2」形式というファイルフォーマットで配信される事が多いです。しかし、このフォーマットはとても複雑で、とても自分でparserを書けるようなものではない事がわかりました。
行いたい事はデータを利用する事ですので、先人の肩にのるべく調べると、 wgrib2というツールで、ダンプできる事がわかりました。下記のような手順でコンパイル可能です。
|
1 2 3 4 5 |
# gfortranがインストールされており、コンパイル環境が整っている前提です。 $ wget ftp://ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/wgrib2.tgz $ tar xzvf wgrib2.tgz $ cd grib2 $ make |
(参考) コンパイル手順
wgrib2/wgrib2が出来上がったツールになります。試しに動かしてみると...
|
1 2 3 4 5 6 7 8 9 10 11 12 |
$ wgrib2/wgrib2 wgrib2 v0.2.0.3 11/2015 Wesley Ebisuzaki, Reinoud Bokhorst, John Howard, Jaakko Hyvätti, Dusan Jovic, Kristian Nilssen, Karl Pfeiffer, Pablo Romero, Manfred Schwarb, Arlindo da Silva, Niklas Sondell, Sergey Varlamov stock build -0xSec inv X Hex dump of section X (0..8) -aerosol_size inv optical properties of an aerosol -aerosol_wavelength inv optical properties of an aerosol -bitmap inv bitmap mode ... -set_ext_name init X X=0/1 extended name on/off -set_regex init X set regex mode X = 0:extended regex (default) 1:pattern 2:extended regex & quote metacharacters -tigge init use modified-TIGGE grib table -transient init X make file X transient, CW2 |
気象データ自体はNOMADSや、学術目的なら地球流体電脳倶楽部あたりで見つける事ができました。
# 業務で気象データを扱う際は、気象業務支援センターより購入する必要があります。
ここではNOMADSから、下記ファイルー条件でデータダウンロードしてみます。grib filterから条件を選択してダウンロードできます。
- 2015/12/11
- 緯度経度の刻みは0.25
- 高度は、10 m above ground
- 風の強さなので、VGRD/UGRD
http://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f000&lev_10_m_above_ground=on&var_UGRD=on&var_VGRD=on&dir=%2Fgfs.2015121100
データを先ほどのwgrib2で見てみると、どんなデータが含まれているかわかります。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
$ ./wgrib2 -V data/2015121100.dat 1:0:vt=2015121100:10 m above ground:anl:UGRD U-Component of Wind [m/s]: ndata=1038240:undef=0:mean=-0.108394:min=-19.33:max=29.75 grid_template=0:winds(N/S): lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48 lat 90.000000 to -90.000000 by 0.250000 lon 0.000000 to 359.750000 by 0.250000 #points=1038240 2:829559:vt=2015121100:10 m above ground:anl:VGRD V-Component of Wind [m/s]: ndata=1038240:undef=0:mean=-0.0255382:min=-22.68:max=27.58 grid_template=0:winds(N/S): lat-lon grid:(1440 x 721) units 1e-06 input WE:NS output WE:SN res 48 lat 90.000000 to -90.000000 by 0.250000 lon 0.000000 to 359.750000 by 0.250000 #points=1038240 |
wgrib2には非常に多くのオプションがありますので、
- データのみ
- 南北方向を北から南に変換
- 1データ4byte浮動小数点(big endian)
としてダンプしてみます。
|
1 2 3 4 5 |
$ ./wgrib2 -d 1 data/2015121100.dat -order we:ns -no_header -ieee data/2015121100-U.ieee $ ./wgrib2 -d 2 data/2015121100.dat -order we:ns -no_header -ieee data/2015121100-V.ieee $ ls -la data/*.ieee -rw-r--r-- 1 kennyj staff 4152960 12 12 15:00 data/2015121100-U.ieee -rw-r--r-- 1 kennyj staff 4152960 12 12 15:00 data/2015121100-V.ieee |
いい感じですね。1440 * 721 * 4 = 4152960なので、ちゃんとダンプできている事がわかります。
ただ、結構ファイルサイズが大きくなるので、精度を落として1データ1byteにするとか、必要な緯度経度内に絞る、gzip圧縮する等、色々工夫はできると思います（ここでは説明を簡単にする為、このまま進めます)。
また〜U.ieeeのあとに〜V.ieeeを連結して１度にダウンロードできるようにします（こちらも簡単の為）。
|
1 |
cat data/2015121100-U.ieee data/2015121100-V.ieee > dest/data/wind.ieee |
風のビジュアル表現に挑戦
さて、ここまでで道具は揃ったので、風のビジュアル表現をしてみようと思います。
先ほど処理した風データは、ある地点でのUV方向の風の強さを表現したデータ、すなわちベクトル場といえる為、そこを漂う粒子として可視化します。
具体的には
- 風データをダウンロードして、javascript上でベクトル場データを準備。
- 多くの粒子を、範囲内にランダムに配置、寿命を設定する。
- 毎描画フレーム毎に
- 粒子座標の風の強さを、ベクトル場データから補間によって求め、次の粒子の座標を再計算する。
- 範囲外に出てしまった、もしくは寿命に達したら、再度ランダムに粒子を配置する。
- 粒子の移動（前回と今回の座標）を線で描画。
というステップで描画してみます。いわゆるパーティクルシステムのアプローチの一つになります。
実際に描画時には幾つか必要なテクニックがありました。
- ベクトル場データの補間は簡単の為、近傍4点での線形補間にしました。
- canvasでの描画の際に、beginPath / moveTo / lineTo / stroke を、1粒子毎に実行すると一気にフレームレートが落ちます。風の強さで色を変える等の機能を実装する際には注意が必要になります。
- 粒子の軌跡を残像として描画すると見た目がぐっと良くなります。こちらは、冒頭の世界風速で利用されているテクニックが大変参考になりました。
|
1 2 3 4 5 6 7 8 |
var ctx = canvas.node().getContext("2d"); // canvasを少し暗く ctx.fillStyle = "rgba(0, 0, 0, 0.90)"; var prev = ctx.globalCompositeOperation; ctx.globalCompositeOperation = "destination-in"; ctx.fillRect(0, 0, canvas.attr("width"), canvas.attr("height")); ctx.globalCompositeOperation = prev; |
詳細は、ソースをご確認ください。当記事作成の為急いで書いたので、不備はご容赦ください(汗)
動かすと下記のような風の表現がアニメーションします。 是非動かしてみてください。
githubからclone後、01~03のshを実行すると、必要なデータは、ダウンロードできると思います（関連するプロダクトは参照ページを参考にご準備下さい)。
世界風速や、DIGITAL EARTHで表現している手法を学ぶ入り口になっている事は、見てとれると思います。
さらに...
今回は、風データを固定された日時のデータを１つ利用して描画しました。
お気づきの方もいらっしゃると思いますが、時系列に並んだ複数のデータを利用して、補間しつつ描画する事によって、「時事刻々と変わる風向き、風の強さ」を表現する事ができます。
今回はコードの整理が間に合わなかったので、またの機会がありましたら投稿してみようと思います。
最後まで読んで頂きありがとうございました。
