Search code examples
rterra

How does the R package `terra` decide on the template raster properties when projecting?


I am trying to understand how the R package terra decides on the properties of the new raster it creates when projecting. I know best practice when projecting is to use a template raster in the new crs, e.g. project(my_raster, my_template_raster). However, if the user does try to project just using a crs e.g. project(my_raster,"ESRI:54009"), it is not clear how terra decides on the dimensions, extent and resolution for the template raster in the new projection. For example, in the reprex below, the projected raster does not cover the extent of the projected version of the polygonized raster which seems strange.

library(terra)
#> terra 1.7.71

r <- rast(xmin = 0, xmax = 5, ymin = 0, ymax = 5, resolution = 1)
values(r) <- 1:ncell(r)

#Mollweide
proj_moll <- "ESRI:54009"

#create vector grid of raster
v_moll <- r |>
  as.polygons(dissolve = FALSE) |>
  project(proj_moll)

#project raster to Molleweide
r_moll <- project(r, proj_moll)

#plot: setting limits manually since ext = argument doesn't seem to function correctly at the moment - see https://github.com/rspatial/terra/issues/1495
plot(r_moll, xlim = c(ext(v_moll)$xmin, ext(v_moll)$xmax), ylim = c(ext(v_moll)$ymin, ext(v_moll)$ymax))
lines(v_moll, col = "red", lwd = 2)

Created on 2024-05-07 with reprex v2.1.0


Solution

  • The algorithm first matches the (average) resolution, and tries to match the extent. Since you can only have entire rows and columns you will always get an under or an overshoot. Adding another row and/or column would lead to a lot of overshoot.

    The effect is hardly visible if you use smaller cells

    library(terra)
    r <- rast(xmin = 0, xmax = 5, ymin = 0, ymax = 5, resolution = 1)
    values(r) <- 1:ncell(r)
    proj_moll <- "ESRI:54009"
    v_moll <- as.polygons(r, dissolve = FALSE) |> project(proj_moll)
    
    # 4x disaggregation to get smaller cells.
    r_moll <- disagg(r, 4) |> project(proj_moll)
    
    plot(r_moll, ext=ext(v_moll)+50000)
    lines(v_moll, col = "red", lwd = 2)
    

    enter image description here