0

I have many point locations with time-date stamps and 18 raster layers each containing the annual NDVI values from the year 2002 till 2019. I want to extract the values for the locations from the raster layer of the year similar to the year of the location time stamp.

Example:

location1 time stamp 12.04.2003 --> extract value from raster layer 2003
location2 time stamp 12.06.2012 --> extract value from raster layer 2012

This is a sample of my location data (EPSG 3035):

  X id areas sex      x_      y_ case                time
1 1  1    13   f 4141879 2762606    1 2010-11-16 00:59:00
2 2  1    13   f 4142395 2759956    1 2010-11-16 21:59:00
3 3  1    13   f 4143634 2761615    1 2010-11-17 04:59:00
4 4  1    13   f 4143171 2761593    1 2010-11-17 11:59:00
5 5  1    13   f 4144488 2762547    1 2010-11-17 18:59:00
6 6  1    13   f 4143885 2761944    1 2010-11-18 08:59:00

First I created a raster stack with all 18 raster layers and added the date for each layer:

#create raster stack NDVI 
setwd("E:/MA_2Try/annual_NDVI/average")
grids <- list.files("E:/MA_2Try/annual_NDVI/average" , pattern = "*.tif$")
NDVI_stack=stack(grids)

#Add years to stack
dt<-as.data.frame(
  as.POSIXct(c('2002-01-01','2003-01-01','2004-01-01','2005-01-01','2006-01-01',
               '2007-01-01','2008-01-01','2009-01-01','2010-01-01','2011-01-01',
               '2012-01-01','2013-01-01','2014-01-01','2015-01-01','2016-01-01',
               '2017-01-01','2018-01-01','2019-01-01'))) 

NDVI_stack <- setZ(NDVI_stack, dt[,1], "SampleDate")

Here is the resulting info for raster [[1]]

> NDVI_stack[[1]]
class      : RasterLayer 
dimensions : 4125, 5503, 22699875  (nrow, ncol, ncell)
resolution : 250, 250  (x, y)
extent     : 3530500, 4906250, 2229250, 3260500  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source     : ndvi_average_2002.tif.tif 
names      : ndvi_average_2002.tif 
SampleDate : 2002-01-01

How can I achieve this?

Nicole
  • 1
  • 1

2 Answers2

0

You can try something like this:

library(raster)
library(sf)
rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
s <- stack(rast,rast/100,rast*100)
s <- setZ(s, c('2010-01-01','2011-01-01','2013-01-01'), "SampleDate")

pt1 = st_point(c(0,1))
pt2 = st_point(c(1,1))
d = data.frame(id = 1:2,sex=c("m","f"),date = c("2010-11-16 00:59:00","2011-11-16 21:59:00"))
d$geom = st_sfc(pt1, pt2)
df = st_as_sf(d)


myex <- function(rst,point){
  ex <- as.list(1:nrow(point))
  for (i in 1:nrow(point)) {
    year_p <- as.integer(lubridate::year(point[i,]$date))
    diff <- abs(year_p-as.integer(lubridate::year(unlist(rst@z))))
    index <- which(diff==min(diff))
    ex[[i]] <- raster::extract(rst[[index]],point[i,],include_cols=colnames(point))
  }
  ex <- as.data.frame(do.call(rbind,ex))
  ex$id <- point$id
  return(ex)
}

myex(s,df)

where the output looks like this:

    V1 id
1 10.0  1
2  0.2  2

The idea is to choose the layer with the minimum difference between the point's date and the layer's date and use it for the extraction. Here I use a for loop, but it is possible to use lapply or one of its parallel versions (such as snowfall::sfLaplly). Moreover, if you create a small buffer around points, you can substitute raster::extract with exactextract::exact_extract, which is way faster

Elia
  • 2,210
  • 1
  • 6
  • 18
0

Here is how you can do that with terra (but it should be easy to adjust for raster)

Example data

library(terra)
s <- rast(system.file("ex/logo.tif", package="terra"))   
names(s) <- 2010:2012
pts <- matrix(c(14.53, 70.41, 12.93, 53.8, 60.21, 47.41), ncol=2, byrow=TRUE)
time <- as.POSIXct(c("2011-11-16 00:59:00", "2010-01-16 00:59:00", "2010-11-16 00:59:00"))

Get the layers (years) that you want from the time stamps and match them to the layers. A simple approach would be to round up the year if the month > June:

y <- format(time, format='%Y') |> as.integer()
m <- format(time, format='%m') |> as.integer()
y <- y + (m > 6)
y
#[1] 2012 2010 2011

Now use extract

e <- extract(s, vect(pts), layer=as.character(y))
e
#  ID layer value
#1  1  2012   253
#2  2  2010    73
#3  3  2011    14

With raster you could this last step like this:

ee <- extract(s, pts)
ee
#  2010 2011 2012
#1  255  255  253
#2   73   74   60
#3    6   14   27

j <- match(y, names(s))
j
#[1] 3 1 2
     
ee[cbind(1:nrow(ee), j)]
#[1] 253  73  14

But note that raster may pre-pend "X"s before the layernames of they start with a character between 0 and 9.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63