memo - GrADS-Note - Seesaa Wiki（ウィキ）

！！このページはただのメモなので、未確認事項が含まれています。ご了承下さい。

非緯度経度座標系のデータの変換

現在の海洋モデルの多くは極問題を避けるために、緯度経度座標系以外の座標系を採用している。プロダクトとしてのデータも緯度経度座標系でないことが多い。例えばCMIP5の海水温や海面水温データは

MIROC5 海洋モデルCOCO4.5によるSST出力データ

MPI-ESM-LR 海洋モデルMPIOMによるSST出力データ

このような座標系で提供されている。これらをGrADSで解析するには、緯度経度との対応関係を用いて緯度経度座標系に変換すればよい。

Step 1. NetCDFファイルから必要な情報を抜き出す

コントロールファイルMIROC5.ctlを作成、変数とlon、latの(x, y)配列データを利用する

DSET    ^tos_Omon_MIROC5_piControl_r1i1p1_230001-239912.nc
DTYPE   netcdf
TITLE   XXX
UNDEF   0
OPTIONS 365_day_calendar
xdef 256 linear 0 1.40625   #適当
ydef 224 linear -89.598214286 0.803571429   #適当
zdef 1 levels 0
tdef 1200 linear 12Z16JAN2300 1mo
VARS    3
tos=>tos 1 t,y,x Temperature
lon=>lont 1 y,x Temperature   #lonとlatを予約済み変数以外の名前に置き換える
lat=>lati 1 y,x Temperature
ENDVARS

SSTを色で、緯度経度を等値線で描く

'open MIROC5.ctl '
'set x 1 256'
'set y 1 224'
'cc'
'sh'
'set mpdraw off'
'set xlab off'
'set ylab off'
'set grads off'
'set grid off'
'set gxout grfill'
'set parea 1 8 1 6'
'd tos'
'xcbar 1 8 0.6 0.75 -edge triangle -line on'
'co'
'set cthick 5'
'set cstyle 1'
'set ccolor 1'
'set clab on'
'set cint 60'
'd lont'
'set cthick 5'
'set cstyle 1'
'set ccolor 1'
'set clab on'
'set cint 20'
'd lati'

上の図ができる。
ctlファイルを用いて

 tos(x, y)
 lont(x, y)
 lati(x, y)

の三つのデータを出力

Step 2.緯度経度座標系に変換

lont, lati, tosの3列データを緯度経度座標系に変換する。
変換方法は何でも良いが、GMTコマンドが使用可能であればnearneighborを使うと楽。

lont, lati, tosの3列asciiデータを作成

program main
implicit none

integer,parameter:: nx=256,ny=224
real(4):: lon(nx,ny),lat(nx,ny),a(nx,ny)
integer:: ix,iy

open(10,file='anom.bin', &
       access='direct',form='unformatted',recl=4*nx*ny)
open(11,file='lon.bin', &
       access='direct',form='unformatted',recl=4*nx*ny)
open(12,file='lat.bin', &
       access='direct',form='unformatted',recl=4*nx*ny)

open(20,file='data1.txt', &
       access='sequential',form='formatted')

read(10,rec=1) ((a(ix,iy),ix=1,nx),iy=1,ny)
read(11,rec=1) ((lon(ix,iy),ix=1,nx),iy=1,ny)
read(12,rec=1) ((lat(ix,iy),ix=1,nx),iy=1,ny)

do iy=1,ny
do ix=1,nx
write(20,*) lon(ix,iy),lat(ix,iy),a(ix,iy)
enddo
enddo

end

※欠損値の処理をすること

nearneighborで1度×1度の緯度経度座標系グリッドデータに変換

nearneighbor data1.txt -Gdata1.bin -I1 -R0/360/-89.5/89.5 -S4 -N8/7 -fg

ヘッダを読み飛ばして、コントロールファイルで読み込み、作図する。

dset ^data1.bin
title XXX
FILEHEADER 5012
undef -9.99e+08
options big_endian
xdef 361 linear 0 1
ydef 178 linear -89.5 1
zdef 1 levels 0
tdef 1 linear 12Z16JAN1979 1mo
vars 1
sst 1 -999 TOA Outgoing Shortwave Radiation
endvars

緯度経度座標系のSSTデータができる（下図はSST上昇量の例）。nearneighborのオプションを使って内挿方法を調整すること。

QuikSCAT(hdfファイル)をGrADSで読む

衛星から海上風を測っているQuikSCAT
http://podaac.jpl.nasa.gov/PRODUCTS/p109.html
というデータを度々見かけて気になるので、GrADSで絵を描いてみようと思い立ったものの、
なかなか一筋縄では行かなかった。

・QuikSCATというデータセット固有の問題
・hdf形式をGrADSで扱うという問題

のおそらく二つ目が大きな障壁だった。

hdf形式のデータは一応GrADSで扱えるようなのだが(1.9以前ではgradshdf)、HDF4
http://www.spherewind.com/id-1174316400.html#hdf
がインストールされていないと無理？なのか、
QuikSCATのhdfファイルに緯度経度情報が含まれていないのが問題なのか、
適切なコントロールファイルを用意していないことが問題なのか、
GrADSでそのまま開くことはできず。

次にQuikSCATの提供データに含まれている、hdfを読むためのソフトウェア
ftp://podaac.jpl.nasa.gov/pub/ocean_wind/quikscat/...

を使ってバイナリ変換を目論むが、
ifortでread_qscat3.fがうまく動かない。
そもそもmakeについて理解が浅いことが問題なのか…

結局、hdfをNetCDFに変換するユーティリティ、hdf2netcdf
ftp://podaac-ftp.jpl.nasa.gov/OceanWinds/quikscat/
を使って(COARSE規約に従っていない)NetCDFファイルに変換し、
特別なctlファイルを使って読む方法をとった。

http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%2...

解決策

ftp://podaac-ftp.jpl.nasa.gov/OceanWinds/quikscat/

から L3/README.txt, data/, doc/ と hdf2netcdf/ をダウンロード。
（jpl、swなどのサブディレクトリに移動しているかも）
以下のシェルスクリプトを使って、hdf形式のL3/の中身をNetCDF形式に変換。オリジナルを解凍したり圧縮したりするので時間がかかる。

 #!/bin/bash
 
 qsdir=・・・/QSCAT
 hdf2ncdir=${qsdir}/hdf2netcdf
 ctgr=L3
 indir=${qsdir}/hdf/${ctgr}/data
 outdir=${qsdir}/nc
 ncdump=${hdf2ncdir}/ncdump/Linux2.4/ncdump
 ncgen=${hdf2ncdir}/ncgen/linux_2.4-i686/ncgen 
 mkdir -p ${outdir}

 for year in 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 ;do 
    filegz_list=`ls ${indir}/${year}/QS_XWGRD3_${year}*.gz`
    for filegz in ${filegz_list} ;do
       echo "gunzip ${filegz} ..."
       gunzip ${filegz}
    done #filegz 
    file_list=`ls ${indir}/${year}/QS_XWGRD3_${year}*`
    for file in ${file_list} ;do
       doy=`expr ${file} : ".*QS_XWGRD3_${year}\(...\)."`
       ${ncdump} ${file} | ${ncgen} -b -o ${outdir}/${year}_${doy}.nc
       echo  'infile=' ${file}
       echo  'outfile=' ${outdir}/${year}_${doy}.nc
       echo "gzip ${file} ..."
       gzip ${file}
    done #file
 done #year

できたNetCDFファイルを読むctlファイルを作成する。

 DSET ^nc/1999_%ch.nc

 dtype netcdf
 TITLE QuikSCAT Level 3 Ocean Wind Vectors in 25km grid
 UNDEF -999
 UNPACK scale_factor add_offset
 options template
 chsub 1 1 001
 chsub 2 2 002
 ・・・
 chsub 365 365 365
 *chsub 366 366 366
 XDEF 1440 linear 0.0 0.25
 YDEF  720 linear -90.0 0.25
 ZDEF 1 levels 0.0
 TDEF 365 linear 1jan1999 1dy
 VARS 16
 asc_avg_wind_speed=>aaws     0 y,x wind speed for the ascending pass [m/s]
 des_avg_wind_speed=>daws     0 y,x wind speed for the descending pass [m/s]
 asc_avg_wind_vel_u=>aawvu    0 y,x u-component of the wind velocity [m/s] 
 des_avg_wind_vel_u=>dawvu    0 y,x u-component of the wind velocity [m/s]
 asc_avg_wind_vel_v=>aawvv    0 y,x v-component of the wind velocity [m/s] 
 des_avg_wind_vel_v=>dawvv    0 y,x v-component of the wind velocity [m/s]
 asc_avg_wind_speed_sq=>aawss 0 y,x square of the wind speed [m2/s2]
 des_avg_wind_speed_sq=>dawss 0 y,x square of the wind speed [m2/s2]
 asc_wvc_count=>awc           0 y,x distinguishes NULL values from zero value  of wind speed [count]
 des_wvc_count=>dwc           0 y,x distinguishes NULL values from zero value of wind speed [count]
 asc_time_frac=>atf           0 y,x time in fraction of a day [fraction of day]
 des_time_frac=>dtf           0 y,x time in fraction of a day [fraction of day]
 asc_rain_prob=>arp           0 y,x probability of columnar rain rate that is greater than 2 km*mm/hr [n/a]
 des_rain_prob=>drp           0 y,x probability of columnar rain rate that is greater than 2 km*mm/hr [n/a]
 asc_rain_flag=>arf           0 y,x the rain flag algorithm detected the presence of rain [n/a]
 des_rain_flag=>drf           0 y,x the rain flag algorithm detected the presence of rain [n/a]
 ENDVARS

GrADSステーションデータ (station data) の取り扱い

試したことがないので、詳しそうなリンクをはっておく。GrADSと地点データは相性が良いとは言えませんが…

http://www.iges.org/grads/gadoc/aboutstationdata.h...

http://hydro.iis.u-tokyo.ac.jp/~kei/?IT%20memo%2FG...

http://mc39.genin.jp/GrADS_Station_Data.html

http://www.rain.hyarc.nagoya-u.ac.jp/laboratory/OB...

http://boa-sorte.net/grads_station_file.html

http://www28.atwiki.jp/fmemo/pages/20.html

http://nattoumaruo.com/blog/archives/32