3

在搜索了很多，询问并编写了一些代码之后，我有点了解在 R 的 gstat 中进行克里金法的最低要求。

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

长话短说……我从 4 个报告降雨数据的气象站进行了测量。这些点的经纬度坐标为：

lat    long
7.16   124.21
8.6    123.35
8.43   124.28
8.15   125.08

通过我之前关于 StackOverflow 的问题可以看出我的克里金之路。

这：在 R 的 gstat 包中创建变异函数

这：在 R 中为 gstat 中的克里金法创建网格

我知道图像中的坐标（至少根据我的估计）：

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

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

4

2 回答 2

2

为了简化您的问题：

  • 您想根据未进行地理参考的图像来描绘区域。
  • 您只想提取该区域的插值结果

需要几个步骤

  1. 您需要使用 QGIS 对您的图像进行地理配准 ( Raster > Georeferencer)。您需要在背景中有地理参考地图来提供帮助。这将创建一个具有空间信息的栅格对象。
  2. 两种可能。
    • 2.a. 图像的中心部分的颜色可以直接用作 R 中的蒙版（例如，红色像素中间的所有绿色像素）。
    • 2.b。如果没有，则需要使用 QGIS 手动描绘区域的多边形 ( Layer > Create Layer > New Shapefile > Polygon)
  3. 在 R 中导入您的栅格或多边形 shapefile
  4. 使用函数raster::mask使用光栅图像或 SpatialPolygon 提取插值值。
于 2017-04-21T07:31:02.620 回答
2

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

“把它想象成一个甜甜圈，你可以计算甜甜圈的整个圆形，但你只需要孔覆盖的区域，这样你就可以删除或至少忽略你从甜甜圈本身获得的值。”

为此，您加载矢量数据：

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)

第一个情节

然后你删除'外环'并保留内部人员。还表明第二个不再是一个洞：

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)

蓝色是新的 shapefile

之后，您检查哪些点在“洞”新的蓝色矢量图层内：

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，则使用此代码可以解决“环”或甜甜圈“洞”。如果您有图像并想对其进行剪辑，请尝试以下操作：

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

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

情节3

“最后，我将图像所覆盖区域的值存储到 csv 或数据库中。”

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

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

于 2017-04-24T12:55:47.727 回答