在 R{Spatstat} 中将 IDW 插值限制为多边形答案

【问题标题】：Confining an IDW interpolation to a polygon in R{Spatstat}在 R{Spatstat} 中将 IDW 插值限制为多边形
【发布时间】：2016-02-25 03:16:52
【问题描述】：

我是一个相当新的 R 用户，使用嵌入在美国都会区的地方的空间数据。我正在为许多变量生成 IDW 地图，但在定义 IDW 的边界时遇到了麻烦。我想只在都会区边界范围内插入数据，而我目前的方法会填充整个地块。

我首先从dataframe 创建一个ppp 对象，如下所示：

x.points<-SpatialPointsDataFrame(coords =x.data[,c("LON10","LAT10")],
  data=x.data[,c("v1","v2")],
  proj4string=CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"))

x.proj<-spTransform(x.points,CRSobj = CRS("+proj=aea
  +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 
  +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))    

x.ppp<-as.ppp(x.proj)

使用plot(x.ppp) 投影此ppp 对象会显示每个变量（上面代码中的“v1”和“v2”）的图。但是，绘图填充了整个矩形框。在绘制我的ppp 对象时，我通过为它分配一个owin 对象来限制这个框架，如下所示：

y.spdf<-spTransform(y.spdf,CRSobj = CRS("+proj=aea 
  +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 
  +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))

y.owin<-as.owin(y.spdf)

x.ppp<-x.ppp[y.owin]

使用plot(x.ppp) 绘制此图现在可以为我绘制一个正确框在其都市区多边形内的数据点图。虽然这是一个好的开始，但我想继续将此空间框架传递到 IDW 地图上。查阅 spatstat 手册让我相信我可以将我的 owin 对象作为 as.mask 参数传递，但每当我这样做时，我都会收到错误消息。示例代码如下：

idw<-idw(x.ppp, power=4, as.mask(y.owin))

这给了我错误“开关错误（在，像素 = {：EXPR 必须是长度为 1 的向量”。使用 at="points" 不起作用，因为我得到的是结果 dataframe 对象而不是 @ 987654335@。如果我省略了 as.mask 参数，我会得到一个正确的 idw，但它会插入到我之前能够使用 owin 对象设置的都会区边界之外。

我认为这是一项相当常见的任务，但查阅 spatstat 文档让我有点不知所措。 as.mask 参数似乎可行，但我不确定我哪里出错了。关于如何解决这个问题的任何意见或建议，或者我可以求助于尝试更好地理解这个问题的来源，都会有很大的帮助。谢谢！

【问题讨论】：

标签： r spatial spatstat

【解决方案1】：

我明白这有点令人困惑。基本上idw 在包含数据的整个矩形框（owin 的框架）中计算结果。这是您执行rslt <- idw(x.ppp) 时得到的结果。要将结果限制为多边形窗口，我们需要通过该窗口对结果图像进行子集化（并设置drop=FALSE 以获取结果不规则图像，而不是多边形区域内的像素值）。代码解释：

简单案例（一个数字标记值）

library(spatstat)
# Fake data with one numeric mark
X <- humberside
marks(X) <- runif(npoints(X))
# Window is a polygonal region:
W <- Window(X)
W
# IDW in the frame of W
rslt <- idw(X)
# Subset of IDW inside W
rslt <- rslt[W, drop = FALSE]
plot(rslt)

多变量大小写（两个或多个数字标记值）

library(spatstat)
# Fake data with more numeric marks
X <- humberside
marks(X) <- cbind(x = runif(npoints(X)), y = runif(npoints(X)))
# Window is a polygonal region:
W <- Window(X)
W
# IDW in the frame of W for each mark variable (imlist object)
rslt <- idw(X)
# Subset of IDW inside W for each mark variable (imlist object)
rslt <- solapply(rslt, "[.im", W, drop = FALSE)
plot(rslt)

【讨论】：