CODEWITHSUNDEEP
rgis
Bastian
Bastian

Reputation: 78

How to obtain SpatialPolygonsDataframe that contains coordinates

I have a SpatialPolygonsDataFrame of the Spanish census tracts and a set of coordinates.

What I would like to do is to obtain the subset of census tract that contains a pair of coordinates (as a new SpatialPolygonsDataFrame).

I came closest using 'over' following this example: https://www.nceas.ucsb.edu/scicomp/usecases/point-in-polygon But 'over' does not return a SpatialPolygonsDataFrame. seccion contains columns "Long" and "Lat" with the coordinates. map is the SpatialPolygonsDataFrame

# turn seccion into a SpatialPoints
coordinates(seccion) <- c("Long", "Lat")
# the coordinates are in the same lat/lon reference system (as is "map")
proj4string(seccion) <- proj4string(map)
# use 'over' with map as a SpatialPolygonsDataFrame
# object, to determine which section (if any) contains the coordinates
sec <- over(seccion, map)

I tried other things and also had a look at the raster library but until now, no luck...

Upvotes: 0

Views: 223

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47696

You can use the intersect function from the raster package

Example data:

library(raster)
pols <- shapefile(system.file("external/lux.shp", package="raster"))
pnts <- SpatialPoints(matrix(c(6, 5.91, 6.13, 6.23, 6.18, 
                       50.08, 49.95, 49.63, 49.64, 49.59),
                       ncol=2), proj4string=crs(pols))

x <- intersect(pols, pnts)

check out the result:

plot(pols, col='light gray')
plot(x, col='blue', add=TRUE)
points(pnts, pch=20, col='red')

It is possible to get there via over (note the order of the arguments)

v <- over(pols, pnts)
xx <- pols[which(!is.na(v)), ]

Upvotes: 1

Related Questions