[Thread Prev][Thread Next][Index]

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

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,

GIF image

GIF image

! 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 
  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)
    let done=1

[Thread Prev][Thread Next][Index]

Contact Us
Dept of Commerce / NOAA / OAR / PMEL / TMAP

Privacy Policy | Disclaimer | Accessibility Statement