在 R 中的空间数据框中查找最近的网格值答案

【问题标题】：Finding nearest grid value in spatial dataframe in R在 R 中的空间数据框中查找最近的网格值
【发布时间】：2015-07-17 18:29:42
【问题描述】：

我正在学习空间数据框。我创建了一个简单的示例，它使用反距离加权 (IDW) 根据一组初始点在网格中插入值。我现在想获得每个初始点的插值。我认为这应该相对容易，但很难找到一种优雅的方式。换句话说，我想在“geog2.idw”中为 40 个点 (x,y) 中的每一个找到最近的网格值。

代码

#adjusted from example created by Brian J. Knaus

#packages
library(gstat)
require(gstat)
require(lattice)
require(raster)
require(sp)

##### ##### ##### ##### #####
# Generate data.
set.seed(1)
x <- runif(40, -123, -110)
set.seed(2)
y <- runif(40, 40, 48)
set.seed(3)
z <- runif(40, 0, 1)

# Create spatial data.frame.
geog2 <- data.frame(x,y,z)

# Assigning coordinates results in spdataframe.
coordinates(geog2) = ~x+y


#create grid
pixels <- 100

geog.grd <- expand.grid(x=seq(floor(min(coordinates(geog2)[,1])),
                              ceiling(max(coordinates(geog2)[,1])),
                              length.out=pixels),
                        y=seq(floor(min(coordinates(geog2)[,2])),
                              ceiling(max(coordinates(geog2)[,2])),
                              length.out=pixels))
#
grd.pts <- SpatialPixels(SpatialPoints((geog.grd)))
grd <- as(grd.pts, "SpatialGrid")


##### ##### ##### ##### #####
# IDW interpolation.

geog2.idw <- idw(z ~ 1, geog2, grd, idp=6)

spplot(geog2.idw["var1.pred"])

##### ##### ##### ##### #####
# Image.

par(mar=c(2,2,1,1))
#
image(geog2.idw, xlim=c(-125, -110), ylim=c(40, 47))
#
contour(geog2.idw, add=T, lwd=0.5)
#
map("state", lwd=1, add=T)
map("county", lty=3, lwd=0.5, add=T)
#
points(x, y, pch=21, bg=c(1:4))

#######

【问题讨论】：

我收到一个错误：“错误：找不到函数“map””
你可以把x和y变成SpatialPoints:SpatialPoints(cbind(x,y))并使用函数over。

标签： r geospatial raster gstat

【解决方案1】：

你可以把x和y变成SpatialPoints:SpatialPoints(cbind(x,y))并使用函数over:

pts <- SpatialPoints(cbind(x, y))
pts.idw <- over(pts, geog2.idw["var1.pred"])

# Writing extracted values on the map:
image(geog2.idw, xlim=c(-125, -110), ylim=c(40, 47))
contour(geog2.idw, add=T, lwd=0.5)
map("state", lwd=1, add=T)
map("county", lty=3, lwd=0.5, add=T)
text(pts, round(pts.idw[,"var1.pred"],4), col = "blue")

【讨论】：