Interpolate temperature values (with associated lat and lon values) at different depth intervals using Krig function in the 'fields' R package

Question

I have temperature data points in different depth intervals with associated lat long values across the study area. I want to make a raster and then interpolate between cells of the raster where there is no data. I can do it using Krig in the fields package, but wonder whether there are better approaches? The data points are irregularly spaced and we want to take space into account. For each depth interval, we want to create separate rasters.

This is an example of what my data looks like:

# A tibble: 21 x 8
   date.time           id                   lon.x    lat.y depthbin1 depthbin2 depthbin3 depthbin4
   <dttm>              <chr>                <dbl>    <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1 2018-12-09 23:09:44 Kopaitic_Inc1_KI04 451144. 2985192.    -0.7      -0.742    -0.838    NA    
 2 2018-12-09 23:12:25 Kopaitic_Inc1_KI04 451076. 2985416.    -0.9      NA        NA        NA    
 3 2018-12-09 23:13:15 Kopaitic_Inc1_KI04 451054. 2985489.    -0.51     -0.546    -0.595    -0.622
 4 2018-12-09 23:16:00 Kopaitic_Inc1_KI04 450985. 2985731.    -0.474    -0.525    -0.575    -0.645
 5 2018-12-09 23:17:56 Kopaitic_Inc1_KI04 450940. 2985903.    -0.6      NA        NA        NA    
 6 2018-12-09 23:18:36 Kopaitic_Inc1_KI04 450926. 2985962.    -0.544    -0.526    -0.592    -0.639
 7 2018-12-09 23:21:39 Kopaitic_Inc1_KI04 450870. 2986226.    -0.6      -0.595    -0.627    -0.665
 8 2018-12-09 23:25:10 Kopaitic_Inc1_KI04 450820. 2986512.    -0.5      -0.526    -0.567    -0.576
 9 2018-12-09 23:29:41 Kopaitic_Inc1_KI04 450777. 2986829.    -0.4      -0.405    -0.512    -0.610
10 2018-12-09 23:32:19 Kopaitic_Inc1_KI04 450763. 2986985.    -0.896    NA        NA        NA    
# ... with 11 more rows
> 

There are date, id, longitude and latitude variables. And then the mean temperature was measured by a device in every depthbin if the animal dived to that depth interval. If the animal didn't dive that deep, the depthbin value is empty.

This is how I am interpolating at the moment:

# Make a raster layer
library(raster)

# projection
utm.prj = " +proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs "   

# create a SpatialPointsDataFrame
coordinates(divetemps) = ~lon.x+lat.y   
proj4string(divetemps) <-CRS(utm.prj)

# create an empty raster object to the extent of the points
rast <- raster(ext=extent(divetemps),crs = CRS(utm.prj), resolution = 500) # 500 m x 500 m
rast

# rasterize your irregular points 
rasOut<-raster::rasterize(divetemps, rast, divetemps$depthbin1, fun = mean) # we use a mean function here to regularly grid the irregular input points
plot(rasOut)

library(fields)
# Function to Krig 
krigR <- function(rast){ 
  xy <- data.frame(raster::xyFromCell(rast, 1:ncell(rast)))
  v <- getValues(rast)
  krg <- fields::Krig(xy, v) 
  ras.int <- raster::interpolate(rast, krg)
  proj4string(ras.int) <- proj4string(rast)
  return(ras.int)
}

surface = krigR(rasOut) 
plot(surface)

This is an example of the plots that I get when using the fields::Krig function to interpolate the temperature values for depth bin 1 across the whole study area: Interpolated temperature values over study area

I am not entirely happy with the plots that I am getting when using the fields::Krig function because I don't know how accurate they are. I know there is not a big difference in temperatures across the study area. But I am sure that my plots can look better than this.

So I would like to try out other R packages and functions to interpolate temperature values across a study area. Does anyone have any suggestions of functions or packages that I can look into and try out that you have used before?

Answer 1

It is very difficult to help if you do not provide example data (as data, not by printing the data). But your approach does not look right to me. Why do you use rasterize? Here is a simplified workflow:

library(raster)
library(fields)

utm.prj = " +proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs "   
xy <- divetemps[, c("lon.x", "lat.y")]
rast <- raster(ext=extent(xy)+1000, crs=utm.prj, resolution = 500)

m <- fields::Krig(xy, divetemps$depthbin1)
surface <- raster::interpolate(rast, m)

plot(surface)

I think this is conceptually better, but you may still not like the output. I would expect that you would prefer the output of at Thin Plate Spline model:

m <- fields::Tps(xy, divetemps$depthbin1)
surface <- raster::interpolate(rast, m)

And you may consider to not use only x and y but also an additional variable such as elevation to predict temperature. See ?interpolate