!!このページはただのメモなので、未確認事項が含まれています。ご了承下さい。
現在の海洋モデルの多くは極問題を避けるために、緯度経度座標系以外の座標系を採用している。プロダクトとしてのデータも緯度経度座標系でないことが多い。例えばCMIP5の海水温や海面水温データは
MIROC5 海洋モデルCOCO4.5によるSST出力データ
MPI-ESM-LR 海洋モデルMPIOMによるSST出力データ
このような座標系で提供されている。これらをGrADSで解析するには、緯度経度との対応関係を用いて緯度経度座標系に変換すればよい。
Step 1. NetCDFファイルから必要な情報を抜き出す
コントロールファイルMIROC5.ctlを作成、変数とlon、latの(x, y)配列データを利用する
SSTを色で、緯度経度を等値線で描く
上の図ができる。
ctlファイルを用いて
の三つのデータを出力
Step 2.緯度経度座標系に変換
lont, lati, tosの3列データを緯度経度座標系に変換する。
変換方法は何でも良いが、GMTコマンドが使用可能であればnearneighborを使うと楽。
lont, lati, tosの3列asciiデータを作成
※欠損値の処理をすること
nearneighborで1度×1度の緯度経度座標系グリッドデータに変換
ヘッダを読み飛ばして、コントロールファイルで読み込み、作図する。
緯度経度座標系のSSTデータができる(下図はSST上昇量の例)。nearneighborのオプションを使って内挿方法を調整すること。
衛星から海上風を測っている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形式に変換。オリジナルを解凍したり圧縮したりするので時間がかかる。
できたNetCDFファイルを読むctlファイルを作成する。
試したことがないので、詳しそうなリンクをはっておく。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
現在の海洋モデルの多くは極問題を避けるために、緯度経度座標系以外の座標系を採用している。プロダクトとしてのデータも緯度経度座標系でないことが多い。例えば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
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と地点データは相性が良いとは言えませんが…
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
コメントをかく