Как извлечь значения из растрового стека в точку SF C на основе атрибута R? - PullRequest
Новая
       43

Как извлечь значения из растрового стека в точку SF C на основе атрибута R?

0 голосов
/ 16 апреля 2020

Я хотел бы извлечь значения из растрового стека в несколько файлов sf-точек. Растровый стек, с которым я работаю, представляет ежедневные данные о погоде, а файлы sf-точек представляют собой следы лесного пожара с атрибутом, соответствующим дню, когда точка сгорела. Я хотел бы извлечь данные о погоде в эти файлы точек, в частности извлекать данные о погоде за день, когда точка была сожжена.

Я добился определенного прогресса в переименовании имен слоев растрового стека в юлианский день, и подмножество растрового стека в день интереса. Как оказалось, мне нужно преобразовать объект sf в пространственный объект, чтобы заставить raster::extract работать. Действительно сложная часть состоит в том, чтобы поднастроить растровый стек на основе дня записи для каждой точки .

. Я могу как-то увидеть, как я мог бы выполнить итерацию через sf точечный объект и rbind результатов, но действительно ли этот колоссальный l oop действительно единственный способ заставить это работать? 

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

Далее я складываю растры, которые являются Ar c 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. Хитрость заключается в том, чтобы назначить fwi каждой точке на основе 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

Итак, как мне go об этом? Я уверен, что это включает raster::subset и raster::extract, может быть что-то вроде: extract(subset(fwi.stack, paste0("X", GRID.PT$DOB[x])), as(GRID.PT, "Spatial"))

Но я должен написать это как функцию и использовать lapply? Или я должен использовать большой l oop? Спасибо за вашу помощь, ТАК.

1 Ответ

1 голос
/ 16 апреля 2020

Я создал для этого простой воспроизводимый пример, надеюсь, он поможет понять, как решить проблему, с которой вы столкнулись. 

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

})
...