- 积分
- 28
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2019-1-17
- 最后登录
- 1970-1-1
|
发表于 2019-3-10 13:04:11
|
显示全部楼层
我想画的图很简单，就是850hPa与海温温差和海温的叠加图，海温我画等值线，温差画阴影，这样的话，陆地的数据就没有意义，容易引起误会，需要屏蔽。具体的gs文件如下：
*------------------------------------------------
* T850-SKT
*------------------------------------------------
'set dfile 1'
'set lon 70 150'
'set lat 10 60'
'c'
'run fang_map.gs'
'color -40 40 4 -kind blue->aqua->lightskyblue->white->khaki->orange->red '
* 'set cint 4'
'set cstyle 1'
'd th(lev=850)-sstk'
'cbar'
'color -gxout contour -levs -32 -28 -24 -20 -16 -12 -8 -4 0 4 8 12 16 20 24 28 32'
'set cthick 10'
'd sstk-273.15'
'basemap l 3 1'
'set string 1 bl 10'
'draw string 3.7 8.2 T850-SST(C) shadow & SST(C) line'
' set string 1 bl 10'
'draw string 1 7.8 Vaild time(UTC):' nt
' set string 1 bl 10'
'draw string 7 7.8 Initial time(UTC):'strfh'Z'
'printim ' 'E_domain1_sst_sfc_'strfh '_' ts'.png x1800 y1200 white '
basemap修改的如下：
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* basemap.gs
*
* This script overlays a land or ocean mask that exactly matches
* the coastal outlines.
*
* Usage:
* basemap L(and)/O(cean) <fill_color> <outline_color> <L(owres)/M(res)/H(ires)>'
*
* Example:
* ga-> set mpdset mres
* ga-> display sst
* ga-> basemap L 7 2 M
*
* Usage Notes:
*
* The default values for the optional arguments are:
* fill_color 15, outline_color 0, and lowres.
*
* The land and ocean masks are composed of hundreds of
* polygons that are specified in accompanying ascii files.
* The ascii files must be downloaded from the GrADS script library:
* ftp://grads.iges.org/grads/scripts/lpoly_lowres.asc
* ftp://grads.iges.org/grads/scripts/lpoly_mres.asc
* ftp://grads.iges.org/grads/scripts/lpoly_hires.asc
* ftp://grads.iges.org/grads/scripts/opoly_lowres.asc
* ftp://grads.iges.org/grads/scripts/opoly_mres.asc
* ftp://grads.iges.org/grads/scripts/opoly_hires.asc
*
* For the low and medium resolution map files, coverage is global.
* For the high resolution map, coverage is limited to North America
* (0-90N, 170W-10W).
*
* Basemap will work with any scaled or latlon map projection.
* If you are using Grads version 1.8 or higher, this script
* will also work properly with the robinson projection and
* polar stereographic projections from 0-90, 15-90, and 20-90
* (North and South). Other projections will work but are not
* guaranteed because GrADS may not clip the basemap properly.
* A solution to this problem is to use "set mpvals" to override
* the dimension environment limits. For example:
* set mproj nps
* set lon -180 180
* set lat 0 90
* set mpvals -180 180 60 90
* display <something>
* basemap L
* The resulting plot will be a properly clipped square.
*
* Special Feature:
* An additional option is to mask out the Mexican and
* Canadian land regions surrounding the US, so that only the
* conterminous states are seen. To use this feature, change
* your land polygon file from lpoly_lowres.asc to lpoly_US.asc:
* ftp://grads.iges.org/grads/scripts/lpoly_US.asc
* Run basemap twice:
* basemap o 0 0 (<- that's O zero zero) ;* mask out ocean
* basemap L 0 0 ;* mask out non-US land
* This will only work properly if your domain is within the
* boundaries 20N-50N, 130W-60W. Low-res maps only.
*
* Written by Jennifer M. Adams, jma@cola.iges.org
* Updated October 2003
*
* Special thanks to Andrew Schneider, Tom Schneider and Emily Straus
* for their painstaking efforts to create the polygon files.
*
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
function main(args)
* There are defaults for the colors and resolution
* User must specify which mask
if (args='')
say 'Usage: basemap L(and)/O(cean) <fill_color> <outline_color> <L(owres)/M(res)/H(ires)>'
return
else
type = subwrd(args,1)
color = subwrd(args,2)
outline = subwrd(args,3)
res = subwrd(args,4)
if (color = '') ; color = 15 ; endif
if (outline = '') ; outline = 0 ; endif
if (res = 'H' | res = 'h') ; hires=1 ; mres=0; lowres=0 ; endif
if (res = 'M' | res = 'm') ; hires=0 ; mres=1; lowres=0 ; endif
if (res = 'L' | res = 'l') ; hires=0 ; mres=0; lowres=1 ; endif
if (res = '') ; hires=0 ; mres=0; lowres=1 ; endif
endif
* Set the polygon data files
if (type = 'L' | type = 'l')
* if (lowres); file = 'lpoly_US.asc'; endif
if (lowres); file = 'C:/OpenGrADS/Contents/Resources/lpoly_lowres.asc'; endif
if (mres) ; file = 'C:/OpenGrADS/Contents/Resources/lpoly_mres.asc' ; endif
if (hires) ; file = 'C:/OpenGrADS/Contents/Resources/lpoly_hires.asc' ; endif
endif
if (type = 'O' | type = 'o')
if (lowres); file = 'C:/OpenGrADS/Contents/Resources/opoly_lowres.asc'; endif
if (mres) ; file = 'C:/OpenGrADS/Contents/Resources/opoly_mres.asc' ; endif
if (hires) ; file = 'C:/OpenGrADS/Contents/Resources/opoly_hires.asc' ; endif
endif
* Make sure there's a plot already drawn
'q gxinfo'
line5 = sublin(result,5)
line6 = sublin(result,6)
xaxis = subwrd(line5,3)
yaxis = subwrd(line5,6)
proj = subwrd(line6,3)
if (xaxis = 'None' | yaxis = 'None')
say 'You must display a variable before using basemap'
return
endif
* See what version of Grads is running
'q config'
line = sublin(result,1)
word = subwrd(line,2)
version = substr(word,2,3)
if (version >= 1.8)
newgrads = 1
else
newgrads = 0
endif
* See if map projection will be supported
if (newgrads = 0)
if (proj != 1 & proj != 2)
say 'Only scaled or latlon projections are supported with GrADS v'version
return
endif
endif
* Get the image edges for clipping
'q gxinfo'
line3 = sublin(result,3)
line4 = sublin(result,4)
x1 = subwrd(line3,4)
x2 = subwrd(line3,6)
y1 = subwrd(line4,4)
y2 = subwrd(line4,6)
'set clip 'x1' 'x2' 'y1' 'y2
* Read the first record from the polygon file
result = read(file)
rc = sublin(result,1)
rc = subwrd(rc,1)
if (rc!=0)
say 'Error reading 'file
return
endif
nwcmd = sublin(result,2)
* Read subsequent records, allowing for read input buffer overflow
flag = 1
while (flag)
ignore = 0
wcmd = nwcmd
while(1)
result = read(file)
rc = sublin(result,1)
rc = subwrd(rc,1)
if (rc!=0)
flag = 0
break
else
nwcmd = sublin(result,2)
if (subwrd(nwcmd,5) != 'draw')
wcmd = wcmd % nwcmd
else
break
endif
endif
endwhile
* Get the lat/lon range of the current dimension environment
'q dims'
line1 = sublin(result,2)
line2 = sublin(result,3)
minlon = subwrd(line1,6)
maxlon = subwrd(line1,8)
minlat = subwrd(line2,6)
maxlat = subwrd(line2,8)
* The range of the polygon is specified in the first four words of the record
minwx = subwrd(wcmd,1)
maxwx = subwrd(wcmd,2)
minwy = subwrd(wcmd,3)
maxwy = subwrd(wcmd,4)
* If the polygon is outside the current dimension, ignore it
if (minwx >= maxlon) ; ignore = 1 ; endif
if (maxwx <= minlon) ; ignore = 1 ; endif
if (minwy >= maxlat) ; ignore = 1 ; endif
if (maxwy <= minlat) ; ignore = 1 ; endif
if (!ignore)
count = 7
nvert = 1
if (newgrads)
cmd = 'draw mappoly '
else
cmd = 'draw polyf '
endif
while (1)
countx = count
county = count + 1
wx = subwrd(wcmd,countx)
wy = subwrd(wcmd,county)
if ((wx = '') | (wy = ''))
break
endif
* Convert world coordinates to screen coordinates if necessary
if (newgrads)
sx = wx
sy = wy
else
'q w2xy 'wx' 'wy
sx = subwrd(result,3)
sy = subwrd(result,6)
endif
* Append the coordinates to the draw command
cmd = cmd%sx' 'sy' '
count = count + 2
endwhile
* Draw the polygon
'set line 'color
cmd
endif
endwhile
* Draw the continental outline in the specified color
'set map 'outline
'draw map'
'set map auto'
* Draw a new square frame around the plot
* If you have 'set frame off' or 'set frame circle' before running basemap,
* you may want to comment out the next 2 lines
'set line 1 1 6'
'draw rec 'x1' 'y1' 'x2' 'y2
* Close the polygon file
rc = close(file)
return
* THE END *
这样的话，图可以画出来，但是屏蔽不了陆地区域。
非常感谢您，不厌其烦地回复，鞠躬。 |
|
[复制链接]
窥视卡
雷达卡
楼主
显身卡