Overlay vector image to mesh for kriging in R
After a lot of searching, querying, and running some code, I got the bare minimum to do kriging in R gstat.
Using 4 points (absolutely bad, I know), I anchored the uncovered points in between. But I really don't need all these points. There is a small sub-area inside this area ... this area I need.
Shortly speaking. I have measurements taken from 4 weather stations that report rain data. Local and long coordinates for these points:
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
You can see my road to krigings from my previous questions on StackOverflow.
These are: Create semivariogram in gstat R package
And this: Create mesh in R for kriging in gstat
I know that the image in has 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
A conversion will be done for this, but it doesn't matter I guess or is more easily dealt with later.
The image is also irregular (but still not all).
Think of it like a donut, you've drawn the whole circular shape of the donut, but you only want the area covered by the hole, so you are deleting or at least ignoring the values you got from the donut itself.
I have an image (.jpg) of an area in question, I will have to convert the image to a shapefile or some other vector format using QGIS or similar software. After that, I will need to insert this vector image inside the 4-point kriged area, so I know which coordinates to consider and which ones to delete.
Finally, I take the values of the area covered by the image and store them in a csv or database.
Does anyone know how I can get started with this? Total noob in R and statistics. Thanks to everyone who answers.
I just want to know if this is possible and if any clues are given. Thanks again.
My script can also be placed:
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 it doesn't matter, I guess I just want to know if this is possible and if someone can provide the simplest example of how this is possible. I can find a solution for my own problem.
source to share
Let me know if you find it helpful:
"Think of it like a donut, you twist the whole round donut shape, but you only want the area covered by the hole, so you remove or at least ignore the values you got from the donut itself."
To do this, you load your vector 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 remove the "outer ring" and keep the insider. Also indicates that the second is no longer a hole:
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)
Then you click on the points inside the 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 "hole" of the donut if you already have the shapefile. If you have an image and want to pin it, try this:
In case of loading a raster (get an image of the basemap from the Internet):
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)
If you are using 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 in a csv or database."
write.csv(as.data.frame(myCropKrig), 'myCropKrig.csv')
I hope you find this useful and I am answering your meaning.
source to share
To simplify your question:
- You want to outline an area based on an image that is not anchored to anchor.
- You want to extract interpolation results for this area only
Several steps required
- You need to use QGIS to link your image (
Raster > Georeferencer
). You need a georeferenced map in the background. This creates a bitmap with spatial information. - Two possibilities.
- 2.a. The center of your image has a color that can be directly used as a mask in R (for example, all green pixels are in the middle of red pixels).
- 2b. If not, you need to use QGIS to manually delineate the area polygon (
Layer > Create Layer > New Shapefile > Polygon
)
- Import a raster or polygon shape file into R
- Use a function
raster::mask
to extract the values of your interpolation using a bitmap or SpatialPolygon.
source to share