I have a raster stack (or brick) that represent spatial locations as x-y and time as z. I have a spatial points shapefile (sf object) representing events for which I have the spatial localisation (also the x-y) and the date of the event (also z). I want to extract the specific raster value associated to the event day.

This is, I think, a relatively frequent problem that I manage to solve with the raster package, the sf package and a bunch of lapply (see my reproducible example below). However, after looking at my code, I got that feeling that there must be another, more official, way of doing this. That maybe it was already coded in the raster, sf, stars or another package in R.

Here is an example to show what I want to do and my code to do it:

Preparing the session


Making the raster

ras <- lapply(1:7, function(xx) runif(100, xx, xx+1)) %>% 
  lapply(matrix, nrow=10, ncol=10) %>% 
  lapply(raster, xmn=0, xmx=10, ymn=0, ymx=10) %>% 
  do.call("stack", .)

Making the point shapefile

  points <- data.frame(x=runif(10,0,10), y=runif(10,0,10), date=sample(1:7,10, replace = T)) %>% 
    st_as_sf(coords = c("x", "y"))

Building a extract function

extract_3d <- function(row, ras){
  row$date %>% 
    extract2(ras, .) %>% 

Running the function and saving the results

points$res <- points %>% 
  split(., 1:nrow(.)) %>% 
  sapply(extract_3d, ras=ras)
Simple feature collection with 10 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.392242 ymin: 0.2399111 xmax: 9.754538 ymax: 8.89191
epsg (SRID):    NA
proj4string:    NA
   date                   geometry      res
1     1  POINT (8.950553 8.346364) 1.579130
2     5    POINT (6.264524 3.0421) 5.793717
3     1   POINT (9.754538 8.89191) 1.033745
4     4  POINT (2.038216 6.965259) 4.299902
5     2  POINT (2.716971 3.044953) 2.633210
6     5 POINT (1.718511 0.2399111) 5.200416
7     2  POINT (3.215259 4.115244) 2.283648
8     6  POINT (0.392242 2.198866) 6.711109
9     7   POINT (7.32996 5.796152) 7.738792
10    7  POINT (7.956075 4.053148) 7.291688

My code is relatively simple, but it do not scale up to big data sets. I ran it on ~4 millions date on a tiled Daymet rasters and it's been running for over 3 days on 24 cores and 128 gb of ram. The idea here is really to do a space-time uneven extract between a raster brick and a spatial point object in R.

Is there a better, more official and efficient way of doing this that my extract_3d function called in a lapply.

1 Answer 1


Your code seems overly complex for a problem that can just be indexed. I would even say that it would be far less confusing if you just extracted all of the values in the stack and then operated on the resulting matrix. All you need are the dates that correspond to the raster stack and match the dates in the sf points data. I have performed operations like this on millions of points and long time-series daily climate data and have not seen performance penalties like you are detailing so, there is something in your tidy approach that is really slowing things down. There is likely some type of object coercion occurring in the background that is causing the issue (one of the disadvantages of pipe based workflows).

Here is code that follows a workflow of 1) extract the stack values for all points; 2) iterate through the matrix and match the date of the event to the associated row (point) and column (date) to assign the correct event value 3) assign the resulting vector back the point data.

Here we add libraries and create some data. I create a dates vector that is associated with the dates of the stack. I then randomly assign a date "event" to each point in the sf points data to represent the event.


r <- do.call(stack, replicate(n,raster(matrix(runif(i*j), i, j))))
( dates <- seq(as.Date("2000/1/1"), by = "month", length.out = n) )
s <- as(sampleRandom(r[[1]], 100, sp=TRUE), "sf")
  s$date <- sample(dates, nrow(s), replace=TRUE)

Now, I extract all of the data, resulting in each row representing a point and each column a raster (or date).

e <- extract(r, s)

Please note that if you are after an increase in speed, I would recommend moving your analysis to the terra package. The only difference would be that you would use the rast function in place of stack and you have to coerce the sf data to a vect class, which you can do on the fly. It would look something like e <- extract(r, vect(s))

Now that we have a matrix representing points as rows and dates as columns, we can iterate through the matrix and extract the correct event value based on matching an index value of the date to the correct column. This could easily be put into an lapply type function, potentially to leverage a multithreaded version like in the future.apply package. But, let's keep it in a for loop for transparency's sake.

event <- vector()
  for(i in 1:nrow(e)) {
    event[i] <- e[i,][which(dates %in% s$date[i])]
s$event <- event

Here is what is going on. Let's look at the 10th point observation, say i=10


Since we have a date field and a vector of dates that match the dates in the raster stack. All we have to do is match the two to find the index of the correct column in the extracted matrix.

Here is the date we are after s$date[10] from the point data, and the dates vector print(dates) that we want to match (representing the correct raster layer for the given event).

Here is what it looks like in practice, first the column index and then the extracted value.

which(dates %in% s$date[10])
e[i,][which(dates %in% s$date[10])]

As you can see, this is a simple matching exercise that can be performed on the extracted data as opposed to subsets of the raster stack. This approach would still work for points with multiple events because, you functionally have the data in a long format so, each new event would be picked up in the for loop.

You can also simply index a raster stack for extraction. Let's create point data with multiple events at a single location.

s <- as(sampleRandom(r[[1]], 100, sp=TRUE), "sf")
  s$ID <- 1:nrow(s)
    s <- rbind(s, s[sample(1:nrow(s), 100, replace=TRUE),])
  s$date <- sample(dates, nrow(s), replace=TRUE)

Now we can subset it and extract the event values using a raster index for just those dates. I will leave it to you on how to omit the diagonal value.

( s.sub <- s[s$ID == s$ID[which(duplicated(s$ID))[1]],] )
extract(r[[which(dates %in% s.sub$date)]], s.sub) 
  • Wow, there is a lot of meat in this answer! I'll test all of that this week and accept the answer if it works. Thanks!
    – Bastien
    Commented Nov 17, 2020 at 11:55
  • I've tested your different solutions. First, terra is amazing, thanks for the link, it is fast. Now, interestingly, your approach do not always beats mine when using the raster package. With my real data and 100 entry point, I'm about 40% faster. However, as the number of data point increase, your's is faster. However, Ram usage increase as well. On 300000, with terra on a raster stack of 7665 layers, your approach take more that 120gb. So more work is needed. Your answer helped me think of other solutions. thanks.
    – Bastien
    Commented Nov 23, 2020 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.