how to overlay a polygon over SpatialPointsDataFrame and preserving the SPDF data?

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
• Thanks for providing a reproducible example. Always helps when trying to understand a problem! – fdetsch Jun 18 '13 at 16:00

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"))
 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
  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
  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
  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)
 100
> pts = pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),]
> nrow(pts)
 54
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:

Two functions - first Jeffrey Evans' method, then my original, then my hacked 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))
 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))
 FALSE
> dim(evans(pts,columbus))
 27  1
> dim(rowlings(pts,columbus))
 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...

• 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
• We didn't know that until I tested it just now. – Spacedman Jun 20 '13 at 17:37

sp provides a shorter form to select features based on spatial intersection, following the OP example:

pts[ply,]

as of:

points(pts[ply,], col = 'red')

Behind the scenes this is short for

pts[!is.na(over(pts, geometry(ply))),]

The thing to note is that there is a geometry method that drops attributes: over changes behaviour if its second argument has attributes or not (this was OP's confusion). This works accross all Spatial* classes in sp, although some over methods require rgeos, see this vignette for details, e.g. the case of multiple matches for overlapping polygons.

• Good to know! I was not aware of the geometry method. – Jeffrey Evans Feb 21 '15 at 15:35
• Welcome to our site, Edzer--it's nice to see you here! – whuber Feb 23 '15 at 19:34
• Thanks Bill -- it's getting more quiet on stat.ethz.ch/pipermail/r-sig-geo , or maybe we should develop software that causes more trouble! ;-) – Edzer Pebesma Feb 23 '15 at 20:01

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)
• 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
• "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

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)
• 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
• @JoshO'Brien how can one use this answer on many polygons? That is pp should have a ID which indicates in which polygon the points are located. – code123 Apr 19 '16 at 4:57

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)[], " ")
pts.intersect.id <- as.numeric(sapply(pts.intersect.strsplit, "[[", 2))

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

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

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
• 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 :-) – fdetsch Jun 21 '13 at 7:57

There is an extremely simple solution using the spatialEco library.

library(spatialEco)

# intersect points in polygon
pts <- point.in.poly(pts, ply)

# check plot
plot(ply)

# convert to data frame, keeping your data
pts<- as.data.frame(pts)

Check the result:

pts

>             x          y       var1 var2 polyvar
> 2  -0.5832050 -0.8777367 0.04001092    v     357
> 3   0.3947471  0.7020481 0.58108350    v     357
> 5   0.7667997 -0.9465043 0.85682609    q     357
> 6   0.3174604  0.6416281 0.13683264    y     357
> 9  -0.4690151  0.4413502 0.13968804    m     357
> 10  0.4765213  0.6068021 0.97144627    o     357