r - 在 R 代码中使用嵌套 for 循环以减少冗余

Question

编辑:栅格是从国家气象局( https://water.weather.gov/precip/download.php )获取的。我的县形状文件适用于肯塔基州，可以从 Esri 的 ArcGIS 在线数据服务获取。

我有 5 个栅格，其中包含不同记录期间的降水数据:每日、前 30 天、60 天、90 天和 365 天。我想执行区域统计来计算代表县边界多边形的 shapefile 内的每个栅格的平均值。我有 120 个县多边形。

我有有效的代码，但我想改进它。目前，它涉及大量复制和粘贴相同的代码。

library(rgdal)
library(raster)
library(sf)
library(maptools)

# County shapefiles
input_path <- "C:/path/to/shapefiles/GIS_data/Counties"
files <- list.files(input_path, pattern="[.]shp$", full.names=TRUE)
allShapes <- lapply(files, readOGR)
filenames <- list.files(input_path, pattern="[.]shp$")

# Rasters
NWS1_daily <- raster(paste(daily_downloads, daily_fname, sep='/'), band = 1)
NWS1_prev30 <- raster(paste(bulk_downloads, prev30_fname, sep='/'), band = 1)
NWS1_prev90 <- raster(paste(bulk_downloads, prev90_fname, sep='/'), band = 1)
NWS1_prev120 <- raster(paste(bulk_downloads, prev120_fname, sep='/'), band = 1)
NWS1_prev365 <- raster(paste(bulk_downloads, prev365_fname, sep='/'), band = 1)

# Perform zonal statistics to acquire mean observation for each county
observations <- vector()
filenames <- vector()
for (i in 1:length(allShapes)){
  observations[i] <- c(extract(NWS1_daily, allShapes[[i]], fun=mean, na.rm=TRUE, weights=T, normalizeWeights=F))
  filenames<-files[i]
}
filenames <- list.files(input_path, pattern="[.]shp$") 
observations<-data.frame(observations,filenames,stringsAsFactors=FALSE)

observations <- vector()
filenames <- vector()
for (i in 1:length(allShapes)){
  observations[i] <- c(extract(NWS1_prev30, allShapes[[i]], fun=mean, na.rm=TRUE, weights=T, normalizeWeights=F))
  filenames<-files[i]
}
filenames <- list.files(input_path, pattern="[.]shp$") 
observations<-data.frame(observations,filenames,stringsAsFactors=FALSE)

observations <- vector()
filenames <- vector()
for (i in 1:length(allShapes)){
  observations[i] <- c(extract(NWS1_prev90, allShapes[[i]], fun=mean, na.rm=TRUE, weights=T, normalizeWeights=F))
  filenames<-files[i]
}
filenames <- list.files(input_path, pattern="[.]shp$") 
observations<-data.frame(observations,filenames,stringsAsFactors=FALSE)

observations <- vector()
filenames <- vector()
for (i in 1:length(allShapes)){
  observations[i] <- c(extract(NWS1_prev120, allShapes[[i]], fun=mean, na.rm=TRUE, weights=T, normalizeWeights=F))
  filenames<-files[i]
}
filenames <- list.files(input_path, pattern="[.]shp$") 
observations<-data.frame(observations,filenames,stringsAsFactors=FALSE)

observations <- vector()
filenames <- vector()
for (i in 1:length(allShapes)){
  observations[i] <- c(extract(NWS1_prev365, allShapes[[i]], fun=mean, na.rm=TRUE, weights=T, normalizeWeights=F))
  filenames<-files[i]
}
filenames <- list.files(input_path, pattern="[.]shp$") 
observations<-data.frame(observations,filenames,stringsAsFactors=FALSE)

我如何将嵌套的 for 循环放入现有代码中以解决一些冗余问题？

我一直在尝试使用一些在线资源作为引用来做到这一点，但一直无法让它发挥作用。

raster_names_band1 <- c(NWS1_daily, NWS1_prev30, NWS1_prev90, NWS1_prev120, NWS1_prev365)

# Perform zonal statistics on each county to acquire mean observation for each countt
observations <- vector()
raster_names_band1 <- vector()
filenames <- vector()
for (i in 1:length(allShapes))
  {
for (j in 1:length(raster_names_band1))
  observations[i] <- c(extract(raster_names_band1[[j]], allShapes[[i]], fun=mean, na.rm=TRUE, weights=T, normalizeWeights=F))
  filenames<-files[i]
}
filenames <- list.files(input_path, pattern="[.]shp$") 
observations<-data.frame(observations,filenames,stringsAsFactors=FALSE)

## Error in raster_names_band1[[j]] : subscript out of bounds

## Error in data.frame(observations, filenames, stringsAsFactors = FALSE) : 
## arguments imply differing number of rows: 0, 120

感谢您提出的任何建议!

Answer 1

因为RasterLayers ( raster() 的返回)是特殊的 S4 类对象，将它们与 c() 绑定(bind)在一起。如 c(NWS1_daily, NWS1_prev30, NWS1_prev90, NWS1_prev120, NWS1_prev365)可能不会产生五个对象的列表，但可能会产生不同维度的简化数组。使用 str(raster_names_band1) 检查此对象。对于复杂对象，请使用 list()在集合中绑定(bind)在一起。

此外，不要嵌套 for需要记录初始化对象并分配给它们的循环，请考虑嵌套的 apply 系列方法。具体来说，运行lapply迭代光栅对象并嵌套 sapply迭代allShapes来构建data.frame中的观察称呼。由于所有文件名都相同，因此将其作为全局对象分配一次。甚至用 setNames 结束通话呈现命名列表。

# LIST OF RASTER OBJECTS (USE list() NOT c())
raster_list <- list(NWS1_daily, NWS1_prev30, NWS1_prev90, NWS1_prev120, NWS1_prev365)
filenames <- list.files(input_path, pattern="[.]shp$"))

# USER-DEFINED FUNCTION
df_build <- function(obj_raster) {
   observations <- sapply(allShapes, function(s) 
                               c(extract(obj_raster, s, fun=mean, na.rm=TRUE, 
                                         weights=TRUE, normalizeWeights=FALSE))
                          )
   data.frame(observations, filenames, stringsAsFactors=FALSE)
} 

# BUILD NAMED LIST OF DFs
df_list <- setNames(lapply(raster_list, df_build),
                    c("NWS1_daily", "NWS1_prev30", "NWS1_prev90", 
                      "NWS1_prev120", "NWS1_prev365")
                   )

# DF OUTPUTS
df_list$NWS1_daily
df_list$NWS1_prev30
df_list$NWS1_prev90
...

---

事实上，您可以运行嵌套 sapply (包装到 lapply )以避免 list光栅对象和 setNames但使用栅格对象名称的字符向量和 get()通过字符串引用对象:

raster_names <- c("NWS1_daily", "NWS1_prev30", "NWS1_prev90",
                  "NWS1_prev120", "NWS1_prev365")
filenames <- list.files(input_path, pattern="[.]shp$"))

df_build <- function(raster_nm) {
   # keep all same but change: extract(obj_raster, s, ... --> extract(get(raster_nm), s, ...
}

# BUILD NAMED LIST OF DFs
df_list <- sapply(raster_names, df_build, simplify = FALSE)

注意:以上是未经测试的代码，可能需要根据全局对象进行各种调整。