Different output between terra and raster packages when generating a raster of a template grid of a specified resolution

[kbondo1]

Hi,

I was trying to update some of my code from the raster to the terra package to generate 5km rasters using a grid template of my study area, and noticed I was obtaining quite different results when resampling and extracting spatial data from these rasters between the two packages.

I am thinking the following issue may be causing the problem and was wondering your thoughts.

To generate 5 km rasters, we made a raster of a grid template of the study area with a specified resolution in each grid cell and then projected the rasters to the grid template and used bilinear interpolation using both packages.

However, when I generated the grid template using the raster package, I obtained 5546 grid cells and when using the terra package I obtained 5555 grid cells. The extents of the grid templates also slightly differ and when generating the raster of the template of the study area using each package. When I zoom in on the grid templates generated using each package, it appears there are differences in how the grid cells between each method are overlapping. Although this difference seems small, when extracting one variable over several rasters, I obtained quite different results using each method.

Is there a way to make this step more comparable between packages, or is one package using a better method than another to generate the raster from the grid template?

I am using macOS Sonomo 14.5 operating system, terra package version 1.7-39 and gdal "3.5.3", proj "9.1.0", and geos "3.11.0".

An R Markdown file showing results between each package is attached.

Raster and Terra.md

[rhijmans]

Thank you. Can you please include the example as R code (not a markdown file)?

[kbondo1]

Below is the R code:

# Install Libraries
library("terra")
library("raster")
library("tigris")
library("sf")

#Raster Package
nad83.2011.pa.north="+proj=lcc +lat_1=41.95 +lat_2=40.88333333333333 +lat_0=40.16666666666666 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs"

# chosen resolution:
resolution=5000

states <- tigris::states(resolution = "5m",
                         class = "sf",
                         progress_bar = FALSE)

PA=states[states$STUSPS=="PA",]

template=as_Spatial(st_buffer(st_transform(PA, crs(nad83.2011.pa.north)), dist = 10000)) # template 10 km out of PA borders

#create raster of specified resolution and projection
grid.template <- raster(extent(template), 
                        resolution = c(resolution,resolution),
                        crs =  crs(template))

#move it to a spatial object (before it was a empty raster)
gridPolygon <- rasterToPolygons(grid.template)

# crop the grid by the box I created around PA. Now template is the gridpolygon 
template=gridPolygon[template,]

# add cell ID
template@data$layer=c(1:nrow(template@data))

colnames(template@data)[1]="cell_number"

nrow(template@data) # 5546 cells

# Check point: plot the template, PA

plot(template, col=NA)
plot(spTransform(as_Spatial(PA), nad83.2011.pa.north), add=T, border="red", lwd=2)

#Terra Package

nad83.2011.pa.north="+proj=lcc +lat_1=41.95 +lat_2=40.88333333333333 +lat_0=40.16666666666666 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs"

# Chosen resolution:
resolution=5000

states <- tigris::states(resolution = "5m",
                         class = "sf",
                         progress_bar = FALSE)

PA=states[states$STUSPS=="PA",]

template1=st_buffer(st_transform(PA, st_crs(nad83.2011.pa.north)), dist = 10000)# template 10 km out of PA borders

# Create raster of specified resolution and projection
grid.template1 <- rast(ext(template1), 
                       res = c(resolution,resolution),
                       crs =  crs(template1))

# Move it to a spatial object (before it was a empty raster)
gridPolygon <- as.polygons(grid.template1)

# Crop the grid by the box I created around PA. Now template is the gridpolygon 
template1=st_as_sf(gridPolygon)[template1,]

# Add cell ID
template1$layer=c(1:nrow(template1))
colnames(template1)[2]="cell_number"
nrow(template1) # 5555 cells.  

# Check point: plot the template, PA
plot(st_geometry(template1), add = TRUE)
plot(st_transform(PA, st_crs(nad83.2011.pa.north)), add=T, border="red", lwd=2)

#Overlay both template grids from raster and terra packages
plot(st_geometry(template1), col = "red")
plot(template, add = TRUE)
#plot(st_transform(PA, st_crs(nad83.2011.pa.north)), add=T, border="red", lwd=2)

ext(grid.template)

ext(grid.template1)

[rhijmans]

Your example does not work anymore because it depends on package rgdal. It also seems unnecessarily complex.

Perhaps this captures it:

pa.north = "+proj=lcc +lat_1=41.95 +lat_2=40.88333333333333 +lat_0=40.16666666666666 +lon_0=-77.75 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs"
ex <- c(352506.5, 868021.9, -59593.0, 272854.8 )
resolution=5000
rtmp <- raster::raster(raster::extent(ex),  resolution = resolution,    crs = pa.north)
ttmp <- terra::rast(terra::ext(ex),   res = resolution,   crs =  pa.north)

rtmp
#class      : RasterLayer 
#dimensions : 66, 103, 6798  (nrow, ncol, ncell)
#resolution : 5000, 5000  (x, y)
#extent     : 352506.5, 867506.5, -57145.2, 272854.8  (xmin, xmax, ymin, ymax)
#crs        : +proj=lcc +lat_0=40.1666666666667 +lon_0=-77.75 +lat_1=41.95 +lat_2=40.8833333333333 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs 

ttmp
#class       : SpatRaster 
#dimensions  : 66, 103, 1  (nrow, ncol, nlyr)
#resolution  : 5000, 5000  (x, y)
#extent      : 352506.5, 867506.5, -59593, 270407  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=lcc +lat_0=40.1666666666667 +lon_0=-77.75 +lat_1=41.95 +lat_2=40.8833333333333 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs 

In both cases, the extent is adjusted. This is necessary since the resolution does not perfectly fit the extent:

(ex[4]-ex[3]) / 5000
#[1] 66.48956
(ex[2]-ex[1]) / 5000
#[1] 103.1031

But while both raster and terra keep xmin, and adjust xmax, raster adjusts xmin, whereas terra adjusts xmax.

[kbondo1]

Thanks so much! That helped resolve the issue. Also, thanks for pointing out the line of code that no longer worked with the raster package. I looked into it further and then realized that I had adapted that line of code incorrectly for the new packages. Once I fixed that, I was able to obtain the same number of grid cells that I had received when using the raster package. After changing that line of code and after adjusting the extent as you have shown in the code above, I was able to obtain similar results when using both packages.