[S] Re: Points in or out of polygon

Mark Bravington FSMG CEFAS (M.V.BRAVINGTON@cefas.co.uk)
Thu, 5 Mar 1998 09:20:40 +0000


[ Kristi Gebhart asked about testing whether points are inside or outside a given polygon:]

Kristi,

Here's a pair of functions I wrote a while back to do the job. I presume the code is quite
efficient, because it is incomprehensibly dense. It uses the "count-the-crossings" approach. NB that the result is not well-defined for points on the edge; it may come out either T or F.

The function 'in.poly' is the main one. The parameter "xy" should be a
data.frame or 2-column matrix of x and y coordinates. "poly" can either be:

(1) a n-vertex polygon in the same format with n rows (last row does not need to
equal first row as this is done automatically).

(2) the result of a previous call to 'poly.for.testing' with the original
polygon as input.

In other words, the following are equivalent:

> points_ data.frame( x=c(0,1), y=c( 0, 1.5))
> polypoints_ data.frame( x=c( 0.5, 2, 2), y=c( 0, 0, 6))
> in.poly( points, polypoints)
[1] F T
> temp.poly_ poly.for.testing( polypoints)
> in.poly( points, polypoints)
[1] F T

The result is, as you see, a logical vector with one element per point.

If you are repeatedly testing lots of points for inclusion within the same
polygon, and you're doing it within a loop rather than with one call to
"in.poly", it is much more efficient to call 'poly.for.testing' first.

Here's the code-- good luck

in.poly_ function( xy, poly)
{
if( ncol( poly)==2)
poly_ poly.for.testing( poly)
n_ nrow( xy); np_ nrow( poly); nnp_ rep( n,np)
check1_ xor( xy[,1]>=rep( poly[,1], nnp), xy[,1]>rep( poly[,3], nnp))
check2_ rep( poly[,2], nnp) + rep( poly[,4], nnp) * xy[,1] > xy[,2]
as.vector( matrix( check1 & check2, n, np) %*% rep.int( 1, np) %% 2 > 0)
}
poly.for.testing_ function( xy)
{
if( is.list( xy))
xy_ as.data.frame(xy)
if( !is.matrix( xy) || ncol( xy)!=2 || !is.numeric( xy[,1]) ||
!is.numeric( xy[,2]))
stop( "xy must by nX2 numeric matrix or data.frame")

# Columns of poly are start.x, a, end.x, b
xy_ as.matrix( xy)
poly_ matrix( c( xy, xy[ c( 2:nrow( xy), 1), ]), ncol=4)
poly_ poly[ poly[,1] != poly[,3], ]
poly[,4]_ (poly[,4]-poly[,2])/(poly[,3]-poly[,1])
poly[,2]_ poly[,2] - poly[,1]*poly[,4]
poly
}

Mark Bravington

-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news