[gradsusr] maskout problem

[rocio vazquez]

You can use grmask.gs script. I used it to mask my country. You have to
create a poligon from your country. Here i send you the script, i hope it
work.


GrADS script to mask out unwanted areas from the plot (Neatly !).
*
*  Rupa Kumar Kolli, January 11, 1998.
*  Please report bugs/suggestions to kolli at 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('/usr/local/grads/dat/'%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='/usr/local/grads/dat/'%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 0 1 1 4'
'set mpt 1 1 1 4'
'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
*

2015-01-20 2:31 GMT-03:00 Mandal Raju <raju.cat at tropmet.res.in>:

>
>  Hi,
>
>  I have a plot of temperature anomaly over the region: lat 0 to 40N
>                                                        lon 60E to 100E
>
>  But I want to show the values only over Indian sub-continent. So, can
> anybody help me how to maskout those outside region from the actual plot?
>
>
>
>  Thanks & Regards....
>
>    Raju Mandal
>    India
>    Office: 020 2590 4513
>
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20150120/0b4d1a5a/attachment.html