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

Question

I'm using this function to calculate the length of linestring per cell by ID and store in a list, convert each element of the list into a RasterLayer and turn that list into a RasterBrick. Notice that inside the function I use "left_join()" so that all rasters have the same number of cells. However, the final raster converts to RasterStack instead of RasterBrick, why does this happen?

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)

I also have a sf grid created from points:

#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

You use the argument bylayer=TRUE asking for a separate file for each layer. A RasterBrick can only be created from a single file; hence you get a RasterStack.

Also, because you have a list of RasterLayers, it is more efficient to use stack here instead of brick.

 s <- stack(list_length_raster2 )

This is because a RasterStack essentially is a list or RasterLayers. With brick you would combine all the values into a new structure, and that can be expensive.

And if you want a RasterBrick, do:

 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 

Or use terra instead of raster --- That is easier because the terra SpatRaster encompasses raster's RasterLayer, Stack and Brick objects.

Below I show how you may do this with terra (although I am not quite following all the logic of your script; my understanding is that you want to measure, for each grid cell, the length of the lines that cross it). You need terra >= 1.5-29 that you can install with install.packages('terra', repos='https://rspatial.r-universe.dev')

Your example data

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)

Create a SpatVector and a 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

Apply the rasterizeGeom method to each line

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