将矢量图像置于网格中以在 R 中进行克里金法

Question

经过大量搜索、询问和编写一些代码后，我得到了在 R 的 gstat 中进行克里金法的最低限度。

我使用 4 个点（我知道，非常糟糕）克里格化了位于它们之间的未采样点。但实际上，我不需要所有这些要点。在那个区域内，还有一个更小的分区……这个区域才是我真正需要的。

长话短说.. 我有 4 个报告降雨数据的气象站的测量结果。这些点的纬度和经度坐标是：

lat    long
7.16   124.21
8.6    123.35
8.43   124.28
8.15   125.08

我的克里金之路可以通过我之前在 Whosebug 上的问题看到。

这个：

还有这个：

我知道图片中的坐标（至少是我估计的）：

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 和统计数据的总菜鸟。感谢所有回复的人。

我只是想知道它是否可行以及是否提供一些提示。再次感谢。

还不如post我的脚本：

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")

编辑：自上次以来代码已更改。但我想这并不重要，我只是想知道它是否可能并且任何人都可以提供它可能的最简单的例子。我可以从那里得出我自己问题的示例的解决方案。

Answer 1

为了简化您的问题：

您想根据没有地理参考的图像描绘一个区域。
您只想提取该区域的插值结果

需要几个步骤

您需要使用 QGIS 对图像进行地理配准 (Raster > Georeferencer)。您需要在后台有一个地理参考地图来提供帮助。这将创建一个具有空间信息的栅格对象。
两种可能。
- 2.a。图像的中心部分有一种颜色，可以直接用作 R 中的遮罩（例如，红色像素中间的所有绿色像素）。
- 2.b。如果没有，你需要用QGIS手动划出区域的多边形(Layer > Create Layer > New Shapefile > Polygon)
在 R 中导入栅格或多边形 shapefile
使用函数 raster::mask 使用栅格图像或 SpatialPolygon 提取插值值。

Answer 2

如果您觉得有用请告诉我：

"Think of it like a doughnut, you krige the the whole circular shape of the doughnut but you only need the area covered by the hole so you remove or at least disregard the values you got from the doughnut itself."

为此，您加载矢量数据：

donut <- rgdal::readOGR('/variogram', 'donut')
day_1@proj4string@projargs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # Becouse donut shape have this CRS

plot(donut, axes = TRUE, col = 3)
plot(day_1, col = 2, pch = 20, add = TRUE)

那你把'external ring'删掉，保留内幕。也表明第二个不再是洞：

hole <- donut # for keep original shape
hole@polygons[1][[1]]@Polygons[1] <- NULL
hole@polygons[1][[1]]@Polygons[1][[1]]@hole <- FALSE

plot(hole, axes = TRUE, col = 4, add = TRUE)

然后检查哪些点在 'hole' 新的蓝色矢量图层中：

over.pts <- over(day_1, hole)
day_1_subset <- day_1[!is.na(over.pts$Id), ]

plot(donut, axes = TRUE, col = 3)
plot(hole, col = 4, add = TRUE)
plot(day_1, col = 2, pch = 20, add = TRUE)
plot(day_1_subset, col = 'white', pch = 1, cex = 2, add = TRUE)

write.csv(day_1_subset@data, 'myfile.csv') # write intersected points table
write.csv(as.data.frame(coordinates(day_1_subset)), 'myfile.csv') # write intersected points coords
writeOGR(day_1_subset, 'path', 'mysubsetlayer', driver = 'ESRI Shapefile') # write intersected points shape

使用此代码，如果您已经拥有 shapefile，则可以解决 'ring' 或甜甜圈 'hole'。如果您有图像并想剪辑它，请尝试以下操作：

如果您加载栅格（从网络获取底图图像）：

coordDf <- as.data.frame(coordinates(day_1)) # get basemap from points
# coordDf <- data.frame(hole@polygons[1][[1]]@Polygons[1][[1]]@coords) # get basemap from hole
colnames(coordDf) <- c('x', 'y') 
imag <- dismo::gmap(coordDf, lonlat = TRUE)
myimag <- raster::crop(day_1.kriged, hole)
plot(myimag)
plot(day_1, add = TRUE, col = 2)

如果您使用 day_1.kriged:

myCropKrig<- raster::crop(day_1.kriged, hole)

  myCropKrig %>% 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) +
  geom_point(data=coordDf[!is.na(over.pts$Id), ], aes(x=x, y=y), color="blue", size=3, shape=20) +
  theme_bw()

和"Finally, I take the values of the area covered by the image and store them into a csv or database."

write.csv(as.data.frame(myCropKrig), 'myCropKrig.csv')

希望你觉得这有用，我会回应你的意思

将矢量图像置于网格中以在 R 中进行克里金法

Underlay a vector image to a grid for kriging in R

r

sp

kriging

gstat

geostatistics