[GIS] Difference between gdalwarp and projectRaster

coordinate systemgdalgdalwarprraster

I am trying to project a Raster. In R there is the projectRaster() function to to this (below a fully reproducibly example) :

# example Raster
require(raster)
r <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, ncols=40, nrows=40)
r <- setValues(r, 1:ncell(r))
projection(r)
# project to
newproj <- "+init=epsg:4714"


# using raster package to reproject
pr1 <- projectRaster(r, crs = CRS(newproj), method = 'bilinear')

Which works fine. However it is quite slow.

In order to increase speed I though to use gdalwarp instead (with a SSD the cost of reading and writing from/to disk/R are not very high).

However, I cannot reproduce the results of projectRaster() using gdalwarp:

# using gdalwarp to reproject
tf <- tempfile(fileext = '.tif')
tf2 <- tempfile(fileext = '.tif')
writeRaster(r, tf)
system(command = paste(paste0("gdalwarp -t_srs \'", newproj, "\' -r bilinear -overwrite"), 
                       tf,
                       tf2))
pr2 <- raster(tf2)

It seems to work, however the results are different:

# Info
system(command = paste("gdalinfo", 
                       tf))
system(command = paste("gdalinfo", 
                       tf2))

# plots
plot(r)
plot(pr1)
plot(pr2)

#extents
extent(r)
extent(pr1)
extent(pr2)

# PROJ4
proj4string(r)
proj4string(pr1)
proj4string(pr2)

# extract value
take <- SpatialPoints(matrix(c(-100, 50), byrow = T, ncol = 2), proj4string = CRS(newproj))
plot(take, add = TRUE)
extract(pr1, take)
extract(pr2, take)

What am I missing / doing wrong?

Are there other (faster) alternatives to projectRaster()?

Best Answer

Nice and reproducible question. Personally, I'd expect that the reason for the difference is in the implementations of the bilinear reprojection. You can obviously look into source code for the two approaches, but I'd expect that to be a vast overkill.
It appears that the R implementation introduces bigger "errors" / "changes" than the raw GDAL version (atleast in my versions & tests - projectRaster introduces changes around +-0.01 while GDAL gives values around +-0.002).

If you compare both approaches using a nearest neighbor reprojection they match as expected.

Related Question