#!/usr/bin/env Rscript

library(raster)

# Create gridded example data
xyz <- expand.grid(x = seq(-112.4, -111.6, by = 0.03),
                   y = seq(40.3, 40.9, by = 0.03))
xyz$z <- with(xyz, abs(x * y))
r_coarse <- rasterFromXYZ(xyz, crs = '+proj=longlat')
r_coarse
# class       : RasterLayer
# dimensions  : 21, 27, 567  (nrow, ncol, ncell)
# resolution  : 0.03, 0.03  (x, y)
# extent      : -112.415, -111.605, 40.285, 40.915  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84
# data source : in memory
# names       : z
# values      : 4498.286, 4597.16  (min, max)

xres <- 0.01
yres <- 0.01
r_interpolated <- projectRaster(r_coarse,
                                res = c(xres, yres),
                                crs = '+proj=longlat',
                                method = 'bilinear')
r_interpolated
# class       : RasterLayer
# dimensions  : 73, 91, 6643  (nrow, ncol, ncell)
# resolution  : 0.01, 0.01  (x, y)
# extent      : -112.465, -111.555, 40.235, 40.965  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84
# data source : in memory
# names       : z
# values      : 4496.767, 4598.693  (min, max)

op <- par(mfrow = c(1, 2))
plot(r_coarse, main = 'Coarse', col = viridis::viridis(256))
plot(r_interpolated, main = 'Interpolated', col = viridis::viridis(256))
par(op)