======================
sp.gs
======================
header='ERA-I_110118-110125_Barents_Sfc'
fac=100
varname='sp'
outdir='Fig'
'sdfopen /work4/data/ERA-Interim/' header '.nc'
'set mproj nps'
'set lat 30 90'
'set lon -90 90'
i=1
while (i <= 32)
'set t ' i
var=sp
'sclvar=' var '/100'
'cc'
'set mpdset mres'
'set gxout shade2'
'color -kind dodgerblue->white -levs 0 0.4 0.9 -xcbar 3 5 2.4 2.6 -fh 0.12 -fw 0.11 -edge triangle -line on'
'd ci'
'basemap L'
levs1='960 965 970 975 980 985 990 995 1000 1005 1010 1015 1020 1025 1030'
'set gxout contour'
'set clevs 'levs1
'set clab off'
'set cthick 10'
'set ccolor 0'
'd msl/100'
'color -gxout contour -kind red->red->orangered->orange->gold->navajowhite -levs ' levs1 ' -xcbar 0.5 8 3 3.2 -fh 0.12 -fw 0.11 -edge triangle -line on'
'set cthick 5'
# 'd maskout(sp/100,sst)'
'd msl/100'
'set clevs 990 1000'
'set clab forced'
'set cstyle 0'
'set ccolor 0'
'set cthick 0'
'set clopts 1 3 0.08'
'd msl/100'
'circlat 20 85 0.1'
'circlon 10 0.1'
'q dims'
# say result
line=sublin(result,5)
datetime=subwrd(line,6)
num=''
if (i<10)
num='00'i
endif
if (i>=10 & i<100)
num='0'i
endif
if (i>100)
num=i
endif
'draw title 'datetime
figfile=outdir '/' header '_' varname '' num '.png'
say 'Output: ' figfile
'gxprint ' figfile
i=i+1
endwhile
'allclose'
----------------------
End of sp.gs
----------------------
======================
circlat.gs
======================
function circlat(args)
* usage: run circlat dlat (strsiz (font))
* mproj=nps もしくは sps において 緯度の数字を描く
* dlat は数字を描く経度の間隔で必須.
* lon は数字を描く経度で必須.
* strsiz は文字の大きさで必須.
* font はフォント番号でオプション
dlat=subwrd(args,1)
if (dlat=0); say '*circlat* dlat=0'; 'quit'; return; endif
lon0=subwrd(args,2)
strsiz=subwrd(args,3)
if (strsiz=""); say '*circlat* strsiz is not defined'; 'quit'; return; endif
if (strsiz=0); say '*circlat* strsiz=0'; 'quit'; return; endif
font=subwrd(args,4)
'q gxinfo'
cline=sublin(result,6)
mproj=subwrd(cline,3)
if (mproj!=3 & mproj!=4)
say ' *circlon* (mproj!=3:nps & mproj !=4:sps) is not supported'
say 'mproj='%mproj
'quit'
endif
if (mproj=3);* mproj=3 is nps
'query w2xy 0 90'
* say 'result='result
endif
if (mproj=4);* mproj=4 is sps
'query w2xy 0 -90'
endif
xpol=subwrd(result,3)
ypol=subwrd(result,6)
'q dims'
lonline=sublin(result,2)
lon1=subwrd(lonline,6)
lon2=subwrd(lonline,8)
latline=sublin(result,3)
lat1=subwrd(latline,6)
lat2=subwrd(latline,8)
lat=math_int(lat1/dlat)*dlat-dlat
'set line 0';* set color for recf
'set string 1 c'
'set strsiz '%strsiz
if (font!=""); ' set font '%font; endif
while(lat<=lat2)
if (lat1<=lat)
'query w2xy '%lon0%' '%lat
x=subwrd(result,3)
y=subwrd(result,6)
if (lat<0); clat=-lat%'S'; endif
if (lat>0); clat=lat%'N'; endif
if (lat=0); clat='Eq'; endif
len=math_strlen(clat)
xwl=strsiz*len*1.1/2
xwr=strsiz*len*1.3/2
ywb=strsiz*1.4/2
ywt=strsiz*1.2/2
'draw recf '%x-xwl%' '%y-ywb%' '%x+xwr%' '%y+ywt
'draw string '%x%' '%y%' '%clat
endif
lat=lat+dlat
endwhile
'set string 1 bl 3 0'
return
----------------------
End of circlat.gs
----------------------
======================
circlon.gs
======================
function circlon(args)
* usage: run circlon dlon (strsiz (font))
* mproj=nps もしくは sps において longitudeの数字を描く
* dlon は数字を描く経度の間隔で必須.
* strsiz は文字の大きさでオプション
* font はフォント番号でオプション
dlon=subwrd(args,1)
strsiz=subwrd(args,2)
font=subwrd(args,3)
'q gxinfo'
cline=sublin(result,6)
mproj=subwrd(cline,3)
if (mproj!=3 & mproj!=4)
say ' *circlon* (mproj!=3:nps & mproj !=4:sps) is not supported'
say 'mproj='%mproj
'quit'
endif
if (mproj=3);* mproj=3 is nps
'query w2xy 0 90'
* say 'result='result
endif
if (mproj=4);* mproj=4 is sps
'query w2xy 0 -90'
endif
xpol=subwrd(result,3)
ypol=subwrd(result,6)
'q dims'
lonline=sublin(result,2)
latline=sublin(result,3)
if (mproj=3);* mproj=nps
lat0=subwrd(latline,6)
else;* mproj=sps
lat0=subwrd(latline,8)
endif
lon1=subwrd(lonline,6)
lon2=subwrd(lonline,8)
'query w2xy '%lon1%' '%lat0
xd=subwrd(result,3)
yd=subwrd(result,6)
say 'xpol='xpol
say 'ypol='ypol
say 'xd='xd
say 'yd='yd
* 'run sqrt '%(xpol-xd)*(xpol-xd)+(ypol-yd)*(ypol-yd); rad=result;* for Grads 1.7 or earlier
rad=math_pow((xpol-xd)*(xpol-xd)+(ypol-yd)*(ypol-yd),0.5);
* say 'rad='rad
* say 'lon1='lon1' lon2='lon2
radd=rad+0.2
if (strsiz!="")
radd=rad+1.5*strsiz
'set strsiz '%strsiz
endif
if (font!=""); ' set font '%font; endif
lonm=(lon1+lon2)/2
lon=0; while (lon<360)
if (lon1<=lon & lon<=lon2)
* say 'lon='lon
anglestr=lon-lonm
angleline=(anglestr-90)/180*3.14
* 'set string 1 c 3 '%anglestr
'set string 1 c'
cosangle=math_cos(angleline);
sinangle=math_sin(angleline);
* 'run cos '%angleline; cosangle=result;* for Grads 1.7 or earlier
* 'run sin '%angleline; sinangle=result;* for Grads 1.7 or earlier
xrim=xpol+rad*cosangle
yrim=ypol+rad*sinangle
* 'draw line '%xpol%' '%ypol%' '%xrim%' '%yrim
xchar=xpol+radd*cosangle
ychar=ypol+radd*sinangle
ew=''
if (0<lon & lon<180); ew='E'; endif
if (180<lon & lon<360); ew='W'; endif
if (lon<180)
clon=lon%ew
else
clon=(360-lon)%ew
endif
'draw string '%xchar%' '%ychar%' '%clon
endif
lon=lon+dlon
endwhile
'set string 1 bl 3 0'
return
lat=lat1; while (lat<=lat2)
endwhile
----------------------
End of circlon.gs
----------------------