How to make/use Station Data †

#ref(): File not found: "amdcomp_wsp_197601.gif" at page "IT memo/GrADS memo/Station data"

↑

Station Data の作り方 †

  • データファイルの構造
    • 同じ時刻のデータについて、1レコード（28バイトヘッダ＋4バイトデータ×要素数）を繰り返す。
      4byte blankstidが入っても良い(?)
      4byte ascii: stidstidは8文字まで。
      4byte real (*2): lat, lon
      4byte real (*1) and integer (*2): time, nlev, nflagそれぞれ通常0.,1,1(?)
      4byte real (*X): rval30個までは動作確認
    • 時刻が変わるときは、rvalなしの空ヘッダ(28バイト)を差し込む。その際nlev=0にする。
    • 上記の繰り返し（最後も28バイトの空ヘッダ←重要）

fortranを使う場合の注意：

  • open文の属性でrecordtype="stream"を使うと便利。しかし、少なくともPGFとGNUは非対応。ifortはOK（コンパイル時の-assume bytereclを忘れずに）

その他注意：

  • コントロールファイルにZDEFが設定できないため、鉛直方向にレコードを切ることができない。時間方向のみ。したがって、変数を増やして対応するのが良い（1000hPaのT、500hPaのTなどを同じレコードに入れる）。（改善してほしい点）
↑

Station Data の読み方 †

  • GRIBデータを読む際、ctlファイルからgribmapを通してidxファイルを作成するように、まず、ctlファイルを作成する。通常のGRIBグリッドデータと違う点は、以下の二項目。 
    dtype station
stnmap XXX.map
  • 次に、stnmapの実行 
    $ stnmap -i XXX.ctl
  • GrADSを起動して、ctlを開く 
    $ grads
> open XXX.ctl
> set digsiz 0.3
> set stid on
> d p
↑

BUFRファイル †

GrADS本家のページによるとBUFRファイルを読み込んで描画できるとのこと。しかし、SIO/GDASのサンプルパッケージについているprepbufr78120100.bufr（いわゆるNCEP PREPBUFRファイル）について、 bufrscanで中身を読み取り、bufr.ctlに内容を書いてみるものの、うまくいかなかった。

>open bufr.ctl
>d temp
Parsing BUFR file prepbufr78120100.bufr
Finished parsing BUFR file

などとでるが、図は出力されず。gxout statにしてみると、データをひとつとして読み取っていないようだ。

上のprepbufrをGrADSで読めるツワモノの人、助けて…。参考までに作ったCTLファイルを下に貼ります。それか、GrADSで読めるBUFRファイルとCTLファイルの組み合わせをコピーさせてくれるだけでも非常にありがたいです。お願いします。

DSET ^prepbufr78120100.bufr
TITLE prepbufr78120100.bufr
UNDEF -9999
DTYPE bufr
OPTIONS big_endian
TDEF 6 linear 0Z1DEC1978 1hr
STID 1,192
XVAR 6,2
YVAR 5,2
VARS 13
slon 0 6,2 Station longitude
slat 0 5,2 Station latitude
selv 0 10,194 Station elevation
pres 0 7,192 Pressure
spfh 0 13,192 Specific humidity
dewpt 0 12,194 Dew point temperature
temp 0 12,192 Temperature
hgt 0 10,195 Height
uwnd 0 11,3 U-component wind
vwnd 0 11,4 V-component wind
wdir 0 11,1 Wind direction
wspd 0 11,191 Wind speed
tpw 0 13,221 Total precipitable water
ENDVARS

（2009/3/4追記）少しだけ進展。別のprepbufr（NCEP 2000/01/01 0Zのもの）を同じコントロールファイルで開いてみると、座標情報と観測属性だけは図化することができた。どうも、Replication Factorという値が-001のものだけは読み込めるようだ。肝心の観測値はReplication Factorがすべて000以降なので、まったく読めず。もう、prepbufrのGrADSによる直接図化は諦めて、Hwangさんがつくったpickupというプログラムを使って書き出すことにする。

↑

描画のコツ †

↑

Station IDの記入 †

> set stid on
> d var
↑

小数点以下の桁数の変更 †

> set dignum 1 ## 小数点以下1桁
> d var
↑

時系列変化 †

任意の一地点を指定した上で、d var(stid=stationname)など。

> set lat 1
> set lon 1
> set t 1 2
> d p(stid=QQQ)
↑

内挿(oacres) †

Reference用に緯度・経度方向に等間隔のデータを用意し、oacres関数を用いる。

> open grid.ctl # 等間隔の格子データ（例えば2.5度再解析大気データなど。）
> open stn.ctl # Station data
> d oacres(gridvar.1,stnvar.2) # gridvarの格子で、変数stnvarを内挿

oacres内挿後、演算可能。従ってRMSなども楽に計算できる。散布図も通常と同様に描ける。

> d sqrt(aave(pow(oacres(gridvar.1,stnvar.2)-gridvar.1,2),lon=-180,lon=180,lat=-90,lat=90))

しかし、内挿後のデータは面的に広がるため、一対一の比較にはならない。できるだけポイントでの比較をしたいときは、

oacres(gridvar.1,stnvar.2,3,1,0.5)-gridvar.1

などとする。要するに、内挿の影響範囲をできるだけ小さくしている。下のgr2stnを用いるとほとんどの空間平均の演算（aaveなど）や散布図の描画ができなくなるため、こうせざるを得ないようだ。（2008/07/17現在。誰かstationデータ使った散布図の描き方、教えてください。）　

↑

内挿(gr2stn) †

oacresとは逆に、grid格子のデータをstn格子点に内挿する。RMSを計算するにはこちらのほうが（おそらく）正確。しかし、一点一点計算するため、遅い。

> d gr2stn(gridvar.1,stnvar.2)

どうも、station dataのままでは、空間平均の算出や散布図描画ができない模様。一番使いたい機能な気がするのだが…。

↑

アメダス年報データの取り扱い †

基本的に、下記リンクにある名大中村さんの言うとおりにやる。ただ、すべての観測点のデータをひとつのStation Dataファイルにまとめ、CTLファイルにoptions templateを記載した。

↑

Link †