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

# 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

# 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

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

surface = krigR(rasOut) 

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?

1 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:


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)


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

  • Hi Robert, I use rasterize to get a mean temperature for each specific raster cell. Because for some raster cells we may have a lot of temperature observations, but we only need one mean temperature per raster cell. And then, if I understand correctly, the interpolate function will interpolate temperatures over the raster cells where we don't have any temperature observations. Commented Feb 1, 2022 at 14:05
  • So far, your method has not worked yet, I keep getting the error extent(xy) : c("x", "y") %in% names(x) are not all TRUE. So would you want to take a look at some data that I can post on google drive? Commented Feb 1, 2022 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.