Plotting spatial model predictions (issues with plot)

Question

I have created the following model and predictions but I'm having trouble with the code to plot the predictions. I think it's a dimensions issue, does anyone know the changes I need to make for this to work?

code used;

#variogram
summer_vario = variog(geo_summer_df2, option = 'bin', estimator.type='modulus', bin.cloud = TRUE)
#fitting a basic parametric model
defult_summer_mod = variofit(summer_vario)

#creating predictions
preds_grid = matrix(c(-5.697, 55.441, -0.807, 51.682, -5.328, 50.218, -2.451, 54.684, -4.121, 50.355, -1.586, 54.768, -0.131, 51.505, -4.158, 52.915, 
                  -0.442, 53.875, -3.413, 56.214, -2.860,   54.076, -3.323, 57.711, 0.566,  52.651, -0.626, 54.481, -1.185, 60.139, -2.643, 51.006,
                  -1.491,   53.381, -1.536, 52.424, -6.319, 58.213, -1.992, 51.503), nrow = 20, byrow = TRUE)
summer_preds = krige.conv(geo_summer_df2, locations = preds_grid, krige = krige.control(obj.model = defult_summer_mod))

#plotting predictions
#mean
image(summer_preds, col = viridis::viridis(100), zlim = c(100, max(c(summer_preds$predict))), 
  coords.data = geo_summer_df2[1]$coords, main = 'Mean', xlab = 'x', ylab = 'y',
  x.leg = c(700, 900), y.leg = c(20, 70))
#variation
image(summer_preds, values = summer_preds$krige.var, col = heat.colors(100)[100:1],
  zlim = c(0,max(c(summer_preds$krige.var))), coords.data = geo_summer_df2[1]$coords,
  main = 'Variance', xlab = 'x', ylab = 'y', x.leg = c(700, 900), y.leg = c(20, 70))

data used;

https://drive.google.com/file/d/1ngwto6hgqCumoDsStOtPoG2J5EbmqxDf/view?usp=sharing https://drive.google.com/file/d/1s9yBHsgaFRlF38CgiXCf_vum1DyhEbz4/view?usp=sharing

data changes made before code at the top of the page

#converting data to long format and combining both dataframes
MaxTemp %>%
 pivot_longer(.,Machrihanish:Lyneham, names_to = "Location") %>%
 full_join(.,metadata) -> MaxTemp_df

#renaming value column to temperature
MaxTemp_df = MaxTemp_df %>%
 rename(Temp = 'value')

#filtering data for summer months
summer_df = MaxTemp_df %>% 
  filter(Date >= 20200701 & Date <=20200731)

 #converting our data to geodata
 geo_summer_df = as.geodata(summer_df, coords.col = 4:5, data.col = 3) 
geo_summer_df2 = jitterDupCoords(geo_summer_df, max = 0.1, min = 0.05)

Answer 1

You're right about the dimensions. The predictions should be made over a regular grid of locations if you want to plot them as an image. Get all the unique x co-ordinates and all the unique y co-ordinates, sort them, then use expand.grid to get x, y co-ordinates for the whole grid. You'll then need to use this for kriging.

When you come to drawing the image, you need to arrange the predictions into a matrix:

xvals <- sort(unique(preds_grid[,1]))
yvals <- sort(unique(preds_grid[,2]))
preds_grid <- as.matrix(expand.grid(xvals, yvals))
colnames(preds_grid) <- NULL

summer_preds = krige.conv(geo_summer_df2, locations = preds_grid, 
                          krige = krige.control(obj.model = default_summer_mod))

image(xvals, yvals, matrix(summer_preds$predict, nrow = length(xvals)),
      col = viridis::viridis(100), main = 'Mean', xlab = 'x', ylab = 'y')

image(xvals, yvals, matrix(summer_preds$krige.var, nrow = length(xvals)), 
      col = heat.colors(100)[100:1], main = 'Variance', xlab = 'x', ylab = 'y')

Note that you will get better images if you use a finely-spaced sequence for x and y:

xvals <- seq(-7, 1, 0.1)
yvals <- seq(50, 62, 0.1)

The plots this produces with the same code otherwise are:

---

Update - using ggplot

The following adds the data to an outline of the British Isles:

devtools::install_github("ropensci/rnaturalearthhires")
library(rnaturalearth)

xvals <- seq(-7, 1, 0.1)
yvals <- seq(50, 62, 0.1)
preds_grid <- as.matrix(expand.grid(xvals, yvals))

summer_preds <- krige.conv(
  geo_summer_df2, locations = preds_grid, 
  krige = krige.control(obj.model = default_summer_mod))

df <- as.data.frame(cbind(preds_grid, 
                          mean = summer_preds$predict, 
                          var = summer_preds$krige.var))

gb <- sf::st_crop(ne_coastline(scale = 10, returnclass = 'sf'), 
                  xmin = -7, xmax = 1, ymin = 50, ymax = 62)

ggplot(gb) + 
  geom_tile(data = df, aes(Var1, Var2, fill = mean),
            width = 0.11, height = 0.11, size = 0) +
  geom_sf() +
  scale_fill_viridis_c() +
  ggtitle('Mean')

ggplot(gb) + 
  geom_tile(data = df, aes(Var1, Var2, fill = var),
            width = 0.11, height = 0.11, size = 0) +
  geom_sf() +
  scale_fill_gradientn(colors = heat.colors(100, rev = TRUE)) +
  ggtitle('Variance')