注:これを行うためのよりエレガントな方法があるかどうかを知るために、これを部分的に投げています。ですから、空間タイプ、改善のための提案を提案してください。
(特に、値が抽出されるポイントでグリッドを設定するステップ2は、"SpatialPoints"
私には常に痛々しいほど低レベルに見えます。)
これは、その目的のために構築されたオブジェクトに含まれる座標でover()
属性を抽出するために使用されます。"SpatialPolygonDataFrame"
"SpatialPoints"
library(rgdal)
## (1) Read in an example shapefile
dsn <- system.file("vectors", package = "rgdal")[1]
scot_BNG <- readOGR(dsn=dsn, layer="scot_BNG")
scot_BNG <- scot_BNG[1:5,] # Let's just use part of it
## (2) Set up a SpatialPoints object with the grid of points
## for which you want to extract values
res <- 10000 ## Distance between grid points (30 in OP's question)
BB <- bbox(scot_BNG)
BB <- res*round(BB/res) ## Pretty up the bounding box
GT <- GridTopology(cellcentre.offset = BB[,1],
cellsize = c(res, res),
cells.dim = (c(diff(BB[1,]), diff(BB[2,]))/res) + 1)
SP <- SpatialPoints(GT, proj4string = CRS(proj4string(scot_BNG)))
## (3) Extract the values
vals <- over(SP, scot_BNG)
res <- cbind(coordinates(SP), vals)
## Finally, have a look at a few of the points.
x <- res[!is.na(res$SP_ID),]
rbind(head(x,3), tail(x,3))[1:10]
# x y SP_ID NAME ID_x COUNT SMR LONG LAT PY
# 4 230000 970000 0 Sutherland 12 5 279.3 58.06 4.64 37521
# 5 240000 970000 0 Sutherland 12 5 279.3 58.06 4.64 37521
# 25 220000 960000 0 Sutherland 12 5 279.3 58.06 4.64 37521
# 425 260000 780000 4 Bedenoch 17 2 186.9 57.06 4.09 27075
# 426 270000 780000 4 Bedenoch 17 2 186.9 57.06 4.09 27075
# 427 280000 780000 4 Bedenoch 17 2 186.9 57.06 4.09 27075