Saturday, September 19, 2015

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)

No comments:

Post a Comment