Hi all,
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
|