Re: [S] Points in or out of polygon

Jeff Longmate (jlongmate@smtplink.coh.org)
Wed, 04 Mar 98 15:29:25 -0800


"Gebhart Kristi" <Gebhart@CIRA.colostate.edu> writes:

> I have a need to determine whether or not points are in or out of
> arbitrarily shaped and located polygons. Has anyone written an S-PLUS
> function that would do something like, given the x-y coordinates of the
> vertices of a polygon, determine whether the x-y coordinates of another
> point are inside or outside of it?

I couldn't resist a try at something a little more general.
The idea is that if you go in any direction, say due north,
and count boundary crossings, an odd count means your in, and even
means you are out, assuming the region is defined by closed curves
(here connect-the-dots loops). Concave is ok. Islands and holes
are ok. Assuming the data are in a list of n_i by 2 matrices, each
with the first point replicated at the end (so a call to lines() draws
a closed loop), then this seems to work:

is.north <- function(b,x) approx(b[c(1,3)],b[c(2,4)],x[1])$y>x[2]
crossings <- function(b,x){
pairs <- cbind(b[-nrow(b),],b[-1,])
dimnames(pairs)[[2]] <- c("x1","y1","x2","y2")
maybe <- pairs[ (pairs[,1] < x[1] & pairs[,3] >= x[1])|
(pairs[,1] >= x[1] & pairs[,3] < x[1]),,drop=F]
if (length(maybe) == 0) result <- 0
else result <- sum(apply(maybe,1,is.north,x))
result
}
inside <- function(x,bound)
as.logical(sum(sapply(bound,crossings,x))%%2)

First get pairs of points. Then collect the pairs (in maybe) whose x
coordinates straddle the test point (x). Then use approx to test if the
candidate segments are to the north. Here's an example of use:

# Draw boundaries on a 100 by 100 field via mouse:
plot(c(0,100),c(0,100),type="n")
get.points <- function(n){
points <- sapply(locator(n),round) # enter via mouse
rbind(points,points[1,]) # Copy first point to end.
}
points1 <- get.points(20)
lines(points1)
points2 <- get.points(10)
lines(points2)
points3 <- get.points(10)
lines(points3)
points <- list(points1,points2,points3)

# Now we have an interesting boundary. Test a point:
inside(unlist(locator(1)),points)

SEEMS to work.

Jeff Longmate
jlongmat@smtplink.coh.org

-----------------------------------------------------------------------
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