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.

pts <- data.frame(x = rnorm(100), y = rnorm(100), var1 = runif(100), var2 = sample(letters, 100, replace = TRUE))
coordinates(pts) <- ~ x + y
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)
1        NA
2       357
3       357
4        NA
5       357
6       357
  • 1
    Thanks for providing a reproducible example. Always helps when trying to understand a problem!
    – fdetsch
    Jun 18, 2013 at 16:00

6 Answers 6


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

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:

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))),]

rowlings <- function(pts,ply){

rowlings2 <- function(pts,ply){
  class(ply) <- "SpatialPolygons"

obrien <- function(pts,ply){

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


Looks good


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

  • Great, this also works for multiple polygons (I've added an example to play with into Joshua's answer). Jun 19, 2013 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. Jun 19, 2013 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, 2013 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? Jun 20, 2013 at 16:59
  • 1
    We didn't know that until I tested it just now.
    – Spacedman
    Jun 20, 2013 at 17:37

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


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. Feb 21, 2015 at 15:35
  • 2
    Welcome to our site, Edzer--it's nice to see you here!
    – whuber
    Feb 23, 2015 at 19:34
  • 1
    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! ;-) Feb 23, 2015 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.


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))),]  

  axis(1); axis(2)
    plot(ply, add=TRUE, border="red")
  • 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, 2013 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. Jun 20, 2013 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, 2013 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. Jun 20, 2013 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, 2013 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.

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

## Confirm that it works
#              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(ply, border="red", add=TRUE)
plot(pp, col="red", add=TRUE)
  • 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, 2013 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. Jun 20, 2013 at 17:27
  • Ah, any. That might be marginally faster than the version I just benchmarked.
    – Spacedman
    Jun 20, 2013 at 17:38
  • @Spacedman -- From my quick tests, it looks like obrien and rowlings2 run neck and neck, with obrien maybe 2% faster. Jun 20, 2013 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, 2016 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

# 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, ]

              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
  • 1
    Should tmp be pts.intersect? Also, parsing the returned dimnames like that is depending on undocumented behaviour.
    – Spacedman
    Jun 20, 2013 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, 2013 at 7:57

There is an extremely simple solution using the spatialEco library.


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

# check plot
  plot(a, add=T)

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

Check the result:


>             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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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