This is a great contribution! There are codes that do point-in-polygon calculations and I believe they do similar calculations, counting how many times a ray crosses the closed. Using Ferret to define the closed curve around a country is something I hadn't thought of. Notice in his figures how he uses the outline of New Zealand to make a polygon at the same resolution as the grid the data is on.
You could use the scripts land or land_detail to make a plot showing boundaries for nations or states or even lakes, and then run polydef.jnl, clicking on those curves to create boundaries for them.
yes? go/help land
yes? go/help land_detail
Mick and I were just talking about this and trying some examples and found a few additions to make. We changed polydef.jnl so it will automatically duplicate the first point at the end to make sure the polygon is closed. I'm attaching a tar file with the scripts:
polydef.jnl ! The Ferret script to use mouse-clicks to define a polygon. Click
! to the left of the x axis to end.
add_vertex.jnl ! called by polydef. This is changed to close the polygon.
get_vertices ! A Unix shell script that writes a script to define a
! grid axis and opens the vertex file. Make this script
! executable; perhaps install it in your bin/ directory.
InsidePolygon.jnl ! Defines the varible INSIDE which is 1 inside the polygon,
! 0 outside it. An example script is included at its start.
Here are the examples from the start of InsidePolygon.jnl
!Examples:
!1) Using polydef to define the vertices vx and vy
yes? use etopo60
yes? shade rose
yes? go polydef ! Click to define a set of points. Or otherwise define vx and vy
yes? let x0 = x[gx=rose,d=1]
yes? let y0 = y[gy=rose,d=1]
yes? go InsidePolygon
yes? shade/over/palette=gray inside
!2) Read vertices vx and vy from an existing file
! say this file contains xloc, yloc defining some boundary.
yes? use etopo60
yes? shade rose
yes? let x0 = x[gx=rose,d=1]
yes? let y0 = y[gy=rose,d=1]
yes? use my_vertices.nc
yes? let vx = zsequence(xloc)
yes? let vy = zsequence(yloc)
yes? go InsidePolygon
yes? shade/over/palette=gray inside
mick spillane wrote:
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 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
