**To**:**Ferret users <oar.pmel.ferret_users@xxxxxxxx>****Subject**:**[ferret_users] Creating and using an irregular mask, as in the recent question regarding filling data over China****From**:**mick spillane <mick.spillane@xxxxxxxx>**- Date: Tue, 26 May 2009 16:56:26 -0700
- Sender: owner-ferret_users@xxxxxxxx
- User-agent: Thunderbird 2.0.0.21 (Macintosh/20090302)

Hello All,

! average elevation of new zealand using polygon definition use etopo20 ; let topo=rose[d=1] shade/x=160:180/y=-50:-30/... rose

The ! illustrates the cells selected let masked_elevation=inside*topo

let land=if(topo gt 0)then topo

let onland=if(topo ge 0)then inside let land_elevation=onland*topo

shade/x=.../y=... topo file/form=free/var=xx,yy PoliticalBoundary.file plot/o/vs/nolab/line=2 xx,yy

go InsidePolygon

list/nohead masked_elevation[x=@ave,y=@ave] ! obtain the answer

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

**Follow-Ups**:

- Previous by thread:
**[ferret_users] ***** Announcing the official release of Ferret version 6.2 ******* - Next by thread:
**Re: [ferret_users] Creating and using an irregular mask, as in the recent question regarding filling data over China**

Contact Us

Dept of Commerce / NOAA / OAR / PMEL / TMAP

Privacy Policy | Disclaimer | Accessibility Statement