How to do a temporal gap filling using spline interpolation for a NDVI RasterStack?

Question

I have a spline interpolation function that I used for interpolating a RasterStack of 31 NDVI layers for the month of June 2018. The function needs to include the requirement of only generating pixel values inside a NDVI interval (-1 to +1).

My problem is that the RasterBrick resulted have pixel values outside the (-1 to +1) interval. How can I fix that?

So, to give a background of my goal:

I want to create a NDVI time series for June of 2018 with no time gaps. I have 5 NDVI raster images for the month:

• 2018-05-09
• 2018-05-14
• 2018-05-19
• 2018-05-24
• 2018-05-29

And the rest of the days I want to create empty rasters with NA pixel values to have then 31 raster images to stack them and then apply the spline function. Here are my steps:

STEP 1

I created an empty raster with NA values for all the missing days of the month, here is an example of one date:

#creating empty raster for the missing dates JUNE2018
#define raster layer
NDVI_20180501 <- raster(ncol=7925, nrow=7218, xmn=602770, xmx=682020, ymn=6590220, ymx=6662400, crs=CRS("+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"))
#set all pixel values to NA
NDVI_20180501[] <- NA
#save the raster (write to disk) as .tif file
writeRaster(NDVI_20180501, "NDVI_20180501.tif", format="GTiff", overwrite=TRUE)

STEP 2

Then I stack all the empty raster with the NDVI rasters:

#list of raster files
file_list <- list.files(pattern = "\\.tif$")

#Stack all rasters
rstack_ndvi_062018 <- stack(file_list)

STEP 3

Then thats the part that I believe lies the problem, I wrote the spline function inspired by this answer https://stackoverflow.com/a/37846446/21455915 but altering so the spline interpolation would generate values that are inside the NDVI interval (from -1 to +1):

f2 <- function(x) {
  z <- which(is.na(x))
  nz <- length(z)
  nx <- length(x)
  if (nz > 0 & nz < nx) { 
    # Find the min and max values within the -1 to +1 range
    min_val <- min(x[x >= -1 & x <= 1])
    max_val <- max(x[x >= -1 & x <= 1])
    
    # If there are no values within the range, set the range to -1 to +1
    if (is.na(min_val) | is.na(max_val)) {
      min_val <- -1
      max_val <- 1
    }
    
    # Interpolate the missing values using the spline function
    x[z] <- spline(x=1:nx, y=pmin(pmax(x, min_val), max_val), xout=z, method="natural")$y
  }
  x
}

STEP 4

Used calc() to apply the above spline function into the RasterStack created on STEP 2:

nogap_rstack_ndvi_062018_f2 <- calc(rstack_ndvi_062018, f2)

But when observing the result by:

nogap_rstack_ndvi_062018_f2

You can see that the pixel values for the interpolated layers are outside the -1 to +1 interval:

class : RasterBrick dimensions : 7218, 7925, 57202650, 31 (nrow, ncol, ncell, nlayers) resolution : 10, 10 (x, y) extent : 602770, 682020, 6590220, 6662400 (xmin, xmax, ymin, ymax) crs : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs source : r_tmp_2023-03-23_170523_8596_34248.grd names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... min values : -3.8575468, -3.5863018, -3.3150567, -3.0438116, -2.7725666, -2.5148830, -2.3982520, -2.2816209, -2.1649899, -2.0483589, -1.9317279, -1.8150969, -1.6984659, -1.5818349, -1.4652039, ... max values : 4.568034, 4.293518, 4.019002, 3.744486, 3.469970, 3.195453, 2.920937, 2.646421, 2.371905, 2.097389, 1.828929, 1.733057, 1.637185, 1.541313, 1.445441 , ...

How can I fix this?

And thanks for helping!