Trying to create buffer around buildings, but terra::buffer fails for larger rasters

Question

I'm trying to create a raster with info on whether or not each point is near a building. To this end, I'm starting with a geotagged raster that contains the locations of buildings (resolution 10 meters) across a region, extracted from categorized satellite data, and using terra's buffer() function to essentially extend the size of each building by a few hundred meters.

The input raster file (Dropbox link) is about 40 Mb, and consists mostly of NA cells, with all buildings represented by the value 1. Here I try it on a subset of the input, and it works fine:

# Load input
library(terra)

test.buildings <- rast("Z_buildings_raw.tif")
    
# Try a subset
test.buildings2 <- crop(test.buildings, ext(c(400000,450000,6900000,6950000)))
test.buffered2 <- buffer(test.buildings2, width = 300)

plot(test.buildings2, main = "Cropped input")
plot(test.buffered2, main = "Cropped output")

But when I try using the full input file... there are no error messages or warnings, but something clearly goes wrong, because I just get a mass of undifferentiated pixels:

test.buffered <- buffer(test.buildings, width = 300)

plot(test.buildings, main = "Full input")
plot(test.buffered, main = "Full output")

I've tried cropping to different extents, and it seems that the problem appears when the input to buffer() is above a certain size. If anyone has an alternative method I can use, I'd be happy to try that instead.

Answer 1

If I had the choice, I'd probably try to get the desired buffers based on the possibly available buildings vector data from which your "Z_buildings_raw.tif" eventually was created and rasterize the result afterwards, since I'd expect this to be more efficient.

However, your question seems related to the memory allocation issues {terra} is having here and there when processing huge (we're dealing with 1 bln cells here) raster datasets.

What I did was to split your raster in n tiles with a little overlap, apply buffer() on the individual tiles and merge them again to a single raster of full extent. Since this was the first time performing such an operation, there probably are more elegant ways. Also, maybe there is a terra built-in function to be used instead of lapply() but haven't used the app-family very often yet.

Anyway, this took around 4 minutes on my machine and the result looks quite promising, but better check this in detail yourself.

library(terra)
#> terra 1.7.71

# read data
r <- rast("Z_buildings_raw.tif")
r
#> class       : SpatRaster 
#> dimensions  : 39613, 25322, 1  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 348210, 601430, 6825750, 7221880  (xmin, xmax, ymin, ymax)
#> coord. ref. : SWEREF99 TM (EPSG:3006) 
#> source      : Z_buildings_raw.tif 
#> name        : Layer_1 
#> min value   :       1 
#> max value   :       1

# target value for tiling is approx. 5000 x 5000 cells, since this worked for you
x <- rast(nrows = 8, 
          ncols = 5,
          crs = crs(r),
          ext = ext(r))
x
#> class       : SpatRaster 
#> dimensions  : 8, 5, 1  (nrow, ncol, nlyr)
#> resolution  : 50644, 49516.25  (x, y)
#> extent      : 348210, 601430, 6825750, 7221880  (xmin, xmax, ymin, ymax)
#> coord. ref. : SWEREF99 TM (EPSG:3006)

# with a res of 10 m and a buffer width of 300 m, we would need 30 extra cols/rows
getTileExtents(r, x, buffer = 30) |> head(5)
#>         xmin   xmax    ymin    ymax
#>  [1,] 348210 399150 7172060 7221880
#>  [2,] 398550 449800 7172060 7221880
#>  [3,] 449200 500440 7172060 7221880
#>  [4,] 499840 551090 7172060 7221880
#>  [5,] 550490 601430 7172060 7221880

# make tiles
filename <- paste0(tempfile(), "_.tif")
ff <- makeTiles(r, x, filename, buffer = 30)
head(ff, 5)
#>  [1] "...\\AppData\\Local\\Temp\\Rtmpchlnxf\\file2a48460f2161_1.tif" 
#>  [2] "...\\AppData\\Local\\Temp\\Rtmpchlnxf\\file2a48460f2161_2.tif" 
#>  [3] "...\\AppData\\Local\\Temp\\Rtmpchlnxf\\file2a48460f2161_3.tif" 
#>  [4] "...\\AppData\\Local\\Temp\\Rtmpchlnxf\\file2a48460f2161_4.tif" 
#>  [5] "...\\AppData\\Local\\Temp\\Rtmpchlnxf\\file2a48460f2161_5.tif" 

# read individual tiles to a list of SpatRaster objects
ll <- lapply(ff, FUN = rast)

# apply buffer on your tiles
result <- lapply(ll, FUN = buffer, width = 300)

# combine (extended) tiles to a full dataset making use of `max()` in overlapping areas
# `merge()` only lets you to choose between first and last value
r2 <- sprc(result) |> mosaic(fun = "max") 
r2
#> class       : SpatRaster 
#> dimensions  : 39613, 25322, 1  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 348210, 601430, 6825750, 7221880  (xmin, xmax, ymin, ymax)
#> coord. ref. : SWEREF99 TM (EPSG:3006) 
#> source      : spat_2a4823112828_10824.tif 
#> varname     : file2a48460f2161_1 
#> name        : Layer_1 
#> min value   :       0 
#> max value   :       1

# inspect region of interest
plot(r2, ext = ext(c(400000, 450000, 6900000, 6950000)))