* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* shp_demo2.gs
*
* Demonstrates a way to use the shapefile interface in GrADS by
* creating a new shapefile based on a precipitation forecast 
* and then drawing the shapefile with the sea level pressure contours
*
* Written by Jennifer Adams, June 2010
*
function main()

'reinit'
'set display color white'; 'clear'

* N.B. Don't forget to change the /path/dataset name for your own data
'use /data/grib/gfs/gfs'

* Create the shapefile with precip > 0.8
'set t 3'
'set lat 22 74'
'set lon 115 253'
'myp=maskout(p,p-0.8)'
'set gxout shp'
'set shp -pt shppt'
'd myp'
num = subwrd(result,1)

* Draw the SLP contours
'set mproj nps'
'set mpvals 140 220 34 65.6'
'set gxout contour'
'set ylint 10'
'set xlint 30'
'set ccolor rainbow'
*'set clab masked'
'set mpdset mres'
'set cthick 6'
'set clevs 988 996 1004 1012 1020 1028'
'set ccols 14  11  13   10   12   2   '
'd smth9(smth9(slp/100))'
'draw map'
'set cthick 6'
'set clab off'
'set clevs 980 984 992 1000 1008 1016 1024 1032'
'set ccols 9   9   4   5    3    7    8    6   '
'd smth9(smth9(slp/100))'

* overlay stippling from shapefile with a few gray shades
'set rgb 20 180 180 180'
'set rgb 25 130 130 130'
'set rgb 30  80  80  80'
'set rgb 35  30  30  30'
ccols = '20 25 30 35'
clevs = '  1  2  5'
c=1; while (c<=18); levs.c = subwrd(clevs,c); c=c+1; endwhile
c=1; while (c<=19); cols.c = subwrd(ccols,c); c=c+1; endwhile

'q dbf shppt'
data = result
n=0
while (n<num)
  db = sublin(data,n+2)
  len = math_strlen(db)
  c=1
* get position of final comma
  while(c<=len)
    char = substr(db,c,1)
    if (char=','); pos=c; endif
    c=c+1
  endwhile
* parse value from end of comma-delimited list of attributes
  val = substr(db,pos+1,len-pos)
* the fill color will depend on the data value
    if (val <= levs.1               ) ; 'set line 'cols.1; endif 
    if (val <= levs.2 & val > levs.1) ; 'set line 'cols.2; endif
    if (val <= levs.3 & val > levs.2) ; 'set line 'cols.3; endif
    if (                val > levs.3) ; 'set line 'cols.4; endif
* draw the point
  'set shpopts -1 5 0.025'
  'draw shp shppt 'n
  n=n+1
endwhile

'draw title Sea Level Pressure Contours with Stippled Precipitation'