如何根据属性R从栅格堆栈中提取值到点SFC？

Question

我想将栅格堆栈中的值提取到多个 sf 点文件中。我正在使用的栅格堆栈代表每日天气数据，而 sf 点文件代表野火足迹，其属性对应于该点燃烧的日期。我想将天气数据提取到这些点文件中，具体提取点燃烧当天的天气数据。

我在将栅格堆栈图层名称重命名为儒略日并将栅格堆栈子集化为感兴趣的日期方面取得了一些进展。事实证明，我需要将 sf 对象转换为空间对象才能使 raster::extract 正常工作。真正棘手的部分是根据 每个点 .

的燃烧日期对栅格堆栈进行子集化

我大概明白了如何遍历 sf 点对象和 rbind 结果，但是这个庞大的循环真的是实现此功能的唯一方法吗？

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

接下来我堆叠栅格，它们是 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, ...

这是 sf 点文件。诀窍是根据 jday attre

将 fwi 分配给每个点

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

那么，我应该怎么做呢？我确定它涉及 raster::subset 和 raster::extract，可能是这样的： extract(subset(fwi.stack, paste0("X", GRID.PT$DOB[x])), as(GRID.PT, "Spatial"))

但是我应该把它写成一个函数并使用lapply吗？还是我应该使用大循环？谢谢你的帮助，SO。

Answer 1

我为此创建了一个简单的可复制示例，希望它能让您深入了解如何解决您面临的问题。

## 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, ])

})

如何根据属性R从栅格堆栈中提取值到点SFC？

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

gis

r

raster

sf