Follow Up on Spatial Overlays with R Getting Altitude for a Set of Points
A short follow up on a previous post on spatial overlays with R.
library(sp)
library(dismo)
# some addresses in Austria
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
# spatial data to extract (altitude)
alt <- getData(alt, country = "AT")
# convert alt from raster to grid (needed for over::sp)
# and assign CRS (which is the same as spdf_pts, see > alt@crs)
# dont mind warning - the CRS is the same..
spdf_alt <- as(alt, SpatialGridDataFrame)
proj4string(spdf_alt) <- crs
# view
plot(alt)
# plot pts on top
plot(spdf_pts, cex = 2, col = 2, add = T)
# check data
str(spdf_pts)
str(spdf_alt)
# get the raster/pixel/grid data (> ?over):
cbind(spdf_pts$interpretedPlace, over(spdf_pts, spdf_alt))
# result:
# spdf_pts$interpretedPlace AUT_msk_alt
# 1 Grubenweg, 6071 Aldrans, Austria 736
# 3 Saint Stephens Vienna, Stephansplatz 1, 1010 Vienna, Austria 183
# 2 Mozartplatz, 5020 Salzburg, Austria 450