循环提取特定的 raster/spatial 点对
Loop to extract specific raster/spatial point pairs
我有许多不重叠的点 shapefile，我想将其归因于类似的栅格，也不重叠。我想用栅格数据来归因于这些点。对于我正在执行此操作的某些栅格数据类型，我能够先合并栅格，然后再合并属性。
但是，我的最后一组栅格数据没有同源，所以我无法 merge/mosaic 它们。我试图在不合并栅格的情况下将这些点归因于栅格。这将要求我在特定空间点 - 栅格对上使用 extract() 。我用一个唯一的 4 字母名称命名了每个空间点文件，这也是我希望 extract() 使用的栅格名称的一部分。
我在下面创建了一个可重现的示例来模拟我的数据和问题。谁能就我如何编写 extract() 的循环代码以获取空间点文件提取到类似命名的栅格提出建议？
或者，如果 better/possible 组合所有空间点并循环提取所有栅格，然后管理数据，以便所有提取的值都在数据帧的一个向量或列中，也许那会更好。
我正在使用 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:
#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 ‘"list"’"
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")
}
你不能通过文件名来匹配它们吗？如果是这样，那就是你应该做的。您提出的建议有不止一种可能，但您似乎正在尝试解决您自己创造的问题 --- 避免所有这些可能更容易。
我会使用
读入文件
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
我有许多不重叠的点 shapefile，我想将其归因于类似的栅格，也不重叠。我想用栅格数据来归因于这些点。对于我正在执行此操作的某些栅格数据类型，我能够先合并栅格，然后再合并属性。 但是，我的最后一组栅格数据没有同源，所以我无法 merge/mosaic 它们。我试图在不合并栅格的情况下将这些点归因于栅格。这将要求我在特定空间点 - 栅格对上使用 extract() 。我用一个唯一的 4 字母名称命名了每个空间点文件，这也是我希望 extract() 使用的栅格名称的一部分。
我在下面创建了一个可重现的示例来模拟我的数据和问题。谁能就我如何编写 extract() 的循环代码以获取空间点文件提取到类似命名的栅格提出建议？
或者，如果 better/possible 组合所有空间点并循环提取所有栅格，然后管理数据，以便所有提取的值都在数据帧的一个向量或列中，也许那会更好。
我正在使用 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:
#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 ‘"list"’"
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")
}
你不能通过文件名来匹配它们吗？如果是这样，那就是你应该做的。您提出的建议有不止一种可能，但您似乎正在尝试解决您自己创造的问题 --- 避免所有这些可能更容易。
我会使用
读入文件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