为什么使用具有相同单元数的 RasterLayers 的 writeRaster(brick()) 会导致 RasterStack？

Question

我正在使用此函数按 ID 计算每个单元格的线串长度并存储在列表中，将列表的每个元素转换为 RasterLayer 并将该列表转换为 RasterBrick。请注意，在我使用“left_join()”的函数内部，所有栅格都具有相同数量的像元。但是，最终的栅格转换为 RasterStack 而不是 RasterBrick，为什么会这样？

library(tidyverse)
library(sf)
library(raster)
library(purrr)


id <- c("844", "844", "844", "844", "844","844", "844", "844", "844", "844",
        "844", "844", "845", "845", "845", "845", "845","845", "845", "845", 
        "845","845", "845", "845")
lat <- c(-30.6456, -29.5648, -27.6667, -31.5587, -30.6934, -29.3147, -23.0538, 
         -26.5877, -26.6923, -23.40865, -23.1143, -23.28331, -31.6456, -24.5648, 
         -27.6867, -31.4587, -30.6784, -28.3447, -23.0466, -27.5877, -26.8524, 
         -23.8855, -24.1143, -23.5874)
long <- c(-50.4879, -49.8715, -51.8716, -50.4456, -50.9842, -51.9787, -41.2343, 
          -40.2859, -40.19599, -41.64302, -41.58042, -41.55057, -50.4576, -48.8715, 
          -51.4566, -51.4456, -50.4477, -50.9937, -41.4789, -41.3859, -40.2536, 
          -41.6502, -40.5442, -41.4057)
df <- tibble(id, lat, long)


#converting to sf
df.sf <- df %>% 
 sf::st_as_sf(coords = c("long", "lat"), crs = 4326)

我还有一个由点创建的 sf 网格：

#creating grid
xy <- sf::st_coordinates(df.sf)

grid = sf::st_make_grid(sf::st_bbox(df.sf),
                       cellsize = .1, square = FALSE) %>%
 sf::st_as_sf() %>%
 dplyr::mutate(cell = 1:nrow(.)) 

#creating raster
r <- raster::raster(grid, res=0.1) 

#creating linestring to each id
df.line <- df.sf %>% 
  dplyr::group_by(id) %>%
  dplyr::summarize() %>%
  sf::st_cast("LINESTRING") 

#build_length_raster <- function(df.line) {

  intersect_list <- by(
    df.line, 
    df.line$id,
    function(id_df) sf::st_intersection(grid, id_df) %>% 
      dplyr::mutate(length = as.numeric(sf::st_length(.))) %>%
      sf::st_drop_geometry()
  ) 
  
  list_length_grid <- purrr::map(intersect_list, function(grid_id) {
    grid_id %>% dplyr::left_join(x=grid, by="cell") %>% 
      dplyr::mutate(length=length) %>%
      dplyr::mutate_if(is.numeric,coalesce,0)
  })

  list_length_raster <- purrr::map(list_length_grid, function(grid_id) {
    raster::rasterize(grid_id, r, field="length", na.rm=F, background=0)
  })

  list_length_raster2 <- unlist(list_length_raster, recursive=F)
  
  raster_brick <- raster::writeRaster(raster::brick(list_length_raster2 ), 
                                      names(list_length_raster2 ), 
                                      bylayer=TRUE, overwrite=TRUE)

#}

Answer 1

您使用参数 bylayer=TRUE 要求为每个图层创建一个单独的文件。只能从单个文件创建 RasterBrick；因此你得到一个RasterStack。

此外，因为您有一个 RasterLayers 列表，所以在这里使用 stack 比 brick 更有效。

 s <- stack(list_length_raster2 )

这是因为 RasterStack 本质上是一个列表或 RasterLayers。使用 brick，您可以将所有值组合到一个新结构中，这可能很昂贵。

如果您想要 RasterBrick，请执行以下操作：

 b <- raster::writeRaster(s, "test.tif", overwrite=T)
 b
 #class      : RasterBrick 
 #dimensions : 87, 120, 10440, 2  (nrow, ncol, ncell, nlayers)
 #resolution : 0.1, 0.1  (x, y)
 #extent     : -52.0787, -40.0787, -31.71421, -23.01421  (xmin, xmax, ymin, ymax)
 #crs        : +proj=longlat +datum=WGS84 +no_defs 
 #source     : test.tif 
 #names      :   test.1,   test.2 
 #min values :        0,        0 
 #max values : 21762.09, 11923.31

或使用 terra 而不是 raster --- 这更容易，因为 terra SpatRaster 包含 raster 的 RasterLayer、Stack 和 Brick对象。

下面我展示了如何使用 terra 执行此操作（尽管我没有完全遵循您脚本的所有逻辑；我的理解是您想要测量每个网格单元格的长度穿过它的线）。您需要可以安装的 terra >= 1.5-29 install.packages('terra', repos='https://rspatial.r-universe.dev')

您的示例数据

id <- c("844", "844", "844", "844", "844","844", "844", "844", "844", "844",
        "844", "844", "845", "845", "845", "845", "845","845", "845", "845", "845","845", "845", "845")
lat <- c(-30.6456, -29.5648, -27.6667, -31.5587, -30.6934, -29.3147, -23.0538, 
         -26.5877, -26.6923, -23.40865, -23.1143, -23.28331, -31.6456, -24.5648, 
         -27.6867, -31.4587, -30.6784, -28.3447, -23.0466, -27.5877, -26.8524, 
         -23.8855, -24.1143, -23.5874)
long <- c(-50.4879, -49.8715, -51.8716, -50.4456, -50.9842, -51.9787, -41.2343, 
          -40.2859, -40.19599, -41.64302, -41.58042, -41.55057, -50.4576, -48.8715, 
          -51.4566, -51.4456, -50.4477, -50.9937, -41.4789, -41.3859, -40.2536, 
          -41.6502, -40.5442, -41.4057)

创建一个 SpatVector 和一个 SpatRaster

library(terra)    
m <- cbind(object=as.integer(as.factor(id)), part=1, x=long, y=lat)
v <- vect(m, type="lines", att=data.frame(id=unique(id)), crs="+proj=longlat")
r <- rast(v, res=1) # lower resolution for example

对每一行应用rasterizeGeom方法

x <- lapply(1:length(v), 
       \(i) rasterizeGeom(v[i], r, "km")
     )
x <- rast(x) 
names(x) <- v$id

x
#class       : SpatRaster 
#dimensions  : 9, 12, 2  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : -51.9787, -39.9787, -31.6456, -22.6456  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : memory  
#              memory  
#names       :      844,      845 
#min values  :        0,        0 
#max values  : 276.6844, 313.1531

为什么使用具有相同单元数的 RasterLayers 的 writeRaster(brick()) 会导致 RasterStack？

Why writeRaster(brick()) using RasterLayers with same number of cells result to RasterStack?

r

layer

geospatial

rasterize

r-raster