Я пытаюсь извлечь значения из растрового стека и добавить их в существующий фрейм данных. Значения представляют собой набор метрик (PatchStat из пакета r SDMtools), которые я могу извлечь в формат списка, но я застрял, пытаясь привязать значения к существующему фрейму данных.
Входные данные:
library(sp)
library(sf)
library(raster)
library(dplyr)
library(SDMTools)
mydata <- read.table(header=TRUE, text = "
animal X Y ord.year
1 pb_20414 157978.9 2323819 2009168
2 pb_20414 156476.3 2325586 2009168
3 pb_06817 188512.0 2299679 2006263
4 pb_06817 207270.9 2287248 2006264")
# add rasters
s <- stack(system.file("external/rlogo.grd", package="raster"))
names(s) <- c('masie_ice_r00_v01_2009168_4km', 'masie_ice_r00_v01_2006263_4km', 'masie_ice_r00_v01_2006264_4km')
# Create sp object
projection <-CRS('+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m + datum=WGS84 +no_defs +towgs84=0,0,0') # matches MASIE raster
coords <- cbind(mydata$X, mydata$Y)
mydata.sp <- SpatialPointsDataFrame(coords = coords, data = mydata, proj4string = projection)
# Create sf object
mydata.sf <- st_as_sf(mydata)
mydata.buf30 <- st_buffer(mydata.sf, 30000)
Моя цель - сопоставить каждую точку GPS (X, Y) с правильным GeoTIFF по дате (mydata$ord.year), обрезать растр в (пространственно явный) 30-километровый буфер, запустить PatchStat в программе SDMtools для R и добавить результаты в исходный фрейм данных. Суть в том, что результаты PatchStat предоставляются в кадре данных, поэтому у меня возникают проблемы с сопоставлением этих результатов с моим существующим кадром данных.
Вот пример результатов, представленных при запуске PatchStat:
patchID n.cell n.core.cell n.edges.perimeter n.edges.internal area core.area perimeter
2 3 73 13 86 206 73 13 86
perim.area.ratio shape.index frac.dim.index core.area.index
2 1.178082 2.388889 1.430175 0.1780822
Вот что я смог сделать до сих пор:
# separate date component of TIF name to correspond to mydata$ord.year
stack <- list()
date<-vector()
for (i in 1:length(rasterlist)) {
stack[[i]]<-raster(rasterlist[i])
tt<-unlist(strsplit(names(stack[[i]]), "[_]"))
date[i]<-tt[which(nchar(tt)==max(nchar(tt)))]
}
st <- stack(stack) # Create rasterstack object
# crop raster to buffer
mydata.sp <- as(mydata.sf, 'Spatial') # back to sp object
# pull raster data from GeoTIFF that corresponds to ordinal date
pat <- list()
for (i in 1:nrow(mydata.sp)) {
st2<-st[[which(date==mydata.sp$ord.year[i])]]
GeoCrop <- raster::crop(st2, mydata.sp[i,])
GeoCrop_mask <- raster::mask(GeoCrop, mydata.sp[i,])
pat[[i]] <- PatchStat(GeoCrop_mask)}
Кроме того, я исключил один из двух типов растительного покрова, чтобы у каждого элемента в списке была только одна строка:
pat2 <- lapply(pat, `[`, -1,) # remove first row in each list element so only one row remains (using program plyr for R)
Теперь я хотел бы сопоставить эти строки с моим исходным фреймом данных, чтобы pat2[[1]] добавлялся к mydata.sp[1,] следующим образом (предполагая, что a, b и c - столбцы метаданных в моем исходном SpatialPointsDataFrame). Я хотел бы, чтобы все столбцы данных из PatchStat были добавлены, но для экономии времени и места я включил только первые три здесь:
a b c PatchID n.cell n.core.cell
1 2 3 3 73 13
Примечание: Если возможно, я бы хотел, чтобы весь этот процесс был включен в цикл for, чтобы минимизировать пространство для ошибок, а также время обработки.
Большое спасибо!