Expanding the extent of a raster while keeping original cell values

Question

Due to the type of models I'm running, I need to expand my rasters beyond their original extent. For example, my raster has

ext(original_raster)
SpatExtent : -2700100, 2750900, -2500900, 3342100 (xmin, xmax, ymin, ymax)

and I want to expand it to +500 cells on each side

extend_cells <- 500
extended_extent <- c(floor(ext(original_raster)[1]) - extend_cells * original_res[1],
                      ceiling(ext(original_raster)[2]) + extend_cells * original_res[1],
                      floor(ext(original_raster)[3]) - extend_cells * original_res[2],
                      ceiling(ext(original_raster)[4]) + extend_cells * original_res[2])

# Create a new raster with nodata values
new_raster <- rast(nrow=original_dims[1] + extend_cells * 2,
                   ncol=original_dims[2] + extend_cells * 2,
                   xmin=extended_extent[1], xmax=extended_extent[2],
                   ymin=extended_extent[3], ymax=extended_extent[4],
                   crs=crs(original_raster))
new_raster[] <- NA

which then gives me

ext(new_raster)
SpatExtent : -3200100, 3250900, -3000900, 3842100 (xmin, xmax, ymin, ymax)

But then I need to overlay the cells that I have from the original raster with their values on the new extended raster, and that's where I run into problems. Mostly because it feels nearly impossible to use any function in terra because my rasters have different extents (which is the whole point). I think the way to go is to replace each cell by calling each cell individually with an ifelse statement. So I did it like this:

new_raster[] <- ifelse(!is.na(original_raster[]), original_raster[], NA)

But then, it changes the values in my raster, as you can see here: (vertices is just a sf object with coordinates from my model)

vertices <- fm_vertices(mesh)
real_values <- terra::extract(original_raster, vertices)
new_values <- terra::extract(new_raster, vertices)
head(real_values)
  ID         slope
1  1            NA
2  2            NA
3  3 -0.0006669735
4  4 -0.0004087704
5  5  0.0003122997
6  6 -0.0003670746
> head(new_values)
  ID         slope
1  1 -0.0004545018
2  2 -0.0535268076
3  3  0.0047381981
4  4  0.0008052177
5  5 -0.0083403271
6  6  0.0005182655

Does anyone have any suggestions on how I can do this? I also accept ideas on how to do this in python (if that would work better) and I have also tried transforming everything into points with sf. I'm running out of ideas...

Answer 1

Since you are assigning NA's to the new extent raster then you can use terra::merge or simplify the workflow and just use terra::extend

library(terra)

or <- rast(ext(-2700100, 2750900, -2500900, 3342100), 
          resolution = 1000) 
  or[] <- runif(ncell(or))

extend_cells <- 500

Using terra::extend

original.extend <- extend(or, c(extend_cells,extend_cells))
  plot(ext(original.extend), lwd=2, col="red")
    plot(original.extend, colNA=NA, add=TRUE)
      plot(ext(or), lwd=2, add=TRUE)

nrow(original.extend) - nrow(or)
ncol(original.extend) - ncol(or)

Compared to your method of expanding the extent and then using terra::merge

extended_extent <- c(floor(ext(or)[1]) - extend_cells * res(or)[1],
                      ceiling(ext(or)[2]) + extend_cells * res(or)[1],
                      floor(ext(or)[3]) - extend_cells * res(or)[2],
                      ceiling(ext(or)[4]) + extend_cells * res(or)[2])

# Create a new raster with nodata values
new_raster <- rast(nrow=nrow(or) + extend_cells * 2,
                   ncol=ncol(or) + extend_cells * 2,
                   xmin=extended_extent[1], xmax=extended_extent[2],
                   ymin=extended_extent[3], ymax=extended_extent[4],
                   crs=crs(or))
  new_raster[] <- NA
  
rextend <- merge(or, new_raster) 
  plot(ext(rextend), lwd=2, col="red")
    plot(rextend, colNA=NA, add=TRUE)
      plot(ext(or), lwd=2, add=TRUE)

length(or[is.na(or)]) 
length(rextend[is.na(rextend)])