Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a SpatialPointsDataFrame with some additional data. I would like to extract those points inside a polygon and at the same time, preserve SPDF object and its corresponding data.

So far I've had little luck and resorted to matching and merging through a common ID, but this works only becaues I have gridded data with individual IDS.

Here's a quick example, I'm looking for points inside the red square.

library(sp)
set.seed(357)
pts <- data.frame(x = rnorm(100), y = rnorm(100), var1 = runif(100), var2 = sample(letters, 100, replace = TRUE))
coordinates(pts) <- ~ x + y
class(pts)
plot(pts)
axis(1); axis(2)

ply <- matrix(c(-1,-1, 1,-1, 1,1, -1,1, -1,-1), ncol = 2, byrow = TRUE)
ply <- SpatialPolygons(list(Polygons(list(Polygon(ply)), ID = 1)))
ply <- SpatialPolygonsDataFrame(Sr = ply, data = data.frame(polyvar = 357))
plot(ply, add = TRUE, border = "red")

The most obvious approach would be to use over, but this returns the data from the polygon.

> over(pts, ply)
    polyvar
1        NA
2       357
3       357
4        NA
5       357
6       357
share|improve this question
    
Thanks for providing a reproducible example. Always helps when trying to understand a problem! –  flowla Jun 18 '13 at 16:00

4 Answers 4

up vote 4 down vote accepted

From the sp:over help:

 x = "SpatialPoints", y = "SpatialPolygons" returns a numeric
      vector of length equal to the number of points; the number is
      the index (number) of the polygon of ‘y’ in which a point
      falls; NA denotes the point does not fall in a polygon; if a
      point falls in multiple polygons, the last polygon is
      recorded.

So if you convert your SpatialPolygonsDataFrame to SpatialPolygons you get back a vector of indexes and you can subset your points on NA:

> over(pts,as(ply,"SpatialPolygons"))
  [1] NA  1  1 NA  1  1 NA NA  1  1  1 NA NA  1  1  1  1  1 NA NA NA  1 NA  1 NA
 [26]  1  1  1 NA NA NA NA NA  1  1 NA NA NA  1  1  1 NA  1  1  1 NA NA NA  1  1
 [51]  1 NA NA NA  1 NA  1 NA  1 NA NA  1 NA  1  1 NA  1  1 NA  1 NA  1  1  1  1
 [76]  1  1  1  1  1 NA NA NA  1 NA  1 NA NA NA NA  1  1 NA  1 NA NA  1  1  1 NA

> nrow(pts)
[1] 100
> pts = pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),]
> nrow(pts)
[1] 54
> head(pts@data)
         var1 var2
2  0.04001092    v
3  0.58108350    v
5  0.85682609    q
6  0.13683264    y
9  0.13968804    m
10 0.97144627    o
> 

For the doubters, here's evidence that the conversion overhead is not a problem:

Some functions - first Jeffrey Evans' method, then my original, then my hacky conversion, then a version based on gIntersects based on Josh O'Brien's answer:

evans <- function(pts,ply){
  prid <- over(pts,ply)
  ptid <- na.omit(prid) 
  pt.poly <- pts[as.numeric(as.character(row.names(ptid))),]
  return(pt.poly)
}

rowlings <- function(pts,ply){
  return(pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),])
}

rowlings2 <- function(pts,ply){
  class(ply) <- "SpatialPolygons"
  return(pts[!is.na(over(pts,ply)),])
}

obrien <- function(pts,ply){
pts[apply(gIntersects(columbus,pts,byid=TRUE),1,sum)==1,]
}

Now for a real-world example, I've scattered some random points over the columbus data set:

require(spdep)
example(columbus)
pts=data.frame(
    x=runif(100,5,12),
    y=runif(100,10,15),
    z=sample(letters,100,TRUE))
coordinates(pts)=~x+y

Looks good

plot(columbus)
points(pts)

Check the functions are doing the same thing:

> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] TRUE

And run 500 times for benchmarking:

> system.time({for(i in 1:500){evans(pts,columbus)}})
   user  system elapsed 
  7.661   0.600   8.474 
> system.time({for(i in 1:500){rowlings(pts,columbus)}})
   user  system elapsed 
  6.528   0.284   6.933 
> system.time({for(i in 1:500){rowlings2(pts,columbus)}})
   user  system elapsed 
  5.952   0.600   7.222 
> system.time({for(i in 1:500){obrien(pts,columbus)}})
  user  system elapsed 
  4.752   0.004   4.781 

As per my intuition, its not a great overhead, in fact it might be less of an overhead than converting all the row indexes to character and back, or running na.omit to get missing values. Which incidentally leads to another failure mode of the evans function...

If a row of the polygons data frame is all NA (which is perfectly valid), then the overlay with SpatialPolygonsDataFrame for points in that polygon will produce an output data frame with all NAs, which evans() will then drop:

> columbus@data[1,]=rep(NA,20)
> columbus@data[5,]=rep(NA,20)
> columbus@data[17,]=rep(NA,20)
> columbus@data[15,]=rep(NA,20)
> set.seed(123)
> pts=data.frame(x=runif(100,5,12),y=runif(100,10,15),z=sample(letters,100,TRUE))
> coordinates(pts)=~x+y
> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] FALSE
> dim(evans(pts,columbus))
[1] 27  1
> dim(rowlings(pts,columbus))
[1] 28  1
> 

BUT gIntersects is faster, even with having to sweep the matrix to check intersections in R rather than in C code. I suspect its the prepared geometry skills of GEOS, creating spatial indexes - yeah, with prepared=FALSE it takes a bit longer, about 5.5 seconds.

I'm surprised there isn't a function to either straight return the indices or the points. When I wrote splancs 20 years ago the point-in-polygon functions had both...

share|improve this answer
    
Great, this also works for multiple polygons (I've added an example to play with into Joshua's answer). –  Roman Luštrik Jun 19 '13 at 10:15
    
With large polygon datasets coercion into a SpatialPolygons object is a lot of overhead and not necessary. Applying "over" to a SpatialPolygonsDataFrame returns the row index which can be used to subset the points. See my example below. –  Jeffrey Evans Jun 19 '13 at 14:58
    
A lot of overhead? Its essentially just taking the @polygons slot from the SpatialPolygonsDataFrame object. You can even 'fake' it by reassigning the class of a SpatialPolygonsDataFrame to be "SpatialPolygons" (although this is hacky and not recommended). Anything that is going to use the geometry is going to have to get that slot at some stage anyway, so relatively speaking its no overhead at all. Its insignificant anyway in any real-world application where you are then going to be doing a load of point-polygon tests. –  Spacedman Jun 19 '13 at 23:20
    
There is more than speed considerations in accounting for overhead. In creating a new object in the R namespace you are using necessary RAM. Where this is not a problem in a small datasets it will effect performance with large data. R does exhibit a linear performance die off. As the data gets larger performance takes a ding. If you do not need to create an additional object why would you? –  Jeffrey Evans Jun 20 '13 at 16:59
    
Good point about NA's in the polygon dataframe. However, you can utilize match or which to return an "honest" vector representing row indexes without worrying about NA's in the polygon dataframe. Besides, as we are arguing over this everybody is ignoring that that most efficient solution proposed was using rgeos and gIntersects. –  Jeffrey Evans Jun 20 '13 at 17:07

Here's a possible approach using the rgeos package. Basically, it makes use of the gIntersection function that allows you to intersect two sp objects. By extracting the IDs of those points that lie within the polygon, you are afterwards able to subset your original SpatialPointsDataFrame, keeping all the corresponding data. The code is almost self-explaining, but if there are any questions, please feel free to ask!

# Required package
library(rgeos)

# Intersect polygons and points, keeping point IDs
pts.intersect <- gIntersection(ply, pts, byid = TRUE)

# Extract point IDs from intersected data
pts.intersect.strsplit <- strsplit(dimnames(pts.intersect@coords)[[1]], " ")
pts.intersect.id <- as.numeric(sapply(pts.intersect.strsplit, "[[", 2))

# Subset original SpatialPointsDataFrame by extracted point IDs
pts.extract <- pts[pts.intersect.id, ]

head(coordinates(pts.extract))
              x          y
[1,] -0.5832050 -0.8777367
[2,]  0.3947471  0.7020481
[3,]  0.7667997 -0.9465043
[4,]  0.3174604  0.6416281
[5,] -0.4690151  0.4413502
[6,]  0.4765213  0.6068021

head(pts.extract)
         var1 var2
2  0.04001092    v
3  0.58108350    v
5  0.85682609    q
6  0.13683264    y
9  0.13968804    m
10 0.97144627    o
share|improve this answer
1  
Should tmp be pts.intersect? Also, parsing the returned dimnames like that is depending on undocumented behaviour. –  Spacedman Jun 20 '13 at 17:25
    
@Spacedman, you're right about tmp, forgot to remove it when finishing the code. Also, you're right about parsing the dimnames. This was kind of a quick solution to provide the questioner with a fast answer, and there surely are better (and more universal) approaches, for example yours :-) –  flowla Jun 21 '13 at 7:57

Is this what you are after?

One note, on edit: The wrapping call to apply() is needed to make this work with arbitrary SpatialPolygons objects, possibly containing more than one polygon feature. Thanks to @Spacedman for prodding me to demonstrate how to apply this to the more general case.

library(rgeos)
pp <- pts[apply(gIntersects(pts, ply, byid=TRUE), 2, any),]


## Confirm that it works
pp[1:5,]
#              coordinates       var1 var2
# 2 (-0.583205, -0.877737) 0.04001092    v
# 3   (0.394747, 0.702048) 0.58108350    v
# 5    (0.7668, -0.946504) 0.85682609    q
# 6    (0.31746, 0.641628) 0.13683264    y
# 9   (-0.469015, 0.44135) 0.13968804    m

plot(pts)
plot(ply, border="red", add=TRUE)
plot(pp, col="red", add=TRUE)
share|improve this answer
    
Fails horribly if ply has more than one feature, because gIntersects returns a matrix with one row for each feature. You probably can probably sweep the rows for a TRUE value. –  Spacedman Jun 20 '13 at 17:20
    
@Spacedman -- Bingo. Need to do apply(gIntersects(pts, ply, byid=TRUE), 2, any). In fact, I'll go ahead and switch the answer to that, since it encompasses the case of a single polygon as well. –  Josh O'Brien Jun 20 '13 at 17:27
    
Ah, any. That might be marginally faster than the version I just benchmarked. –  Spacedman Jun 20 '13 at 17:38
    
@Spacedman -- From my quick tests, it looks like obrien and rowlings2 run neck and neck, with obrien maybe 2% faster. –  Josh O'Brien Jun 20 '13 at 17:43

You were on the right track with over. The rownames of the returned object correspond to the row index of the points. You can implement your exact approach with just a few addition lines of code.

library(sp)
set.seed(357)

pts <- data.frame(x=rnorm(100), y=rnorm(100), var1=runif(100), 
                  var2=sample(letters, 100, replace=TRUE))
  coordinates(pts) <- ~ x + y

ply <- matrix(c(-1,-1, 1,-1, 1,1, -1,1, -1,-1), ncol=2, byrow=TRUE)
  ply <- SpatialPolygons(list(Polygons(list(Polygon(ply)), ID=1)))
    ply <- SpatialPolygonsDataFrame(Sr=ply, data=data.frame(polyvar=357))

# Subset points intersecting polygon
prid <- over(pts,ply)
  ptid <- na.omit(prid) 
    pt.poly <- pts[as.numeric(as.character(row.names(ptid))),]  

plot(pts)
  axis(1); axis(2)
    plot(ply, add=TRUE, border="red")
      plot(pt.poly,pch=19,add=TRUE) 
share|improve this answer
    
Wrong - the rownames of the returned object corresponds to the row index in_this_case - in general the row names seem to be the row names of the points - which might not even be numeric. You could modify your solution to do a character match which might make it a bit more robust. –  Spacedman Jun 19 '13 at 23:25
    
@Sapcedman, Don't be so dogmatic. The solution is not incorrect. If you want to subset points to a set of polygons or assign polygon values to points the over function works without coercion. There are multiple was to perform the crosswalk once you have the resulting over object. You solution of coercing to a SpatialPolygon object creates considerable necessary overhead because this operation can be performed directly on the SpatialPolygonDataFrame object. By the way before you edit a post be sure that you are right. The term library and package are used interchangeably in R. –  Jeffrey Evans Jun 20 '13 at 1:03
    
I've added some benchmarks to my post, and spotted another problem with your function. Also "Packages are collections of R functions, data, and compiled code in a well-defined format. The directory where packages are stored is called the library" –  Spacedman Jun 20 '13 at 7:57
    
While you are technically correct regarding "package" vs. "library" your are arguing semantics. I just had an Ecological Modelling editor request that we change our use of "package" (which is actually my preference) to "library". My point is that they are becoming interchangeable terms and a matter of preference. –  Jeffrey Evans Jun 20 '13 at 16:54
1  
"technically correct" as Dr Sheldon Cooper once remarked, "is the best sort of correct". That editor is technically wrong, which is the worst sort of wrong. –  Spacedman Jun 20 '13 at 17:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.