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

# example Raster
r <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, ncols=40, nrows=40)
r <- setValues(r, 1:ncell(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"), 
pr2 <- raster(tf2)

It seems to work, however the results are different:

# Info
system(command = paste("gdalinfo", 
system(command = paste("gdalinfo", 

# plots



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

  • 1
    No one? I supplied a fully reproducible example (should work with Linux or Mac)...
    – EDi
    Commented Jul 20, 2015 at 10:56
  • What are you expecting? Do both options use the same proj.4?
    – user10353
    Commented Jul 20, 2015 at 18:36
  • I expect that both methods yield the same re-projected raster, the same extent and the same value at (-100, 50). However, they apparently don't :(
    – EDi
    Commented Jul 20, 2015 at 21:59
  • 1
    The two programs are creating different grids to warp onto. Even if the bilinear sampling was exactly the same, the points that are being interpolated are in different places, and you'd have different answers. The origins and the pixel sizes are different. You could set some flags in gdalwarp(-te, -tr, etc.) to try and reproduce the R version, and then compare the pixel values and see how different they are.
    – user10353
    Commented Jul 21, 2015 at 19:27
  • I found on multiple occasions that using the -order flag (the "order of polynominal used for warping") on gdalwarp even without using GCPs produced more accurate results.
    – christoph
    Commented Jul 21, 2015 at 21:58

1 Answer 1


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.

  • Thanks for this hint with the projection methods! If I find time I'll take a look deeper into those (However, I'm more familiar with R then with C).
    – EDi
    Commented Jul 21, 2015 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.