I have a set of points (not evenly distributed across space), which have location (lat and long), temperature, and elevation. I want to interpolate/create a raster from those points that have temperature in order to have a continuous temperature surface. I have read through dozens of vignettes and used several different packages (but primarily
gstat). However, I'm still unsure about my methods. I think that the thin plate spline (
fields::Tps) is probably the best, but still don't quite understand the difference between
predictSurface. I'll illustrate this below. The results seem to be fairly similar between the methods, but just don't understand the difference between the two.
library(fields) library(dplyr) library(elevatr) library(sf) library(raster) # Create a bounding box for cropping elevation data. swcolorado <- c( "xmin" = -105, "ymin" = 36, "xmax" = -109.060253, "ymax" = 39 ) %>% sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_as_sf(crs = 4326) %>% sf::st_transform(crs = 4326) # Get elevation data, then crop with the bounding box, and aggregate (so the code runs quicker). elev.raster <- elevatr::get_elev_raster(swcolorado, z = 5) %>% raster::crop(elevation, swcolorado) %>% raster::aggregate(fact = 12) # Temperature dataframe with temperature, longitude, latitude, and elevation. temp.points <- structure(list(long = c(-108.160375, -107.808675, -106.510525, -107.81571, -105.514165, -106.45052, -105.740635, -105.63166, -106.843695, -107.61537, -107.682735, -107.6975, -108.503889, -105.724665, -108.10256), lat = c(37.473105, 37.647645, 36.047665, 37.60664, 37.569525, 37.021745, 37.71109, 38.08806, 37.611035, 37.902105, 37.74759, 37.737778, 37.473889, 37.67773, 37.46906 ), temperature = c(14.6972137463016, 14.8096938066668, 14.9004909660972, 14.8110845322457, 9.7248919765495, 13.271146269724, 18.0169186561757, 10.5835671771822, 10.8751447949415, 8.98234964596861, 11.8225704515042, 11.982697080223, 19.8293764051047, 18.0406584631011, 12.8451985711037 ), elev.locs = c(3024, 2694, 2932, 2706, 3589, 3054, 2293, 3449, 3471, 3739, 3205, 3329, 2111, 2290, 3291)), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 12L, 13L, 14L, 15L, 16L), class = "data.frame") # Need to extract the xy grid and put in ascending order, as the fields package expects that. elev.raster.long <- raster::xFromCol(elev.raster) elev.raster.lat <- raster::yFromRow(elev.raster) %>% sort() elev.raster.elev <- as.matrix(elev.raster) elev.raster.elev <- elev.raster.elev # Transpose, so that rows and columns will match the long lat lists. Then, mirror the columns so that the latitude is ascending. elev.raster.elev <- t(elev.raster.elev) %>% as.data.frame() elev.raster.elev <- elev.raster.elev[, order(ncol(elev.raster.elev):1)] %>% as.matrix() # Put long, lat, and elevation into 1 list. Then, rename to x, y, and z. elev.raster.list <- list(elev.raster.long, elev.raster.lat, elev.raster.elev) names(elev.raster.list) <- c("x", "y", "z") #Create the grid list, which will be used in the prediction below. grid.list <- list(x = elev.raster.list$x, y = elev.raster.list$y) # First, run the model using Tps with elevation as an independent covariate. obj <- fields::Tps( # Accepts points but expects them as matrix. x = as.matrix(temp.points[, c("long", "lat")]), # The dependent variable. Y = temp.points$temperature, # Elevation as an independent covariate. Z = temp.points$elev.locs, miles = TRUE ) # Use predictSurface on the model. out.p <- predictSurface(obj, grid.list = grid.list, ZGrid = elev.raster.list, extrap = TRUE) # Then, convert to raster. out.raster <- raster(out.p) crs(out.raster) <- CRS('+init=EPSG:4326') out.raster <- crop(out.raster, sf::st_bbox(c( "xmin" = -109, "ymin" = 36.0, "xmax" = -105, "ymax" = 38.25 ))) plot(out.raster)
Now using the same data, but using
interpolate() and adding elevation as another independent variable (rather than as a linear covariate).
tps2 <- Tps(as.matrix(temp.points[, c("long", "lat", "elev.locs")]), temp.points$temperature) p2 <- interpolate(elev.raster, tps2, xyOnly = FALSE) p2 <- crop(p2, sf::st_bbox(c( "xmin" = -109, "ymin" = 36.0, "xmax" = -105, "ymax" = 38.25 ))) plot(p2)
Next, I use
interpolate() again, but adding elevation as a linear covariate (rather than as an independent variable).
tps3 <- Tps(as.matrix(temp.points[, c("long", "lat")]), temp.points$temperature, Z = temp.points$elev.locs) p3 <- interpolate(elev.raster, tps3, xyOnly = FALSE, fun = pfun) p3 <- crop(p3, sf::st_bbox(c( "xmin" = -109, "ymin" = 36.0, "xmax" = -105, "ymax" = 38.25 ))) plot(p3)
Also, sorry for the additional complicated code in creating the grid list. If you know of a better way to do this quicker and cleaner, then that would be great too. I adapted the grid list conversion method from a vignette.
predictSurface portions of the code above take about 45 seconds each to run on my Mac.