Diferencia entre gdalwarp y projectRaster

Question

Estoy intentando proyectar un Raster. En R existe la función projectRaster() a esto (abajo un ejemplo totalmente reproducible) :

# 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')

Que funciona bien. Sin embargo, es bastante lento.

Para aumentar la velocidad he pensado en utilizar gdalwarp en su lugar (con un SSD el coste de lectura y escritura desde/hacia el disco/R no es muy alto).

Sin embargo, no puedo reproducir los resultados de projectRaster() utilizando 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)

Parece que funciona, pero los resultados son diferentes:

# 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)

¿Qué me estoy perdiendo / haciendo mal?

¿Existen otras alternativas (más rápidas) a projectRaster() ?

Answer 1

Bonita y reproducible pregunta. Personalmente, supongo que la razón de la diferencia está en las implementaciones de la reproyección bilineal. Obviamente, se puede buscar en el código fuente de los dos enfoques, pero yo esperaría que eso sea una gran exageración.
Parece que la implementación de R introduce mayores "errores" / "cambios" que la versión cruda de GDAL (al menos en mis versiones y pruebas - projectRaster introduce cambios alrededor de +-0,01 mientras que GDAL da valores alrededor de +-0,002).

Si se comparan ambos enfoques utilizando una reproyección del vecino más cercano, coinciden como se esperaba.