Как замаскировать или обрезать результаты IDW

Мой вопрос: как я могу замаскировать или обрезать результаты интерполяции IDW с использованием R так, чтобы они показывали только ту область, которая содержит оригинальный набор точек данных?

В приведенном ниже примере 20 случайных точек используются для интерполяции поверхности с помощью функции IDW в gstat. Конвексная оболочка получается из набора точек и отображается на интерполированной карте. Но я хочу обрезать карту так, чтобы на ней отображались только области внутри облака точек.

Спасибо за любые предложения!

###
#                             IDW в 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)
#

Вот что бы я сделал:

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

# Добавление растра IDW
# Первый шаг для включения результатов IDW на интерактивную карту – преобразовать его в пространственный объект
idw_raster <- rasterFromXYZ(idw_output)

# Мы хотим показывать результаты интерполяции только в области, для которой у нас есть точки данных
# Для этого мы обрежем растр по региону данных, создав конвексную оболочку, буфер и маску (обрезку) растра
# Выполните операции конвексной оболочки и буфера на проецированных данных UTM
xyz_hull <- gConvexHull(xyz)
xyz_buff <- gBuffer(xyz_hull, width = 25000) # произвольный буфер 1 км

# Маскируем растр IDW, чтобы ограничить его регионами, для которых у нас есть информация
idw_raster_crop <- mask(idw_raster, xyz)

spplot(idw_raster_crop)

# Создаем контуры IDW
idw_contour <- rasterToContour(idw_raster_crop, nlevels = 15)

Когда вы будете использовать свои реальные данные, вы можете использовать mapview, т.е. mapview(idw_raster_crop), чтобы увидеть интерактивную карту.

Ссылка: https://rstudio-pubs-static.s3.amazonaws.com/212793_f130ecc723da4ed98680d6d3d5c4aff9.html

Здесь я отвечаю на свой собственный вопрос:

Чтобы получить интерполяцию IDW поверхности, ограниченную областью, содержащей оригинальный набор точек данных xyz:

• Создайте прямоугольную сетку на основе диапазонов точек xy с использованием функции expand.grid {base}
• Создайте конвексную оболочку с использованием функции convexhull.xy{spatstat}
• Вырежьте неправильную сетку, используя функцию inout{splancs} на
  прямоугольной сетке с маской, заданной конвексной оболочкой
• Выполните интерполяцию, используя функцию idw{gstat} с новой сеткой.
• Не забудьте отключить spatstat, потому что у него есть другая функция для
  IDW с тем же именем

Код:

###
# IDW в R
#
# Загрузите gstat:
library(gstat)
# 
# Пример набора данных XYZ:
set.seed(1234)
x <- rnorm(20)*1000+3000
y <- rnorm(20)*1000+3000
z <- rnorm(20)+10
xyz <- data.frame(x,y,z)
# 
# Рассчитайте диапазоны и расстояние:
xr <- round(range(x),1);xr
yr <- round(range(y),1);yr
b <- round((xr[2]-xr[1])/60,2);b
#
# Создайте сетку:
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)
#
# Загрузите spatstat и splancs
library(spatstat)
library(splancs)
#
# Извлечение конвексной оболочки 
w2 <- convexhull.xy(x,y)
#
# отключите spatstat, потому что у него есть другая функция "idw":
detach("package:spatstat", unload=TRUE)
#
# Создайте новую обрезанную сетку:
grd2 <- grd1[inout(grd1,w2$bdry[[1]]),]
plot(grd2)
#
# Преобразуйте сетки и набор данных в sp-объекты:
coordinates(grd1) <- ~x+y
gridded(grd1)   <- TRUE
coordinates(grd2) <- ~x+y
gridded(grd2)   <- TRUE
coordinates(xyz) <- ~x+y
#
# Запустите интерполяции:
idw1<-idw(formula=z~1, locations=xyz, newdata=grd1)
idw2<-idw(formula=z~1, locations=xyz, newdata=grd2)
#
# Постройте результаты:
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)
#