*-------------------------------------------------------------------------------
*  GrADS script to mask out unwanted areas from the plot (Neatly !).
*
*  Rupa Kumar Kolli, January 11, 1998.
*  Please report bugs/suggestions to kolli@tropmet.ernet.in.
*-------------------------------------------------------------------------------
*  Reads polygon definitions in world coordinates from a text file
*  and draws filled polygons in black colour.
*
*  Note:
*
*    -  Should be called ONLY after displaying a plot in GrADS window.
*       use 'set mpdraw off' for the plot; map is drawn in this script.
*
*    -  This script erases grid lines within the polygons; so if you
*       want full gridlines, give xlint and ylint as arguments for
*       this function.  Further, use 'set grid off' before calling this.
*
*    -  Edit the line 'myfile=' below to change your data file containing
*       polygon definitions; see data file 'indmask.txt' for format.
*       Also, edit the line 'set mpdset rupres' to specify the GrADS
*       map data set on which the polygon definitions are based.
*
*    -  For better results, extend your polygons to outside the frame
*       of interest; clipping to frame is done in the script.
*
*    -  If xlint/ylint are not defined, no grid lines are drawn in the
*       erased areas.  If xstart/ystart are not defined, it is assumed that
*       the tick marks start at the axis minima.
*
*  Syntax (all gridline arguments optional):
*
*    grmask maskfile <mapfile> <xlint> <ylint> <xstart> <ystart>
*-------------------------------------------------------------------------------
function grmask(args)
*
*  This is the input file containing Polyfill Definitions in world coords
*
myfile=subwrd(args,1)
if(myfile='');say 'FATAL ERROR: Please specify maskfile !';return;endif
mapfile=subwrd(args,2)
if(mapfile='');mapfile='rupres';endif
*
*  This is the map data set on which the above Polyfill Definitions are based
*  You can create your own GrADS map data sets from ASCII data
*  Contact me for a FORTRAN program that does this.
*
'set mpdset 'mapfile
*
*
xlint=subwrd(args,3)
ylint=subwrd(args,4)
xstart=subwrd(args,5)
ystart=subwrd(args,6)
*
*  Get the physical dimensions of the plot frame
*
'query gxinfo'
rec=sublin(result,3)
xlo=subwrd(rec,4)
xhi=subwrd(rec,6)
rec=sublin(result,4)
ylo=subwrd(rec,4)
yhi=subwrd(rec,6)
*
*  Set fill color to 0 (black) and clip to the current plot frame
*
'set line 0'
'set clip 'xlo' 'xhi' 'ylo' 'yhi
*
*  Get number of Polyfill definition sets from file
*
ret=read(myfile)
if(sublin(ret,1) != 0)
  ret=read('d:/proj/da/'%myfile)
  if(sublin(ret,1) != 0)
    say 'Error opening Polyfill Define File !'
    say 'The file causing error is: '%myfile
    say 'Masking Aborted.'
    return
  endif
  myfile='d:/proj/da/'%myfile
endif
nsets=sublin(ret,2)
*
*  Loop over the Polyfill definition sets
*
isets=1
while(isets<=nsets)
*
*  Get number of lines for this set
*
  ret=read(myfile)
  rec=sublin(ret,2)
  num=subwrd(rec,2)
*
*  Read each line and concatenate all lines into one string
*
  i=1
  s=''
  while(i<=num)
    ret=read(myfile)
    rec=sublin(ret,2)
    s=s%rec
    i=i+1
  endwhile
*
*  Convert all pairs of coordinates in the concatenated string from
*  world coordinates to XY coordinates and concatenate the
*  XY coordinates into a new string
*
  iw=1
  m=''
  wx=subwrd(s,iw)
  while(wx !='')
    wy=subwrd(s,iw+1)
    'q w2xy 'wx' 'wy
    x=subwrd(result,3)
    y=subwrd(result,6)
    m=m%' '%x%' '%y
    iw=iw+2
    wx=subwrd(s,iw)
  endwhile
*
*  draw a filled polygon
*
  'draw polyf'm
  isets=isets+1
endwhile
*
*  End loop on sets
*  Draw the map and frame
*
'set mpdraw on'
'set mpt * -1'
'set mpt 3 15 1 2'
'set mpt 2 14 1 4'
'set mpt 1 4 1 8'
'draw map'
'set line 1 1 3'
'draw rec 'xlo' 'ylo' 'xhi' 'yhi
*
*  Get world coordinate limits of plot frame to draw grid lines
*
if(xlint = ''); return; endif
'query dims'
rec=sublin(result,2)
wxlo=subwrd(rec,6)
wxhi=subwrd(rec,8)
rec=sublin(result,3)
wylo=subwrd(rec,6)
wyhi=subwrd(rec,8)
*
*  Set line style to dotted line (5)
*
'set line 1 5 1'
*
*  Draw vertical grid lines from first argument
*
if(xstart = '' | wxlo=xstart)
  wgx=wxlo+xlint
else
  wgx=xstart
endif
while (wgx < wxhi)
  'q  w2xy 'wgx' 'wylo
  gx=subwrd(result,3)
  'draw line 'gx' 'ylo' 'gx' 'yhi
  wgx=wgx+xlint
endwhile
*
*  Draw horizontal grid lines from second argument
*
if (ylint != '')
  if(ystart = '' | wylo=xstart)
    wgy=wylo+ylint
  else
    wgy=ystart
  endif
  while (wgy < wyhi)
    'q  w2xy 'wxlo' 'wgy
    gy=subwrd(result,6)
    'draw line 'xlo' 'gy' 'xhi' 'gy
    wgy=wgy+ylint
  endwhile
endif
*
*  That's it, return to caller
*
return