[ferret_users] Creating and using an irregular mask, as in the recent question regarding filling data over China

[mick spillane <mick.spillane@xxxxxxxx>]

Hello All,
   A recent question dealt with how to fill (or otherwise analyze) data 
from within an irregular region.  It arose in the context of China but 
has cropped up several times in the archives and the answer of course is 
to define a mask variable to isolate the region in question.  This is 
easier said than done, when the boundary of the region is convoluted, 
but the attached script may be of some help.  It is based on deciding 
which of the deciding which cells of the underlying grid of the field 
being analyzed lie within a polygon defining the region 
(InsidePolygon.jnl).  Its application is aided by using a Ferret-based 
polygon generator (polydef.jnl) using mouse clicks that works on linux 
machines.

   As an illustration consider the question -- "What is the mean 
elevation of New Zealand (based on the etopo20 topography file)?" (This 
was chosen because it can be verified by a simple calculation not 
involving polygons.)

!  average elevation of new zealand using polygon definition
use etopo20 ; let topo=rose[d=1]
shade/x=160:180/y=-50:-30/... rose
go polydef                                        ! this allows user to 
define the polygon by a series of mouse-clicks
                                                        ! (see the 
little x's in the attached graphic) Because of the two islands
                                                        ! a land bridge 
was used to have a single polygon circumscribe both.
go InsidePolygon                             ! this identifies etopo20 
cells that lie within the polygon through a mask
                                                        ! variable 
"inside" that is 1 for inside cells and missing elsewhere
The   ! illustrates the cells selected
let masked_elevation=inside*topo
list/nohead masked_elevation[x=@ave,y=@ave]  ! which in this case gives 
317.8 meters
list/nohead inside[x=@ngd,y=@ngd]                    ! for the 267 cells 
inside the polygon

! check -- "simple method" using the fact that new zealand is surrounded 
by water
let land=if(topo gt 0)then topo
list/nohead/x=160:180/y=-50:-30 land[x=@ave,y=@ave] ! this gives 
320.5meters based on the 265 land cells.

! NOTE: elimination of the 2 "land bridge" cells in the the polygon 
solution ...
let onland=if(topo ge 0)then inside
let land_elevation=onland*topo
list/nohead land_elevation[x=@ave,y=@ave]  ! ... gives the same answer 
320.5m but this is beside the point.

   For the real problem of say the average elevation of a land-locked 
country (where the "simple method" is not available) one would procede 
as follows:
1) obtain a file of the political boundaries -- if the country's border 
is described by a continuous series of lon,lat values that might be the 
polygon to use. however take care that it is not made up of 
non-continuous segments, or repeated or too closely spaced points that 
might foul up the InsidePolygon calculation.  Note that 
InsidePolygon.jnl assumes that "vertices.xy" is the filename of the 
polygon and that the polygon is closed by repeating the first point at 
the end.
2) if the political boundary file is not to be used directly as a 
polygon file do an overlay plot showing (in this case) the topography 
and the political boundary:
   shade/x=.../y=... topo
   file/form=free/var=xx,yy PoliticalBoundary.file
   plot/o/vs/nolab/line=2 xx,yy
3) now use the "polydef.jnl" tool to construct a closed polygon suited 
to the topography file; close the polygon by clicking on the initial 
point before the final click to the left of the plot area that signifies 
the polygon has been completed.  Then, as before,
   go InsidePolygon
   shade/o/nolab/pal=red/pat=dark_vertical inside     ! to verify that 
the polygon matches expectation
   let masked_elevation=inside*topo                         ! or 
whatever masked variable is desired
   list/nohead masked_elevation[x=@ave,y=@ave]  ! obtain the answer

The second attached image illustrates the use of InsidePolygon.jnl to 
get the average elevation of Switzerland using the smith-sandwell data 
(the colored background) and the political boundary provided by Ferret 
(go land_detail " " " " 4).  The x's are the points used in digitizing 
the border and the red hatching the polygon that results.  The average 
elevation of the hatched area is 1278 meters.  I've had a bit of trouble 
in finding online a value to check this against.  The min and max of 
masked_elevation are 193m and 3797m respectively.  The minimum matches 
exactly the online value (Lake Maggiore), but the max is somewhat less 
that the 4634m of Mount Monte Rosa - (smith sandwell too coarse).  
However the area returned (inside[x=@din,y=@din] = 4.124E+10 m^2 agrees 
well with the online value of 41285 km². So all-in-all I think the 
validity and reative ease of use of the InsidePolygon method has been 
demonstrated.

Hope the scripts prove useful; and the documentation is adequate.  I  
would put the scripts in ~/ferret subdirectory or wherever your 
my_ferret_paths points to. Unfortunately at the moment for me, "polydef" 
only works on Linux, not on Mac Ferret.  If someone sees why this is so 
and has a fix for Mac, I'd be delighted to learn.

Good luck,
   Mick

! InsidePolygon : Tests if a point (X0,Y0) is inside a closed polygon defined
!                 by a set of vertices (VX,VY),K=1,NV whose last point duplicates
!                 its first. It is based on the number of intersections
!                 between the NV-1 edges of the polygon and a line from (X0,Y0)
!                 to (X0,YTOP) where YTOP is the upper limit of the domain (90
!                 when X,Y are lon/lat coords).
!                 The point is INSIDE if the #intersections is ODD.

! Written 26-May-2009 by Mick.Spillane@xxxxxxxx

let YTOP=90 ; let NEDGE=VX[k=@ngd]-1

! An edge is a candidate if X0 lies between VX and VX[k=@shf] ...
let XWORKS=if((X0-VX)*(X0-VX[k=@shf]) lt 0)then 1

! ... but the Y-value of the edge, at X=X0, must also be between Y0 and YTOP
let YPRIME=VY+XWORKS*(VY[k=@shf]-VY)*(X0-VX)/(VX[k=@shf]-VX)
let ITCUTS=if((YPRIME-Y0)*(YPRIME-YTOP) lt 0)then 1
let INSIDE=if(mod(ITCUTS[K=1:`NEDGE`@ngd],2) eq 1)then 1

! polydef : use mouse click to define polygon vertices

can mode verify
let done=0 ; sp rm -f vertices.xy

say "****************************************************"
say "*                                                  *"
say "* Add polygon vertices by mouse clicks.  Terminate *"
say "*     by clicking to the left of the plot area.    *"
say "*                                                  *"
say "****************************************************"

! add new vertices to the file vertices.xy
repeat/range=1:1000 go add_vertex

! then read in the resulting file 
sp get_vertices ; go get_vertices

set mode verify

! add a new vertex to the vertex file while done=0

if `done eq 0` then 
  where
  if `($XMOUSE) gt ($XAXIS_MIN)` then
    list/nohead/app/form=(2f12.6)/file=vertices.xy ($XMOUSE),($YMOUSE)
! if visible marks at the vertices are not desired drop the next line
    plot/o/nolab/vs/sym=1 ($XMOUSE),($YMOUSE)
  else
    let done=1
  endif
endif