Suite

Différence entre gdalwarp et projectRaster


J'essaie de projeter un Raster. Dans R il y a leprojetRaster()fonction à cela (ci-dessous un exemple entièrement reproductible):

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

Ce qui fonctionne bien. Cependant c'est assez lent.

Afin d'augmenter la vitesse, j'ai pensé à utilisergdalwarpà la place (avec un SSD, le coût de lecture et d'écriture de/vers le disque/R n'est pas très élevé).

Cependant, je ne peux pas reproduire les résultats deprojetRaster()en utilisantgdalwarp:

# utiliser gdalwarp pour reprojeter tf <- tempfile(fileext = '.tif') tf2 <- tempfile(fileext = '.tif') writeRaster(r, tf) system(command = paste(paste0("gdalwarp -t_srs '" , newproj, "' -r bilinéaire -écraser"), tf, tf2)) pr2 <- raster(tf2)

Cela semble fonctionner, mais les résultats sont différents :

# 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) # extraire la valeur prendre <- SpatialPoints(matrix(c(-100, 50), byrow = T, ncol = 2), proj4string = CRS (newproj)) plot(take, add = TRUE) extract(pr1, take) extract(pr2, take)

Qu'est-ce que j'oublie/je fais mal ?

Existe-t-il d'autres alternatives (plus rapides) àprojetRaster()?


Question agréable et reproductible. Personnellement, je m'attendrais à ce que la raison de la différence soit dans les implémentations de la reprojection bilinéaire. Vous pouvez évidemment examiner le code source pour les deux approches, mais je m'attendrais à ce que ce soit une grande exagération.
Il semble que l'implémentation R introduit des "erreurs" / "changements" plus importants que la version brute de GDAL (au moins dans mes versions et tests - projectRaster introduit des changements autour de +-0,01 tandis que GDAL donne des valeurs autour de +-0,002).

Si vous comparez les deux approches à l'aide d'une reprojection du voisin le plus proche, elles correspondent comme prévu.