从栅格中提取值并附加到现有数据框
Extract values from raster and append to existing dataframe
我正在尝试从 rasterstack 中提取值并将其附加到现有数据框。这些值是一组指标（来自 r 包 SDMtools 的 PatchStat），我可以将其提取为列表格式，但我无法尝试将这些值绑定到我现有的数据框。
输入数据：
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)
我的目标是按日期 (mydata$ord.year) 将每个 GPS 点 (X,Y) 与正确的 GeoTIFF 相匹配，将栅格裁剪到（空间明确的）30 公里缓冲区，运行 R 程序 SDMtools 中的 PatchStat，并将结果附加到原始数据帧。问题是 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 的所有数据列，但为了节省时间和 space，我只在此处包括前三个：
a b c PatchID n.cell n.core.cell
1 2 3 3 73 13
注意：如果可能的话，我希望将整个过程包含在 for 循环中，以最大限度地减少错误空间和处理时间。
非常感谢！
好吧，我做了这件非常丑陋的事情并得到了我想要的。但我不喜欢它。如果有人有更好的主意，我很想听听！
# Change objects to df
pat2 <- lapply(pat, `[`, -1,) # remove first row in each list element
library(plyr) # ldply command
pat3 <- ldply (pat2, data.frame)
pat4 <- bind_cols(pb, pat3)
感谢您努力提供示例数据。但是还是不完整（指的是我们没有的文件，你可以这样
library(raster)
library(SDMTools)
s <- stack(system.file("external/rlogo.grd", package="raster"))
s <- round(s / 50) # to have fewer patches
names(s) <- c('masie_ice_r00_v01_2009168_4km', 'masie_ice_r00_v01_2006263_4km', 'masie_ice_r00_v01_2006264_4km')
df <- data.frame(ord.year=c("2009168", "2009168", "2006263", "2006264"))
pts <- SpatialPoints(cbind(c(20,40,60,80), c(20,40,60,20)))
crs(pts) <- crs(s)
pts <- SpatialPointsDataFrame(pts, df)
制作缓冲区
b <- buffer(pts, 15, dissolve=FALSE)
获取匹配的名称
nms <- names(s)
nms <- gsub('masie_ice_r00_v01_', '', nms)
nms <- gsub('_4km', '', nms)
循环匹配姓名，并将结果放入列表
p <- list()
for (i in 1:length(b)) {
j <- which(b$ord.year[i] == nms)
r <- s[[j]]
z <- crop(r, b[i,])
z <- mask(z, b[i,])
p[[i]] <- PatchStat(z)
}
请注意，p 的每个元素都有一个包含多行和多列的 data.frame。
p[[1]]
#patchID n.cell n.core.cell n.edges.perimeter n.edges.internal area core.area perimeter perim.area.ratio shape.index frac.dim.index core.area.index
#1 1 53 5 68 144 53 5 68 1.2830189 2.266667 1.427207 0.09433962
#2 2 123 8 182 310 123 8 182 1.4796748 3.956522 1.586686 0.06504065
#3 3 149 31 190 406 149 31 190 1.2751678 3.800000 1.543074 0.20805369
#4 4 54 2 114 102 54 2 114 2.1111111 3.800000 1.679578 0.03703704
#5 5 337 206 146 1202 337 206 146 0.4332344 1.972973 1.236172 0.61127596
如果您只想要第一行
pp <- t(sapply(p, function(i) i[1,]))
将它与原来的 data.frame 结合起来现在变得微不足道
dfpp <- cbind(df, pp)
我正在尝试从 rasterstack 中提取值并将其附加到现有数据框。这些值是一组指标（来自 r 包 SDMtools 的 PatchStat），我可以将其提取为列表格式，但我无法尝试将这些值绑定到我现有的数据框。
输入数据：
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)
我的目标是按日期 (mydata$ord.year) 将每个 GPS 点 (X,Y) 与正确的 GeoTIFF 相匹配，将栅格裁剪到（空间明确的）30 公里缓冲区，运行 R 程序 SDMtools 中的 PatchStat，并将结果附加到原始数据帧。问题是 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 的所有数据列，但为了节省时间和 space，我只在此处包括前三个：
a b c PatchID n.cell n.core.cell
1 2 3 3 73 13
注意：如果可能的话，我希望将整个过程包含在 for 循环中，以最大限度地减少错误空间和处理时间。
非常感谢！
好吧，我做了这件非常丑陋的事情并得到了我想要的。但我不喜欢它。如果有人有更好的主意，我很想听听！
# Change objects to df
pat2 <- lapply(pat, `[`, -1,) # remove first row in each list element
library(plyr) # ldply command
pat3 <- ldply (pat2, data.frame)
pat4 <- bind_cols(pb, pat3)
感谢您努力提供示例数据。但是还是不完整（指的是我们没有的文件，你可以这样
library(raster)
library(SDMTools)
s <- stack(system.file("external/rlogo.grd", package="raster"))
s <- round(s / 50) # to have fewer patches
names(s) <- c('masie_ice_r00_v01_2009168_4km', 'masie_ice_r00_v01_2006263_4km', 'masie_ice_r00_v01_2006264_4km')
df <- data.frame(ord.year=c("2009168", "2009168", "2006263", "2006264"))
pts <- SpatialPoints(cbind(c(20,40,60,80), c(20,40,60,20)))
crs(pts) <- crs(s)
pts <- SpatialPointsDataFrame(pts, df)
制作缓冲区
b <- buffer(pts, 15, dissolve=FALSE)
获取匹配的名称
nms <- names(s)
nms <- gsub('masie_ice_r00_v01_', '', nms)
nms <- gsub('_4km', '', nms)
循环匹配姓名，并将结果放入列表
p <- list()
for (i in 1:length(b)) {
j <- which(b$ord.year[i] == nms)
r <- s[[j]]
z <- crop(r, b[i,])
z <- mask(z, b[i,])
p[[i]] <- PatchStat(z)
}
请注意，p 的每个元素都有一个包含多行和多列的 data.frame。
p[[1]]
#patchID n.cell n.core.cell n.edges.perimeter n.edges.internal area core.area perimeter perim.area.ratio shape.index frac.dim.index core.area.index
#1 1 53 5 68 144 53 5 68 1.2830189 2.266667 1.427207 0.09433962
#2 2 123 8 182 310 123 8 182 1.4796748 3.956522 1.586686 0.06504065
#3 3 149 31 190 406 149 31 190 1.2751678 3.800000 1.543074 0.20805369
#4 4 54 2 114 102 54 2 114 2.1111111 3.800000 1.679578 0.03703704
#5 5 337 206 146 1202 337 206 146 0.4332344 1.972973 1.236172 0.61127596
如果您只想要第一行
pp <- t(sapply(p, function(i) i[1,]))
将它与原来的 data.frame 结合起来现在变得微不足道
dfpp <- cbind(df, pp)