如何掩蔽或裁剪IDW结果

Question

我的问题是:我如何能够掩蔽或裁剪结果，使用R，IDW插值的结果，只包含原始的一组数据点？

在下面的示例中，使用20个随机点在gstat中使用IDW函数插值曲面。从点集得到凸包，并在插值映射上绘制。但是，我想要裁剪地图，所以只有点云内的区域显示插值结果。

谢谢你的建议！

###
#                             IDW in R
#
library(gstat)
#             
set.seed(1234)
x <- rnorm(20)+10
y <- rnorm(20)+10
z <- rnorm(20)+10
#
xyz <- data.frame(x,y,z)
coordinates(xyz) <- ~x+y
xr <- range(x)
yr <- range(y)
b <- round((xr[2]-xr[1])/100,2)
#
grd <- expand.grid(x=seq(from=xr[1],to=xr[2],by=b),
y=seq(from=yr[1],to=yr[2],by=b))
#
coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE
#
idw<-idw(formula=z~1, locations=xyz, newdata=grd)
#
image(idw, col=terrain.colors(16))
contour(idw, add=T)
points(x,y, col=4, pch=19)
#
library(spatstat)
conv<-convexhull.xy(x,y)
plot(conv,add=TRUE)
#

Answer 1

我在这里回答我自己的问题：

要获得仅限于包含原始xyz数据点集的区域的IDW曲面插值：

• 使用expand.grid {base}根据xy点范围创建矩形网格
• 使用函数convexhull.xy{spatstat}创建凸包
• 在凸包蒙面的矩形网格上使用inout{splancs}函数切割不规则网格
• 使用idw{gstat}函数对新网格执行插值。
• 一定要分离spatstat，因为它对于同名的IDW具有不同的功能。

代码：

###
# IDW in R
#
# Load gstat:
library(gstat)
# 
# An example XYZ data set:
set.seed(1234)
x <- rnorm(20)*1000+3000
y <- rnorm(20)*1000+3000
z <- rnorm(20)+10
xyz <- data.frame(x,y,z)
# 
# Calculate ranges and spacing:
xr <- round(range(x),1);xr
yr <- round(range(y),1);yr
b <- round((xr[2]-xr[1])/60,2);b
#
# Create a grid:
grd1 <- expand.grid(x=seq(from=xr[1],to=xr[2],by=b),
    y=seq(from=yr[1],to=yr[2],by=b))
plot(grd1)
#
# Load spatstat and splancs
library(spatstat)
library(splancs)
#
# Extract convex hull 
w2 <- convexhull.xy(x,y)
#
# detach spatstat because it has a different "idw" function:
detach("package:spatstat", unload=TRUE)
#
# Create a new cropped grid:
grd2 <- grd1[inout(grd1,w2$bdry[[1]]),]
plot(grd2)
#
# Convert grids and dataset to sp objects:
coordinates(grd1) <- ~x+y
gridded(grd1)   <- TRUE
coordinates(grd2) <- ~x+y
gridded(grd2)   <- TRUE
coordinates(xyz) <- ~x+y
#
# Run the interpolations:
idw1<-idw(formula=z~1, locations=xyz, newdata=grd1)
idw2<-idw(formula=z~1, locations=xyz, newdata=grd2)
#
# Plot the results:
image(idw1, col=terrain.colors(16))
contour(idw1, add=T)
points(x,y, col=4, pch=19)
image(idw2, col=terrain.colors(16))
contour(idw2, add=T)
points(x,y, col=4, pch=19)
#

Answer 2

以下是我要做的事：

idw_output <- as.data.frame(idw)[, 1:3]

# Adding the IDW raster
# The first step to incorporate the IDW results into the interactive map is to turn it into a spatial object
idw_raster <- rasterFromXYZ(idw_output)

# We only want to show the interpolation results within the area of which we have data points
# To do this, we will clip the raster to data region by creating convex hull, buffer, and mask (clip) raster
# Perform convex hull and buffer operations on projected UTM data
xyz_hull <- gConvexHull(xyz)
xyz_buff <- gBuffer(xyz_hull, width = 25000) # arbitrary 1km buffer

# Mask IDW raster to confine to regions that we have information for
idw_raster_crop <- mask(idw_raster, xyz)

spplot(idw_raster_crop)

# Create contours of the IDW
idw_contour <- rasterToContour(idw_raster_crop, nlevels = 15)

使用真实数据时，可以使用mapview ( mapview(idw_raster_crop) )查看交互式地图。

引文：https://rstudio-pubs-static.s3.amazonaws.com/212793_f130ecc723da4ed98680d6d3d5c4aff9.html