使用 shapefile 提取多个栅格像素值并使用半径列表循环
extracting multiple raster pixel values with shapefile and looping it with a list of radiuses
我要解决的问题是我想使用包含多个点位置的 shapefile 提取多个栅格的像素值,对于每个点我想提取多个栅格“时间序列”的值,并循环遍历要在提取函数中输入的缓冲区半径列表,我想要的最终结果是每个缓冲区的每个站点的 csv
到目前为止,我确实设法在 R 下用栅格库编写了一个程序,我的代码由两个循环组成,第一个循环进入每个栅格并提取与一个中多个位置的位置相对应的像素值shapefile,第二个包含要循环提取每个像素半径的均值、最大值和最小值的半径列表
但是此时代码过于松散,没有有效整合radius循环
library(raster)
library(rgdal)
library(stringr)
library(ggplot2)
library(lubridate)
raster_files <- list.files('J:/nvc',pattern="*.tif",full.names = T)
r_name<- list.files(path='J:/nvc',pattern="*.tif",full.names = F)
shape<-readOGR('J:/nvc/locations/locations.shp')
radius <- read.csv('J:/nvc/radius.csv')
rList <- list()
statList <- list()
for (j in 1:nrow(radius)){
for(i in 1:length(raster_files)){
ras <- raster(raster_files[i])
stack <- stack(ras)
crs(stack) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"
mean <- raster::extract(stack,shape,buffer=j,fun=mean,df=TRUE)
min <- raster::extract(stack,shape,buffer=j,fun=max,df=TRUE)
max <- raster::extract(stack,shape,buffer=j,fun=min,df=TRUE)
Name <- r_name[i]
statList[[i]] <- data.frame(max)
df <- do.call(rbind.data.frame,statList)
}
}
N.B : 我只是命名 3 个地点来尝试,alpha bravo 和 charlie
半径为 1,2,4,6,8,10,20 公里
栅格命名为 1h-2020042601 , 1h-2020042602,1h-2020042603 ..... 例如:每小时 1 个栅格
这是一个最小的、可重现的、独立的例子。
您不需要遍历栅格。制作一个 RasterStack
library(raster)
#raster_files <- list.files('J:/nvc',pattern="*.tif",full.names = T)
raster_files <- system.file("external/rlogo.grd", package="raster")
stack <- stack(raster_files)
points <- cbind(x=c(48, 53, 7), y=c(63, 43, 28))
bufs <- c(5, 10)
一种循环方式是 lapply
x <- lapply(bufs, function(b) extract(stack, points, buffer=b))
即returns列表列表。在这种情况下
# 2 buffers
length(x)
#[1] 2
# 3 points
sapply(x, length)
#[1] 3 3
计算所需的统计数据
r <- rapply(x, function(i) c(min(i, na.rm=TRUE), mean(i, na.rm=TRUE), max(i, na.rm=TRUE)))
并很好地组合起来
r <- matrix(r, ncol=3, byrow=TRUE)
colnames(r) <- c("min", "mean", "max")
r <- cbind(buf=rep(bufs, each=nrow(points)),
point=rep(1:nrow(points), length(bufs)), r)
r
# buf point min mean max
#[1,] 5 1 0 119.4375 255
#[2,] 5 2 0 132.6292 255
#[3,] 5 3 81 182.5958 255
#[4,] 10 1 0 156.1361 255
#[5,] 10 2 0 154.1814 255
#[6,] 10 3 74 189.7354 255
我要解决的问题是我想使用包含多个点位置的 shapefile 提取多个栅格的像素值,对于每个点我想提取多个栅格“时间序列”的值,并循环遍历要在提取函数中输入的缓冲区半径列表,我想要的最终结果是每个缓冲区的每个站点的 csv
到目前为止,我确实设法在 R 下用栅格库编写了一个程序,我的代码由两个循环组成,第一个循环进入每个栅格并提取与一个中多个位置的位置相对应的像素值shapefile,第二个包含要循环提取每个像素半径的均值、最大值和最小值的半径列表
但是此时代码过于松散,没有有效整合radius循环
library(raster)
library(rgdal)
library(stringr)
library(ggplot2)
library(lubridate)
raster_files <- list.files('J:/nvc',pattern="*.tif",full.names = T)
r_name<- list.files(path='J:/nvc',pattern="*.tif",full.names = F)
shape<-readOGR('J:/nvc/locations/locations.shp')
radius <- read.csv('J:/nvc/radius.csv')
rList <- list()
statList <- list()
for (j in 1:nrow(radius)){
for(i in 1:length(raster_files)){
ras <- raster(raster_files[i])
stack <- stack(ras)
crs(stack) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"
mean <- raster::extract(stack,shape,buffer=j,fun=mean,df=TRUE)
min <- raster::extract(stack,shape,buffer=j,fun=max,df=TRUE)
max <- raster::extract(stack,shape,buffer=j,fun=min,df=TRUE)
Name <- r_name[i]
statList[[i]] <- data.frame(max)
df <- do.call(rbind.data.frame,statList)
}
}
N.B : 我只是命名 3 个地点来尝试,alpha bravo 和 charlie
半径为 1,2,4,6,8,10,20 公里
栅格命名为 1h-2020042601 , 1h-2020042602,1h-2020042603 ..... 例如:每小时 1 个栅格
这是一个最小的、可重现的、独立的例子。 您不需要遍历栅格。制作一个 RasterStack
library(raster)
#raster_files <- list.files('J:/nvc',pattern="*.tif",full.names = T)
raster_files <- system.file("external/rlogo.grd", package="raster")
stack <- stack(raster_files)
points <- cbind(x=c(48, 53, 7), y=c(63, 43, 28))
bufs <- c(5, 10)
一种循环方式是 lapply
x <- lapply(bufs, function(b) extract(stack, points, buffer=b))
即returns列表列表。在这种情况下
# 2 buffers
length(x)
#[1] 2
# 3 points
sapply(x, length)
#[1] 3 3
计算所需的统计数据
r <- rapply(x, function(i) c(min(i, na.rm=TRUE), mean(i, na.rm=TRUE), max(i, na.rm=TRUE)))
并很好地组合起来
r <- matrix(r, ncol=3, byrow=TRUE)
colnames(r) <- c("min", "mean", "max")
r <- cbind(buf=rep(bufs, each=nrow(points)),
point=rep(1:nrow(points), length(bufs)), r)
r
# buf point min mean max
#[1,] 5 1 0 119.4375 255
#[2,] 5 2 0 132.6292 255
#[3,] 5 3 81 182.5958 255
#[4,] 10 1 0 156.1361 255
#[5,] 10 2 0 154.1814 255
#[6,] 10 3 74 189.7354 255