使用 sf::st_crop() 和 raster::crop() 裁剪光栅堆栈时出错

Posted 2023-02-25

【中文标题】使用 sf::st_crop() 和 raster::crop() 裁剪光栅堆栈时出错

【英文标题】：errors when clipping raster stacks with sf::st_crop() and raster::crop()

【发布时间】：2021-08-14 08:11:08

【问题描述】：

我的代码从高光谱反射光栅堆栈中选择的三层光栅堆栈以及多边形生成地图。

my_rfl <- raster::stack("./filepath/my_raster.bsq")
my_pdk <- read_sf("./data/shp/my_pdk.shp") 

my_stack <- subset(my_rfl, c(18, 30, 68))

my_map <- ggplot() +
  geom_spatial_rgb(data = my_stack, mapping = aes(x = x, 
                                                  y = y,
                                                  r = red,
                                                  g = green,
                                                  b = blue), alpha = 0.6) +
  geom_sf(data = pdk, colour = "blue", fill = "NA")
my_map

我只想查看多边形内的栅格区域，所以我这样做了：

my_stack <- raster::crop(my_stack, my_pdk)

或

my_stack <- st_crop(my_stack, my_pdk)

但是raster::crop()好像丢了红色层，抛出：Error in FUN(X[[i]], ...) : object 'red' not found，当我尝试打印地图的时候。我也试过raster::mask()，结果一样。

st_crop 抛出 Error in UseMethod("st_crop") : no applicable method for 'st_crop' applied to an object of class "c('RasterBrick', 'Raster', 'RasterStackBrick', 'BasicRaster')"

【参考方案1】：

首先：下次做一个可重现的例子，这样更容易帮助你；第二：始终包含您的library，以便清楚您使用的功能来自何处。 但是，我认为问题在于，当您将crop raster::stack 转换为raster::brick 时。请看下面的代码：

library(sf)
library(raster)
library(terrain)
#generate example data
r <- raster(ncol=36, nrow=18, vals=1:(18*36))
s <- stack(r,r*0.5,r/3)
names(s) <- c("red","green","blue")
cds1 <- rbind(c(-150,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-150,-20))
polys <- spPolygons(cds1)
polys <- st_as_sf(polys)

 ggplot() +
  geom_spatial_rgb(data = s, mapping = aes(x = x, 
                                                  y = y,
                                                  r = red,
                                                  g = green,
                                                  b = blue), alpha = 0.6) +
  geom_sf(data = polys, colour = "blue", fill = "NA")

#first crop to the extantion of the polygon
my_stack <- stack(raster::crop(s, polys))
#then mask out (i.e. assign NA) the values outside the polygon
my_stack <- stack(raster::mask(my_stack, polys))

ggplot() +
  geom_spatial_rgb(data = my_stack, mapping = aes(x = x, 
                                           y = y,
                                           r = red,
                                           g = green,
                                           b = blue), alpha = 0.6) +
  geom_sf(data = polys, colour = "blue", fill = "NA")

情节现在看起来像这样 如您所见，raster::stack 的范围与多边形相同

class      : RasterStack 
dimensions : 6, 10, 60, 3  (nrow, ncol, ncell, nlayers)
resolution : 10, 10  (x, y)
extent     : -160, -60, -60, 0  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names      :   red, green,  blue 
min values : 327.0, 163.5, 109.0 
max values : 507.0, 253.5, 169.0 

【讨论】：

感谢您的回答以及为我的 reprex 生成栅格和多边形的方法 - 我不知道该怎么做。