将矢量图像置于网格中以便在R中进行克里金法(Underlay a vector image to a grid for kriging in R)

将矢量图像置于网格中以便在R中进行克里金法(Underlay a vector image to a grid for kriging in R)

经过大量搜索，询问并做了一些代码之后，我在R's 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点kriged区域内插入该矢量图像，因此我知道要实际考虑哪些坐标以及要删除哪些坐标。

最后，我获取图像所覆盖区域的值，并将它们存储到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")

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

After searching around a lot, asking, and doing some code, I kinda got the bare minimum for doing kriging in R's gstat.

Using 4 points (I know, totally bad), I kriged the unsampled points located between them. But in actuality, I don't need all of those points. Inside that area, there is a smaller subarea... this area is the one I actually need.

Long story short.. I have measurements taken from 4 weather stations that report rainfall data. The lat and long coordinates for these points are:

lat long 7.16 124.21 8.6 123.35 8.43 124.28 8.15 125.08

My road to kriging can be seen through my previous questions on StackOverflow.

This: Create variogram in R's gstat package

And this: Create Grid in R for kriging in gstat

I know that the image in has the coordinates (at least according to my estimates):

Leftmost: 124 13ish 0 E(DMS) Rightmost : 124 20ish 0 E Topmost corrdinates: 124 17ish 0 E Bottommost coordinates: 124 16ish 0 E

Conversion will take place for that but that doesn't matter I think, or easier to deal with later.

The image is also irregular (but aren't they all though).

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.

I have an image (.jpg) of the area in question, I will have to convert the image into a shapefile or some other vector format using QGIS or similar software. After that, I will have to insert that vector image inside the 4 point kriged area, so I know which coordinates to actually consider and which ones to remove.

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

Anybody know how I can start with this? Total noob at R and statistics. Thanks to anyone who responds.

I just want to know if its possible and if it is provide some tips. Thanks again.

Might as well also post my script:

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

EDIT: The code has changed since the last time. But that doesn't matter I guess, I just want to know if its possible and can anybody provide the simplest example of it being possible. I can derive the solution to the example to my own problem from there.

最满意答案

如果您发现这有用，请告诉我：

“把它想象成一个甜甜圈，你可以看到甜甜圈的整个圆形形状，但你只需要这个洞所覆盖的区域，这样你就可以移除或者至少忽略你从甜甜圈本身获得的价值。”

为此，您加载矢量数据：

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)

在那之后你会在'洞'新的蓝色矢量图层中找到whicch点：

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'或donut'hole'。 如果您有图像并想剪辑它，请尝试以下操作：

在您加载栅格的情况下（从Web获取底图图像）：

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

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

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

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

Let me know if you find this useful:

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

For this you load your vectorial data:

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)

Then you delete the 'external ring' and keep the insider. Also indicates that the second isn't a hole anymore:

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)

After that you chek whicch points are inside 'hole' new blue vector layer:

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

With this code you can solve the 'ring' or doughnut 'hole' if you already have the shapefile. If you have an image and want to clip it try the follow:

In the case you load a raster (get basemap image from web):

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)

In case you use 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()

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

Hope you find this useful and I respond your meaning