【问题标题】:Extract values from raster and append to existing dataframe从栅格中提取值并附加到现有数据框
【发布时间】:2019-09-01 12:17:39
【问题描述】:

我正在尝试从栅格堆栈中提取值并将其附加到现有数据框。这些值是一组指标(来自 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 &lt;- 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 循环中,以最大程度地减少错误空间和处理时间。

非常感谢!

【问题讨论】:

  • 这很简单,但与所有 R 问题一样,您需要提供一些代码和简单的示例数据(由代码或从 R 包创建)才能回答。
  • facepalm,这里是菜鸟。我已经进行了相应的编辑。

标签: r spatial r-raster sp


【解决方案1】:

感谢您提供示例数据。但它仍然不完整(它指的是我们没有的文件。你可以到这个

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)

【讨论】:

  • 非常感谢,@RobertHijmans!第一部分运作良好。现在我只需要将 PatchStat 结果附加到我现有的数据框 (mydata),最好是在 for 循环中
  • 添加了一点
【解决方案2】:

好吧,我做了一件非常丑陋的事情,得到了我想要的。但我不喜欢它。如果有人有更好的主意,我很想听听!

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

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2014-04-15
    • 1970-01-01
    • 1970-01-01
    • 2018-10-31
    • 2016-08-11
    • 2019-12-31
    • 2019-04-24
    • 1970-01-01
    相关资源
    最近更新 更多