在R中的每个多多边形特征上循环使用 "栅格 "包函数。

Question

我试图循环使用 "栅格 "包中的几个函数，即crop()、mask()、reclassify()和unstack()as.list()。我有十个栅格图层，它们共享相同的范围和数据类型；它们对应于10个时间点的土地覆盖。我想为crop()->mask()->reclassify()->as.list()的每个输出创建单独的列表变量。我能够对1个多边形特征进行管道处理，但我需要能够对存储在multipolygon Shapefile中的10个多边形特征中的每一个进行循环，这样我就可以根据指定的命名惯例保存每个输出列表。

谢谢你，请你提供建议。下面我分享一下我的代码。

EDIT: 我想知道for-loop是否是正确的方法，还是lapply方法更好？

# Load libraries
library(raster)             # for raster processing
library(rgdal)              # for raster/vector processing
library(sf)                 # for Shapefile processing

# Stack 10 rasters together
raster.stack = stack( 
  raster("path/raster1.tif"),
  raster("path/raster2.tif"),
  raster("path/raster3.tif"),
  raster("path/raster4.tif"),
  raster("path/raster5.tif"),
  raster("path/raster6.tif"),
  raster("path/raster7.tif"),
  raster("path/raster8.tif"),
  raster("path/raster9.tif"),
  raster("path/raster10.tif")
)

# Prepare reclassification codes from 9-class raster to 3-class raster
reclasscodes = c(
  0,0, # no data 
  1,1,
  2,1,
  3,1,
  4,1,
  5,2,
  6,2,
  7,3,
  8,3,
  9,3
)

# Convert reclass codes list into n x 2 matrix
reclassmatrix = matrix(reclasscodes, ncol=2, byrow = T)

# Load multipolygon vector Shapefile
multipolygon = shapefile("path/multipolygon.shp") # Shapefile is made of n polygons

# Example subset Shapefile to polygon_1 using attribute "ID"
polygon_1 = subset(multipolygon,ID=="D-4")

# Create output for polygon_1
list_polygon_1 =
  raster.stack %>%
  crop(y = polygon_1) %>%              # crop to bounds
  mask(mask = polygon_1) %>%           # mask to polygon cutline
  reclassify(rcl = reclassmatrix) %>%  # reclassify to 3-class
  as.list() # functions the same as unstack() where raster brick is converted to list of raster layers
# I use %>% because I do not want to save any of the intermediate outputs.
# Resultant output is a variable list for polygon_1 named 'list_polygon_1' which is exactly what I want.
# Worked perfectly.

# How do I repeat this process for polygon_1 to polygon n?

# My attempt
for (i in 1:nrow(multipolygon)) {
  raster.stack %>%
  crop(y = multipolygon[i,]) %>%
  mask(mask = multipolygon[i,]) %>%
  reclassify(rcl = reclassmatrix) %>%
  as.list() %>% # up till here it is the same steps as before for polygon_1
  # now I want to save each list output as a separate variable according to i, e.g. list_polygon_2, list_polygon_3 etc.
  assign(paste(multipolygon$ID, i, sep = '_')) # assign a naming convention for each output variable 
}
# Does not work. Even without the last line of code "..assign(paste(...))" there is no output variable from the as.list() line.

Answer 1

请总是包含一个最小的自足的可重复的例子。

示例数据

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster")) 

xy1 <- xy2 <- xy3 <- matrix(c(10,17, 6,10,71,60,62,71), ncol=2)
xy2[,1] <- xy2[,1] + 30
xy3[,2] <- xy3[,2] - 30
p <- spPolygons(xy1, xy2, xy3)
#plot(r, 1)
#lines(p)

您的需求

rm = matrix(c(0,100,0,100,150,2,150,255,3), ncol=3, byrow=TRUE)

out <- list()
for (i in 1:length(p)) {
    x <- crop(s, p[i,])
    x <- mask(x, p[i,])
    out[[i]] <- reclassify(x, rm)
}

你说的unstack没有意义（unlist也没有用）。我建议不要这样做，但是你可以做的是

out2 <- lapply(out, unstack)

我不确定你真正想要的是什么。如果你想知道单元格的值，你可以做得更简单一些（不需要循环），然后做

r <- reclassify(s, rm)
e <- extract(r, p)

对于你关于lapply和循环的问题。lapply可以简明扼要，但在这种情况下，写一个循环更好，因为它更容易读，也更容易写，特别是当你不使用 %>%.