Expanding the extent of a raster while keeping original cell values

rrasterterra

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…

Best Answer

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)])
Related Question