Matching coordinates to the closest raster in a stack by date in R

Question

I have downloaded NDVI rasters from luna::getModis, and assigned each raster layer the appropriate date. This was then converted to a rasterstack. I also have a dataframe with coordinates and date.

I want to extract only the raster that is closest in date to the each coordinate and add that value to the dataframe.

My current code looks like:

#Donwload rasters for study area and time
product <- "MYD13Q1" #Create product
start <- "2020-11-01" 
end <- "2020-12-01"
aoi <- c(-15.90077, -13.66578, 14.53391, 16.75448) #Coordinates of area - can use shp file but wasn't working so found this easier
user <- "X" #Can use mine or someone elses 
pwd <- "Y"
path_ndvi <- file.path("~", "" ,"" fsep="/") #Where to save files 
setwd(path_ndvi)

mf_download <- luna::getModis(product, start, end, aoi, download=TRUE,
                          path=path_ndvi, username=user, password=pwd) #Download data 

#Create stack of only NDVI files 

mf_all = paste0(mf[grep("MYD13Q1", mf)]) #Create list of all files 
ndvi = lapply(mf_all, function(x) rast(x)[[1]]) #Extract files that are only layer 1 = NDVI
ndvi_stack = c(ndvi) #Change this into a list of files (terra)
ndvi_stack = rast(ndvi_stack) #Change the ndvis into a stack 
plot(ndvi_stack)

#Reproject
ndvi_stack <- project(ndvi_stack, "EPSG:32628")

#Add dates to stack
( dates <- seq(as.Date("2020/10/23"), by = 16, length.out = 3) ) #Create list of dates for rasters 
time(ndvi_stack_qc_repro) <- dates #Add dates to rasters 
terra::time(ndvi_stack_qc_repro) #Check dates added

But after this, I am stuck on how to extract by corresponding dates.

I look forward to your answers.

Answer 1

First of all, the code you provided does not run seamlessly on someone elses machine. Please make sure that all variables are properly defined so that it is less work for others to reproduce your problem. Secondly, you do not give an example of your data.frame which you say contains coordinates and dates. So, in the example below I created one myself, which is based on a random sample of 10 points from non-NA SpatRaster cells.

Once you have the SpatRaster and a data.frame with coordinates and dates, you can convert to a SpatVector, extract values for all layers and then find the column index closest to the specific date. The output below will also include the originally extracted coloumns so you can make sure that the code does what is expected:

# prepare MODIS data
library(terra)
#> terra 1.7.46

outdir <- "./modis-data"
epsg <- "EPSG:32628"
(files <- list.files(outdir, full.names = TRUE))
#> [1] "./modis-data/MYD13Q1.A2020297.h16v07.061.2020353025751.hdf"
#> [2] "./modis-data/MYD13Q1.A2020313.h16v07.061.2020353215522.hdf"
#> [3] "./modis-data/MYD13Q1.A2020329.h16v07.061.2020363095444.hdf"
dates <- sapply(files, function(x) strsplit(basename(x), "\\.")[[1]][2], USE.NAMES = FALSE)
dates <- as.Date(dates, "A%Y%j") # translate dates using DOY symbol
ndvi_stack <- rast(lapply(files, function(x) rast(x, lyrs = 1)))
time(ndvi_stack) <- dates
ndvi_stack <- project(ndvi_stack, epsg)
names(ndvi_stack) <- paste0("NDVI-", dates)
ndvi_stack
#> class       : SpatRaster 
#> dimensions  : 4657, 4954, 3  (nrow, ncol, nlyr)
#> resolution  : 240.1828, 240.1828  (x, y)
#> extent      : -158310.6, 1031555, 1105331, 2223862  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 28N (EPSG:32628) 
#> source(s)   : memory
#> names       : NDVI-2020-10-23, NDVI-2020-11-08, NDVI-2020-11-24 
#> min values  :       -20000000,       -20000000,       -20000000 
#> max values  :        99740000,        99628584,        98764727 
#> time (days) : 2020-10-23 to 2020-11-24

# prepare a sample data.frame
set.seed(152)
sample <- spatSample(ndvi_stack, size = 10, xy = TRUE, na.rm = TRUE, values = FALSE)
sample <- as.data.frame(sample)
sample$date <- time(ndvi_stack)[1]  + runif(10) * (time(ndvi_stack)[3] - time(ndvi_stack)[1])
sample
#>           x       y       date
#> 1  320253.5 1590380 2020-11-13
#> 2  930798.1 1953777 2020-11-11
#> 3  610634.5 1178467 2020-10-30
#> 4  771316.8 1974192 2020-10-27
#> 5  503753.2 1469088 2020-11-20
#> 6  498949.5 2134875 2020-10-26
#> 7  682209.0 1823117 2020-10-25
#> 8  445869.1 1627128 2020-11-08
#> 9  398312.9 1862507 2020-11-12
#> 10 958659.3 2125267 2020-11-08

# first, convert to SpatVector, then do the extraction
sample_vect <- vect(sample, geom = c("x", "y"), crs = epsg)
sample_ext <- extract(ndvi_stack, sample_vect, bind = TRUE)
(df_ext <- values(sample_ext))
#>          date NDVI-2020-10-23 NDVI-2020-11-08 NDVI-2020-11-24
#> 1  2020-11-13        39727939        35156511        31914680
#> 2  2020-11-11        17648983        15990054        16292090
#> 3  2020-10-30        57035264        58819043        53931318
#> 4  2020-10-27        15820405        15092227        15943798
#> 5  2020-11-20        75032705        64615513        50219214
#> 6  2020-10-26        12648618        11888546        12277130
#> 7  2020-10-25        17971768        18015492        18374478
#> 8  2020-11-08        41420464        33183057        30414272
#> 9  2020-11-12        26822114        23814448        23446323
#> 10 2020-11-08        11084589        11183971        11237823

# now, we pre-calculate the index with the smallest distance and extract the values
index_date <- sapply(df_ext$date, function(date) which.min(abs(dates - date)) + 1) # + 1 for column index
closest_ndvi <- sapply(1:nrow(df_ext), function(i) df_ext[i , index_date[i]])
closest_date <- as.Date(names(df_ext)[index_date], "NDVI-%Y-%m-%d")
ndvi <- data.frame(closest_date = closest_date, closest_ndvi = closest_ndvi)
# bind it all together for inspection
(ndvi <- cbind(df_ext, ndvi))
#>          date NDVI-2020-10-23 NDVI-2020-11-08 NDVI-2020-11-24 closest_date
#> 1  2020-11-13        39727939        35156511        31914680   2020-11-08
#> 2  2020-11-11        17648983        15990054        16292090   2020-11-08
#> 3  2020-10-30        57035264        58819043        53931318   2020-10-23
#> 4  2020-10-27        15820405        15092227        15943798   2020-10-23
#> 5  2020-11-20        75032705        64615513        50219214   2020-11-24
#> 6  2020-10-26        12648618        11888546        12277130   2020-10-23
#> 7  2020-10-25        17971768        18015492        18374478   2020-10-23
#> 8  2020-11-08        41420464        33183057        30414272   2020-11-08
#> 9  2020-11-12        26822114        23814448        23446323   2020-11-08
#> 10 2020-11-08        11084589        11183971        11237823   2020-11-08
#>    closest_ndvi
#> 1      35156511
#> 2      15990054
#> 3      57035264
#> 4      15820405
#> 5      50219214
#> 6      12648618
#> 7      17971768
#> 8      33183057
#> 9      23814448
#> 10     11183971

^{Created on 2023-10-19 with reprex v2.0.2}