经过很多搜索、询问和编码，我在R的gstat中实现了Kriging算法的基本功能。
使用4个点（是的，非常少），我对它们之间未采样的点进行了Kriging插值。但实际上，我并不需要所有这些点。在这个区域内，有一个更小的子区域…… 这才是我真正需要的。
长话短说……我有从4个气象站收集到的降雨数据。这些点的纬度和经度坐标如下：
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
我在StackOverflow上的之前的问题可以看出我走向Kriging的路程。
这个：在R的gstat包中创建变异图 还有这个：在gstat中为Kriging创建R网格 我知道图片中至少有以下坐标（根据我的估计）：
Leftmost: 124 13ish 0 E(DMS)
Rightmost : 124 20ish 0 E
Topmost corrdinates: 124 17ish 0 E
Bottommost coordinates: 124 16ish 0 E
转换将会进行，但我认为这并不重要，或者可以在以后更容易处理。
图像也是不规则的（但它们都是这样的）。
把它想象成一个甜甜圈，你克里格了整个甜甜圈的形状，但你只需要孔洞覆盖的区域，所以你删除或至少忽略了从甜甜圈本身得到的值。
我有一个关于该区域的图像（.jpg），我需要使用QGIS或类似软件将图像转换成shapefile或其他矢量格式。之后，我将不得不将该矢量图像插入4个点克里格区域内，以便我知道实际要考虑哪些坐标和哪些要删除。
最后，我将取出图像覆盖区域的值，并将其存储到csv或数据库中。
有人知道我如何开始吗？我对R和统计学完全不懂。感谢任何回答的人。
我只想知道它是否可能，如果可能，请提供一些提示。再次感谢。
也许我应该发布我的脚本：
suppressPackageStartupMessages({
library(sp)
library(gstat)
library(RPostgreSQL)
library(dplyr) # for "glimpse"
library(ggplot2)
library(scales) # for "comma"
library(magrittr)
library(gridExtra)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
})
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname="Rainfall Data", host="localhost", port=5432,
user="postgres", password="postgres")
day_1 <- dbGetQuery(con, "SELECT lat, long, rainfall FROM cotabato.sample")
coordinates(day_1) <- ~ lat + long
plot(day_1)
x.range <- as.integer(c(7.0,9.0))
y.range <- as.integer(c(123.0,126.0))
grid <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.05),
y=seq(from=y.range[1], to=y.range[2], by=0.05))
coordinates(grid) <- ~x+y
plot(grid, cex=1.5)
points(day_1, col='red')
title("Interpolation Grid and Sample Points")
day_1.vgm <- variogram(rainfall~1, day_1, width = 0.02, cutoff = 1.8)
day_1.fit <- fit.variogram(day_1.vgm, model=vgm("Sph", psill = 8000, range = 1))
plot(day_1.vgm, day_1.fit)
plot1 <- day_1 %>% as.data.frame %>%
ggplot(aes(lat, long)) + geom_point(size=1) + coord_equal() +
ggtitle("Points with measurements")
plot(plot1)
############################
plot2 <- grid %>% as.data.frame %>%
ggplot(aes(x, y)) + geom_point(size=1) + coord_equal() +
ggtitle("Points at which to estimate")
plot(plot2)
grid.arrange(plot1, plot2, ncol = 2)
coordinates(grid) <- ~ x + y
############################
day_1.kriged <- krige(rainfall~1, day_1, grid, model=day_1.fit)
day_1.kriged %>% as.data.frame %>%
ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
scale_fill_gradient(low = "yellow", high="red") +
scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
theme_bw()
write.csv(day_1.kriged, file = "Day_1.csv")
编辑：自上次以来，代码已经发生了变化。但我想这并不重要，我只是想知道是否可能，并且有人能否提供最简单的示例来证明它是可能的。我可以从中推导出解决方案，以解决自己的问题。
