# How to use a temperature raster (e.g., PRISM) to constrain temperature values in a thin plate spline regression and interpolation

I have point data with temperature, latitude, longitude, and elevation. I am interpolating across space to the extent of those points, and have been using elevation as a covariate in the fields::Tps. However, this is for paleoclimate reconstructions.

So, I know what the modern day temperature surface looks like (i.e., the distribution of temperature across space). It seems reasonable to use this temperature surface rather than elevation in some way in order to constrain the interpolation. In other words, as interpolation takes place between a point (with a known temperature) and another random point in space, is it possible to say if between these two points is extremely different than the modern surface, then it is assigned a cutoff value (almost like a threshold)? I would like to do this in a way that would be quantitative, and not me just choosing a threshold. This would almost be like having a calibrated dataset and then predictions from the point data.

I had also been thinking about using some kind of similarity/dissimilarity measure, but unsure how to do that in this context. My other thought was to use the maximum amount of change between two pixels from the modern temperature dataset (PRISM) to constrain the interpolation. I attempt to illustrate this below:

In the illustration above, the maximum change in the modern temperature dataset between 2 pixels is 3 (shown by the 2 way arrow). Therefore, there should not be a change in the interpolation greater than or equal to 3 (using absolute value) between any 2 connected pixels in the interpolation. So, you can see how the cells would change with this constraint.

Here, I provide some sample code for what I have done, and also include the PRISM data.

library(fields)
library(dplyr)
library(elevatr)
library(sf)
library(raster)

# Create a bounding box for cropping elevation data.
"xmin" = -109,
"ymin" = 36,
"xmax" = -105,
"ymax" = 38.25
)  %>%
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::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)


To get PRISM data, for the area above:

library(prism)

# Need to first create a directory for the PRISM data to download and set as working directory.
setwd("~/prism_set_dl_dir")

# Then download data (only doing one layer so that it isn't a large file).
get_prism_annual("tmean", years = 2010, keepZip = FALSE, keep_pre81_months = FALSE)

# Load the PRISM data and I crop to the same extent as above.
PRISM <- raster("~/prism_set_dl_dir/PRISM_tmean_stable_4kmM3_2010_bil/PRISM_tmean_stable_4kmM3_2010_bil.bil") %>%
crop(sf::st_bbox(c(
"xmin" = -109,
"ymin" = 36.0,
"xmax" = -105,
"ymax" = 38.25
)))