今更ですが、気象データを解析してみた

いろいろありまして、気象データを使う必要があるかもしれないので、ちょっとやってみました。

気象データ

気象業務支援センターから取得する気象データは国際気象機関が定めるフォーマットに従っています（GRIB2形式）。
もらえるデータはバイナリーデータなので、普通に開いてもわかりません。
なので、フォーマットに従ってデータを読める形に解析する必要があります。

wgrib2

解析部分を自前で作ってもいいのですが、アメリカ海洋大気庁（NOAA）さんが解析ツールを公開してるみたいなので使用させてもらいます。
　→ http://www.ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/

MacOS 10.12 (Sierra) にインストールしてみる

今回は最新版の wgrib2.tag.v2.0.7 をダウンロードしてインストールしてみました。

$ curl -o wgrib2_2.0.7.tgz http://www.ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/wgrib2.tgz.v2.0.7
$ tar zxvf wgrib2_2.0.7.tgz

$ cd grib2　← ディレクトリ名注意！ wgrib2 という名前じゃないので一瞬 あれっ？ ってなりました。。
$ make

ですが、なんかエラーが！
調べてみると FORTRAN がいるっぽい。。。＃FORTRAN...
というわけでインストール。

$ brew install gcc

＃brew で gcc をインストールするともれなく fortran がついてくるらしい、、、
＃そしてあとで brew でインストールした gcc を使うので、インストールします。

インストールが完了したら、確認のために README.Mac を見てみる。
＃よいこははじめに確認しましょう :)

＜略＞
… setting the following environment variables before doing the "gmake -f makefile"

環境変数を指定して、ということなので指定する。
CC と CXX と F77 を先ほど brew でインストールした gcc, g++, gfortran にする。
＃xcode の gcc とかだとエラーになるのですよ。

$ export CC=/usr/local/bin/gcc-7
$ export CXX=/usr/local/bin/c++-7
$ export F77=/usr/local/bin/gfortran-7

あとはおまけ :)
$ export CFLAGS="-O2 -m64"
$ export CXXFLAGS="-O2 -m64"
$ export FFLAGS="-O2 -m64"
$ export wCPPFLAGS="-O2 -m64"

で再度 make 。＃あ、gmake -f makefile はしません :p
エラー。。。
f77 というコマンドがないとかって怒られたので、リンクを作る。

# ln -s /usr/local/bin/gfortran-7 /usr/local/bin/f77

で再度 make 。

できましたー！
確認しましょう。

$ ./wgrib2/wgrib2 --version
v0.2.0.7 12/2017  Wesley Ebisuzaki, Reinoud Bokhorst, John Howard, Jaakko Hyvätti, Dusan Jovic, Daniel Lee, Kristian Nilssen, Karl Pfeiffer, Pablo Romero, Manfred Schwarb, Gregor Schee, Arlindo da Silva, Niklas Sondell, Sam Trahan, Sergey Varlamov

パスの通ってるところにリンクしておきましょう。

# ln -s ./wgrib2/wgrib2 /usr/local/bin/wgrib2

サンプルデータの解析

ここから「メソ数値予報モデルGPV（MSM）」を取得します。今回は地上面の情報にしてみます。
　→ http://www.jmbsc.or.jp/jp/online/c-onlineGsample.html

では解析

$ wgrib2 
Z__C_RJTD_20171210120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin

すると以下の結果が出てきます。

1.1:0:d=2017121012:PRMSL:mean sea level:anl:
1.2:0:d=2017121012:PRES:surface:anl:
1.3:0:d=2017121012:UGRD:10 m above ground:anl:
1.4:0:d=2017121012:VGRD:10 m above ground:anl:
1.5:0:d=2017121012:TMP:1.5 m above ground:anl:
1.6:0:d=2017121012:RH:1.5 m above ground:anl:
1.7:0:d=2017121012:LCDC:surface:anl:
1.8:0:d=2017121012:MCDC:surface:anl:
1.9:0:d=2017121012:HCDC:surface:anl:
1.10:0:d=2017121012:TCDC:surface:anl:
1.11:0:d=2017121012:PRMSL:mean sea level:1 hour fcst:
1.12:0:d=2017121012:PRES:surface:1 hour fcst:
1.13:0:d=2017121012:UGRD:10 m above ground:1 hour fcst:
1.14:0:d=2017121012:VGRD:10 m above ground:1 hour fcst:
1.15:0:d=2017121012:TMP:1.5 m above ground:1 hour fcst:
1.16:0:d=2017121012:RH:1.5 m above ground:1 hour fcst:
….
1.188:0:d=2017121012:TCDC:surface:15 hour fcst:
1.189:0:d=2017121012:APCP:surface:14-15 hour acc fcst:
1.190:0:d=2017121012:DSWRF:surface:14-15 hour ave fcst:

ってな感じで解析されます。コロン（：）が区切りで、左側から
　ID: ? : d=基準日時 : 気象項目 : 条件 : 初期値(anl) or 予測(fcst) :

気象項目は以下 　参考） http://www.hysk.sakura.ne.jp/data_list/JRA25_memo

PRMSL: Pressure reduced to MSL 海面更正気圧
PRES: Pressure 地上気圧
UGRD: u wind 東西風
VGRD: v wind 南北風
TMP: Temperature 気温
RH: Relative humidity 相対湿度
LCDC: Low level cloud cover 下層雲量
MCDC: Mid level cloud cover 中層雲量
HCDC: High level cloud cover 上層雲量
TCDC: Total cloud cover 全雲量
APCP: Convective precipitation ?降水量?
DSWRF: Downward short wave flux 日射量（下向き短波放射フラックス）

条件は

mean sea level 平均海面
surface 地面
10m above ground 地上10m
1.5m above ground 地上1.5m

上記のファイルは190項目のデータがある、ということがわかります。
例えば、気温の初期値が見たい場合は、気温の気象項目は TMP で、初期値は anl なので、ID は 1.5 になります。

$ wgrib2 filename -d 1.5 -csv outfile.csv

値が見たい場合は上記のコマンドで書き出せます。結構な量になるので注意しましょう。
メソ数値予測モデルでは 北緯22.4度 東経120度 から 北緯47.6度 東経150度 までを 緯度方向 505点 × 経度方向 481点 のメッシュに分けて予測しているので、全部で 242905点の予測地点（レコード）になります。
　→ http://www.data.jma.go.jp/add/suishin/jyouhou/pdf/205.pdf

Ruby から利用する

ここからが本題です。
上記でも計算しましたが、なにしろ一つの気象項目でも点数が多いので大変です。
とりあえず、一つの気象項目に絞って処理をします。

前段処理としてバイナリデータに変換します。（ダイレクトアクセスバイナリ形式）

$ wgrib2 filename -d ID -no_header -ieee outfile.bin
＃ファイル名は考えましょう。 → yyyymmddHH_ID_[anl or HHfcst].bin とか

これをすると範囲の左下（東経120度、北緯22.4度）から右上（東経150度、北緯47.6度）に向かって予測ポイントの値のみが羅列されたデータ（単精度浮動小数点：4bytes/1点）として保存されます。
なので、左下を起点に座標をステップ間隔（経度方向に 0.0625度、緯度方向に 0.05度の範囲）で計算すると何ブロック目かわかるので、それで値を取ってきて使用することができます。

例）松江城（東経 133.050746、北緯 35.475115）の場合
　　　(133.050746 - 120) / 0.0625 = 208
　　　(35.475115 - 22.4) / 0.05 = 261

　　　ブロック番号 = 208 + (481 * 261) = 125749番目
　　　ファイルの先頭から 125749 * 4 バイト目 の値（foat） -> 282.898

というわけで、この流れを ruby で書いてみる。

ソースコード（サンプル）

今回は ruby2.4.1 を使用しました。

grib2test.rb

# -*- coding: utf-8 -*-
require "open3"

class Grib2

  CRRPATH = File.expand_path("../", __FILE__) + "/"
  VARPATH = CRRPATH + "var/"
  TMPPATH = CRRPATH + "tmp/"
  
  CONF = {
    "surf" => { "lat" => { "min" => 120, "max" => 150,
                           "res" => 0.0625, "num" => 481 },
                "lng" => { "min" => 22.4, "max" => 47.6,
                           "res" => 0.05, "num" => 505 } },
    "pall" => { "lat" => { "min" => 120, "max" => 150,
                           "res" => 0.125, "num" => 241 },
                "lng" => { "min" => 22.4, "max" => 47.6,
                           "res" => 0.1, "num" => 253 } }  
  }

  def initialize
    Dir.mkdir(VARPATH) unless Dir.exist?(VARPATH)
    Dir.mkdir(TMPPATH) unless Dir.exist?(TMPPATH)
  end
  
  # --------------------------------------------------------------------------
  # datestr -> '20171205'
  # inittime -> 0, 3, 6, 9, 12, 15, 18, 21
  # modeltype -> 'surf' or 'pall'
  # forecasttime -> 0 - 39  (0: nal, 1-39: forecast)
  def create_filename(datestr, inittime, modeltype, forecasttime)
    fhead = "Z__C_RJTD_"
    fmodel = "_MSM_GPV_Rjp_L"
    filename = fhead + datestr + ("%02d" % inittime.to_i) + "0000" + fmodel
    mstr = "surf"
    mstr = "-pall" if modeltype == "pall"
    fstr = "_FH00-15"
    fstr= "_FH34-39" if forecasttime > 33
    fstr= "_FH16-33" if forecasttime > 15
    filename += mstr + fstr + "_grib2.bin"
    return filename
  end

  # modeltype -> 'surf' or 'pall'
  # lat, lng -> 133.05, 35.47
  def calculate_fposition(modeltype, lat, lng)
    mcnf = CONF[modeltype]
    latnum = ((lat - mcnf["lat"]["min"]) / mcnf["lat"]["res"]).to_i
    lngnum = ((lng - mcnf["lng"]["min"]) / mcnf["lng"]["res"]).to_i
    blocknum = latnum + (mcnf["lat"]["num"] * lngnum)
    fpos = blocknum * 4
    return fpos
  end

  # itemname -> "PRMSL", "PRES", "UGRD", "VGRD", "TMP", "RH"
  #             "LCDC", "MCDC", "HCDC", "TCDC"
  #             "APCP", "DSWRF"
  # filename -> "Z__C_RJTD_....._grib2.bin"
  # forecasttime -> 0 - 39 (0: nal, 1-39: forecast)
  def get_surf_data_id(itemname, filename, forecasttime)
    Open3.pipeline_r("wgrib2 #{VARPATH}#{filename}",
                     "grep #{itemname}"){|so, sss|
      if sss.last.value.to_i.zero?
        sos = so.read.split("\n")
        sos.each do |s|
          prm = s.split(":")
          prmfcst = prm[5].scan(/(\d{1,2}) hour/).flatten.last.to_i
          return prm[0] if prmfcst == forecasttime
        end
      else
        puts "can not open grib2.bin..."
      end
    }
    return nil
  end

  # itemname -> "HGT", "UGRD", "VGRD", "TMP", "VVEL", "RH"
  # filename -> "Z__C_RJTD_....._grib2.bin"
  # forecasttime -> 0 - 39 step by 3 hour (0: nal, 1-39: forecast)
  # pressur -> 1000, 975, 950, 925, 900, 850, 800, 700,
  #            600, 500, 400, 300, 200, 150, 100
  def get_pall_data_id(itemname, filename, forecasttime, pressur)
    Open3.pipeline_r("wgrib2 #{VARPATH}#{filename}",
                     "grep #{itemname}"){ |so, sss|
      if sss.last.value.to_i.zero?
        sos = so.read.split("\n")
        sos.each do |s|
          prm = s.split(":")
          prmfcst = prm[5].scan(/(\d{1,2}) hour/).flatten.last.to_i
          prmpres = prm[4].scan(/(\d{3,4}) mb/).flatten.last.to_i
          return prm[0] if prmfcst == forecasttime && prmpres == pressur 
        end
      else
        puts "can not open grib2.bin..."
      end
    }
    return nil
  end

  # data_id -> "1.12"
  # lat, lng -> 133.05, 35.4
  # filename -> "Z__C_RJT_....._grib2.bin"
  # fpos -> 1234
  def get_data(data_id, lat, lng, filename, fpos)
    tmpfilename = "#{$0}_" + Time.now.to_i.to_s + "_#{$$}.bin"
    val = nil
    so, se, ss = Open3.capture3("wgrib2 #{VARPATH}#{filename} -d #{data_id} -no_header -ieee #{TMPPATH}#{tmpfilename}")
    if ss.to_i.zero?
      File.open("#{TMPPATH}#{tmpfilename}", "rb") do |f|
        f.pos = fpos
        stval = f.read(4)
        val = stval.unpack("g")[0]
      end
      File.unlink("#{TMPPATH}#{tmpfilename}")
    else
      puts "wgrib2 parse error..."
    end
    return val
  end
  
end


# --------------------------------
# run

grib2 = Grib2.new

datestr = "20171205"
inittime = 0
forecasttime = 3

lat = 133.050746
lng = 35.475115

# 地上面の情報 → Z__C_RJTD_20171205000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin
modeltype = "surf"
itemname = "TMP"
filename = grib2.create_filename(datestr, inittime, modeltype, forecasttime)
gid = grib2.get_surf_data_id(itemname, filename, forecasttime)
fpos = grib2.calculate_fposition(modeltype, lat, lng)
val = grib2.get_data(gid, lat, lng, filename, fpos)
p "basetime: #{datestr} modeltype: #{modeltype}"
p "lat: #{lat} lng: #{lng} #{itemname} forecast: #{forecasttime} -> #{val}"

# 気圧面の情報 → Z__C_RJTD_20171205000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin
modeltype = "pall"
itemname = "TMP"
pressure = 850
filename = grib2.create_filename(datestr, inittime, modeltype, forecasttime)
gid = grib2.get_pall_data_id(itemname, filename, forecasttime, pressure)
fpos = grib2.calculate_fposition(modeltype, lat, lng)
val = grib2.get_data(gid, lat, lng, filename, fpos)
p "basetime: #{datestr} modeltype: #{modeltype} pres: #{pressure}"
p "lat: #{lat} lng: #{lng} #{itemname} forecast: #{forecasttime} -> #{val}"

気象業務支援センターからダウンロードしたデータは ./var/ ディレクトリ以下に入れておきます。
＃気圧面の情報も取得しておいてください。

で、実行すると

$ ruby gribtest.rb
"basetime: 20171205 modeltype: surf"
"lat: 133.050746 lng: 35.475115 TMP forecast: 3 -> 278.688720703125"
"basetime: 20171205 modeltype: pall pres: 850"
"lat: 133.050746 lng: 35.475115 TMP forecast: 3 -> 265.2142333984375"

な感じです。
あとは各項目で単位があるので、わかりやすいように計算して表示すればいいかと。