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