9

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

  • No one? I supplied a fully reproducible example (should work with Linux or Mac)... – EDi Jul 20 '15 at 10:56
  • What are you expecting? Do both options use the same proj.4? – user10353 Jul 20 '15 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 Jul 20 '15 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 Jul 21 '15 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 Jul 21 '15 at 21:58
10
+50

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 Jul 21 '15 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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