20 Jun 2013

Spatial Overlays with R - Retrieving Polygon Attributes for a Set of Points

A short tutorial for spatial overlays using R-GIS..

library(sp)
library(dismo)
 
# spatial data (political districts of Austria)
gadm <- getData('GADM', country = "AT", level = 2)

# view
plot(gadm)
 
# some addresses
pts <- geocode(c("Aldrans, Grubenweg", "Wien, Stephansdom", "Salzburg, Mozartplatz"))
 
# make pts spatial
coords <- SpatialPoints(pts[, c("longitude", "latitude")])
spdf_pts <- SpatialPointsDataFrame(coords, pts)

# assign CRS/projection (which is WGS 1984)
crs <- CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") 
proj4string(spdf_pts) <- crs
 
# check data
str(spdf_pts)
str(gadm)
 
# plot pts on top
plot(spdf_pts, cex = 2, col = 2, add = T)
 
# do an intersection (points in polygon)
# yielding the polygon's attribute data
over(spdf_pts, gadm)

1 comment :

  1. FYI I need to change the code to
    coords <- SpatialPoints(pts[, c("lon", "lat")])
    to make it work. Very cool though!

    ReplyDelete