Gradsでregrid

GrADSにおいて、二つのデータセットのグリッドをそろえて、差を見てみることをおこなった。

OpenGrADSのre関数を用いた。
下のソース例、
define sstre=re(sst.1,180,linear,0,2,89,linear,-88,2,ba)
では、変数sst.1を
経度方向に、180グリッド、0度から2度刻みで、線形
緯度方向に、89グリッド、南緯88度から2度刻みで、線形
box averageでグリッドに落としている。

例ではNOAA OISST version 2(1)とNOAA Extended Reconstructed SST version 3b(2)という二つの１月の海面水温データセットから２０１１年１月の値をそれぞれプロットし、3番目の図で、その差(2-1)を引いている。データはNOAAのサイトからOPeNDAPで取得している。

下がgradsのスクリプト。
"grads -p"でportraitモードで使う。
他に、mul.gs, color.gs, draws.gs という Chihiro Kodama氏のgradsスクリプトを用いている。

"reinit"
*NOAA OISST V2
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.oisst.v2/sst.m
nmean.nc"
"set time jan2011"
"mul 1 3 1 3 -yoffset -0.5"
"set grads off"
"color 0 30 2 -kind grainbow -gxout grfill"
"d sst"
"set strsiz 0.15"
"draws 1) NOAA OISST V2"

*NOAA Extended Reconstructed SST V3b
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.ersst/sst.mnme
an.nc"
"set time jan2011"
"mul 1 3 1 2  -yoffset -0.5 "
"set grads off"
"color 0 30 2 -kind grainbow -gxout grfill"
"d sst.2"
"set strsiz 0.15"
"draws 2) NOAA Extended Reconstructed SST V3b"
"cbarn 0.7 1 7.6 7.5"

*diff
"define sstre=re(sst.1,180,linear,0,2,89,linear,-88,2,ba)"
"mul 1 3 1 1  -yoffset -0.5"
"set grads off"
"color -levs -2 -1.5 -1 -0.5 0.5 1.0 1.5 2  -kind blue->white->red -gxout grfill"
"d sst.2-sstre"
"set strsiz 0.15"
"draws 2)-1)"
"cbarn 0.7 1 1. 2."
"printim sst_regrid_compare.png white"

GaLab in Python Interface to GrADS

Python Interface to GrADSのGaLabを用いてプロットしてみた。
プロットしたデータはEarth System Research Laboratoryで配布されているNCEP/NCAR再解析の2011年１月の海水面気圧をOPeNDAPで取得している。

import matplotlib.pyplot as plt
from grads.galab import GaLab  # for python interface for grads

ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)

ga.open("http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/slp.mon.mean.nc")

ga.blue_marble("on")
ga("set time jan2011")
ga.basemap('ortho',opts=(130.,40.))

ga.contour("slp")
del ga
plt.title("Sea Level Pressure (Jan2011)")
plt.savefig("galab_slp.png")

Matplotlib Basemap Toolkitよりも手軽にプロットできる半面、緯度・経度線を取り除いたりとかの自由はききにくいようだ。