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

Question

我的代码从从高光谱反射比光栅堆栈中选择的三个层的光栅堆栈中生成地图，外加多边形。

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')"

Answer 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