Underlay a vector image to a grid for kriging in R

Question

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.

Answer 1

To simplify your question:

• You want to delineate an area based on an image that is not georeferenced.
• You want to extract results of a interpolation only for this area

Few steps are required

1. You need to use QGIS to georeference your image (Raster > Georeferencer). You need to have a georeferenced map in background to help. This creates a raster object with spatial information.
2. Two possibilities.
   • 2.a. The central part of your image has a color than can be directly used as a mask in R (Ex. All green pixels in middle of red pixels).
   • 2.b. If not, you need to use QGIS to delineate manually a Polygon of the area (Layer > Create Layer > New Shapefile > Polygon)
3. Import your raster or polygon shapefile in R
4. Use function raster::mask to extract values of your interpolation using the raster image or the SpatialPolygon.

Answer 2

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