时间戳

首页 > 解决方案 > 循环提取特定的栅格/空间点对

r - 循环提取特定的栅格/空间点对

问题描述

我有许多不重叠的点 shapefile，我想将其归因于类似栅格，也不重叠。我想将这些点归因于栅格数据。对于我正在使用的某些栅格数据类型，我能够先合并栅格，然后再合并属性。但是，我最后一组栅格数据的来源不同，所以我无法合并/拼接它们。我正在尝试将点归因于栅格而不合并栅格。这将要求我在特定的空间点 - 栅格对上使用 extract()。我用唯一的 4 字母名称命名了每个空间点文件，这也是我希望 extract() 使用的栅格名称的一部分。

我在下面创建了一个可重现的示例来模仿我的数据和问题。任何人都可以就我如何编码 extract() 的循环以获取空间点文件以提取到类似命名的栅格提出建议吗？

或者，如果将所有空间点组合起来并循环提取所有栅格，然后管理数据以使所有提取的值都位于数据帧的一个向量或列中，那可能会更好/可能会更好。

我正在使用 RStudio 1.2.1335

注意：我将此问题发布到 GIS Stack Exchange，但没有收到任何答案，希望交叉发布可以。

library(raster)
library(sp)   

#create point shapefiles
loc1 <- data.frame(x = c(-100,-90, -80, -70),
                  y = c(-100,-90, -80, -70))
loc2 <- data.frame(x = c(-100,-90, -80, -70),
                   y = c(100,90, 80, 70))
loc3 <- data.frame(x = c(100,90, 80, 70),
                   y = c(100,90, 80, 70))
coordinates(loc1) <- ~x+y
sp.loc1<- SpatialPoints(loc1,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
coordinates(loc2) <- ~x+y
sp.loc2<- SpatialPoints(loc2,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
coordinates(loc3) <- ~x+y
sp.loc3<- SpatialPoints(loc3,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
plot(sp.loc1)

#create rasters which have a common part of the naming convention of point shapefiles targeting for attribution
loc1_blablabla<- raster(xmn=-200, xmx=0, ymn=-200, ymx=0)
loc1_blablabla
ncell(loc1_blablabla)
#it has 64800 cells
values(loc1_blablabla)<-1:64800
plot(loc1_blablabla, add=TRUE)
plot(sp.loc1, add=TRUE)
loc2_blablabla<- raster(xmn=-200, xmx=0, ymn=0, ymx=200)
values(loc2_blablabla)<-1:64800
loc3_blablabla<- raster(xmn=0, xmx=200, ymn=0, ymx=200)
values(loc3_blablabla)<-1:64800

#plot all, but first create extent poly
#borrowed some code from here:https://gis.stackexchange.com/questions/206929/r-create-a-boundingbox-convert-to-polygon-class-and-plot/206952
library(rgeos)
ext.box<-rgeos::bbox2SP(n=250, s=-250, w=-250, e=250)
plot(ext.box)
plot(loc1_blablabla, add=TRUE)
plot(loc2_blablabla, add=TRUE)
plot(loc3_blablabla, add=TRUE)
plot(sp.loc1, add=TRUE)
plot(sp.loc2, add=TRUE)
plot(sp.loc3, add=TRUE)

#now attempt to extract raster values at points for multiple non-overlapping point and raster files ie. extract(loc1_blablabla, loc1)
#try lapply as in: https://stackoverflow.com/questions/59164538/create-a-loop-to-extract-data-from-multiple-raster
#create list as below, however, with my real data I would use list.files()
loc.list<-list(loc1,loc2,loc3)
rast.list<-list(loc1_blablabla,loc2_blablabla,loc3_blablabla)
attr.data<-lapply(rastlist.extract,loc.list)
#this doesn't work -- i think I need coordinates as the last term, not a list of spatial points

#also tried a for loop, but this gave error:  "Error in (function (classes, fdef, mtable): unable to find an inherited method for function ‘shapefile’ for signature ‘&quot;list"’&quot;
for (i in 1:length(loc.list)) {
  #this reads in each spatialpoint file and assigns each sp's name to the variable 'sp.name'
  sp.name<-shapefile(loc.list[i]) 
  attr.data<-data.frame(coordinates(sp.name),
                             extract(rast.list[grep(sp.name,rast.list)],sp.name))
  #eventually add this line to affix the new vector attr.data to the coordinates for each plot:  names(attr.data)<-c("x","y","raster.value")
}

标签： rextractraster

解决方案

你不能通过他们的文件名来匹配它们吗？如果是这样，那是你应该做的。您提出的建议不止一种方式，但似乎您正在尝试解决您自己创造的问题——避免所有这些可能会容易得多。

我会使用

library(raster)
r <- lapply(raster_file_names, raster)
s <- lapply(shapefile_names, shapefile)

然后循环这些

使用您的示例数据

library(raster)
p1 <- SpatialPoints(cbind(x = c(-100,-90, -80, -70), y = c(-100,-90, -80, -70)))
p2 <- SpatialPoints(cbind(x = c(-100,-90, -80, -70), y = c(100,90, 80, 70)))
p3 <- SpatialPoints(cbind(x = c(100,90, 80, 70), y = c(100,90, 80, 70)))
pts <- list(p1, p2, p3)

r1 <- raster(xmn=-200, xmx=0, ymn=-200, ymx=0, vals=1:64800)
r2 <- raster(xmn=-200, xmx=0, ymn=0, ymx=200, vals=1:64800)
r3 <- raster(xmn=0, xmx=200, ymn=0, ymx=200, vals=1:64800)
ras <- list(r1, r2, r3)

e <- list()
for (i in 1:length(ras)) {
    e[[i]] <- extract(ras[[i]], pts[[i]])
}
e
#[[1]]
#[1] 32581 29359 26137 22915
#[[2]]
#[1] 32581 35839 39097 42355
#[[3]]
#[1] 32581 35803 39025 42247

推荐阅读