Finds points in a polygon

Description

Given a polygon and a set of points this function returns the subset of points that are within the polygon.

Usage

1

Arguments

poly

an nx2 matrix of polygon vertices. Matrix columns correspond to vertices' x and y coordinates, respectively.

points

an mx2 matrix of points. Matrix columns correspond to points' x and y coordinates, respectively.

...

currently no additional arguments.

Details

It is assumed that the polygon is to be closed by joining the last vertex to the first vertex.

Value

If points are found with the polygon, then a vector is returned with elements corresponding to the row indices of points, otherwise NA is returned.

Author(s)

Andrew O. Finley finleya@msu.edu,
Sudipto Banerjee sudiptob@biostat.umn.edu,

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
## Not run: 
##Example 1
points <- cbind(runif(1000, 0, 10),runif(1000, 0, 10))

poly <- cbind(c(1:9,8:1), c(1,2*(5:3),2,-1,17,9,8,2:9))

point.indx <- pointsInPoly(poly, points)

plot(points, pch=19, cex=0.5, xlab="x", ylab="y", col="red")
points(points[point.indx,], pch=19, cex=0.5, col="blue")
polygon(poly)

##Example 2
##a function to partition the domain
tiles <- function(points, x.cnt, y.cnt, tol = 1.0e-10){

  x.min <- min(points[,1])-tol
  x.max <- max(points[,1])+tol
  y.min <- min(points[,2])-tol
  y.max <- max(points[,2])+tol
  
  x.cnt <- x.cnt+1
  y.cnt <- y.cnt+1
  
  x <- seq(x.min, x.max, length.out=x.cnt)
  y <- seq(y.min, y.max, length.out=y.cnt)

  tile.list <- vector("list", (length(y)-1)*(length(x)-1))

  l <- 1
  for(i in 1:(length(y)-1)){
    for(j in 1:(length(x)-1)){
      tile.list[[l]] <- rbind(c(x[j], y[i]),
                              c(x[j+1], y[i]),
                              c(x[j+1], y[i+1]),
                              c(x[j], y[i+1]))
       l <- l+1
    }

  }

  tile.list
}

n <- 1000
points <- cbind(runif(n, 0, 10), runif(n, 0, 10))

grd <- tiles(points, x.cnt=10, y.cnt=10)

plot(points, pch=19, cex=0.5, xlab="x", ylab="y")

sum.points <- 0
for(i in 1:length(grd)){
  polygon(grd[[i]], border="red")
  
  point.indx <- pointsInPoly(grd[[i]], points)
  
  if(!is.na(point.indx[1])){
    sum.points <- length(point.indx)+sum.points
    
    text(mean(grd[[i]][,1]), mean(grd[[i]][,2]), length(point.indx), col="red")
  }
}
sum.points


## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.