Sunday 22 December 2019

raster - Difference between gdalwarp and projectRaster


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 comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...