How to extract values from raster stack to point SFC based on attribute R?

Question

I would like to extract values from a raster stack to a number of sf point files. The raster stack I'm working with represents daily weather data, and the sf point files represent wildfire footprints, with an attribute corresponding to the day that the point burned. I'd like to extract weather data to these point files, specifically extracting the weather data for the day that the point burned.

I've made some progress in renaming the raster stack layer names to julian day, and subset the raster stack to the day of interest. As it turns out, I need to convert the sf object to a spatial object to get raster::extract working. The really tricky part is to subset the raster stack based on the day of burning for each point.

I can sort of see how I could iterate through the sf point object and rbind the results, but is this colossal loop really the only way to get this working?

library(raster)
library(sf)
library(rgdal)
library(lubridate)

Next I stack the rasters, which are Arc Grids.

fwi.list <- list.files(path = "C:/Example/", pattern="fwi")
fwi.stack <- stack(fwi.list)
crs(fwi.stack) <- "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

print(fwi.stack)
    class      : RasterStack 
dimensions : 1600, 1900, 3040000, 153  (nrow, ncol, ncell, nlayers)
resolution : 3000, 3000  (x, y)
extent     : -2600000, 3100000, -885076, 3914924  (xmin, xmax, ymin, ymax)
crs        : +proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
names      : X121, X122, X123, X124, X125, X126, X127, X128, X129, X130, X131, X132, X133, X134, X135, ... 
min values :    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0, ... 
max values : 1138,  681,  690,  662,  873,  618,  417,  893,  440,  511,  805,  522,  575,  543,  540, ... 

Here is the sf pointfile. The trick will be to assign fwi to each point based on the jday attre

GRID.PT <- st_read("C:/Example/470_2015.shp")
GRID.PT

Simple feature collection with 2126 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: -1039086 ymin: 2078687 xmax: -1015336 ymax: 2095187
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
First 10 features:
                   geometry DOB
1  POINT (-1026836 2095187) 183
2  POINT (-1026336 2095187) 183
3  POINT (-1026086 2095187) 183
4  POINT (-1025836 2095187) 183
5  POINT (-1027336 2094937) 183
6  POINT (-1027086 2094937) 183
7  POINT (-1026836 2094937) 183
8  POINT (-1026586 2094937) 183
9  POINT (-1026086 2094937) 183
10 POINT (-1025836 2094937) 183

So, how should I go about this? I'm sure it involves raster::subset and raster::extract, maybe something like: extract(subset(fwi.stack, paste0("X", GRID.PT$DOB[x])), as(GRID.PT, "Spatial"))

But should I write this as a function and use lapply? Or should I use a big loop? Thanks for your help, SO.

Answer 1

I've created a simple replicable example for this, hopefully it will give some insight into how to solve the issue you're facing.

## Set up the raster
test.r <- raster( matrix( nrow=5,
                          ncol=5,
                          data = rep(1,25) 
                          ) 
                  )
extent(test.r) <- c(0,5,0,5)
test.r[] <- 5

## Set up a stack
test.r <- stack(test.r,
                test.r+1,
                test.r+2,
                test.r+3,
                test.r+4)
## Name them in a similar fashion
names(test.r) <- paste0("X",1:5)

## Generate some point data
pts <- SpatialPoints(coords = matrix(ncol=2,
                                     c(rep(seq(0.5,4.5,1),6)
                                       )
                                     )
                     )

## Make it 'sf' for applicability
pts <- as(pts,"sf")

pts$val <- rep(1:5,3)

## Perform lapply that assigns values within
lapply( unique( pts$val ), function(x){

  pt <- pts[ pts$val == x, ]
  rast <- test.r[[ grep( x, names( test.r ) ) ] ]

  pts[ pts$val == x, "rast_val" ] <<- extract( rast, pt )

})

## Or in 1 line

lapply( unique( pts$val ), function(x){

  pts[ pts$val == x, "rast_val" ] <<- extract( test.r[[ grep( x, names( test.r ) ) ] ], 
                                               pts[ pts$val == x, ])

})