Différence entre gdalwarp et projectRaster

9

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

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

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

Afin d'augmenter la vitesse, je dois utiliser à la gdalwarpplace (avec un SSD, le coût de lecture et d'écriture depuis / vers le disque / R n'est pas très élevé).

Cependant, je ne peux pas reproduire les résultats de l' projectRaster()utilisation 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)

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)

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

Qu'est-ce que je manque / fais mal?

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

raster gdal r coordinate-system gdalwarp EDi
la source
Personne? J'ai fourni un exemple entièrement reproductible (devrait fonctionner avec Linux ou Mac) ...
EDi
Qu'attendez-vous? Les deux options utilisent-elles le même projet.4?
Je m'attends à ce que les deux méthodes produisent le même raster re-projeté, la même étendue et la même valeur à (-100, 50). Cependant, ils ne le font apparemment pas :(
EDi
1
Les deux programmes créent des grilles différentes sur lesquelles se déformer. Même si l'échantillonnage bilinéaire était exactement le même, les points qui sont interpolés sont à des endroits différents et vous auriez des réponses différentes. Les origines et les tailles de pixels sont différentes. Vous pouvez définir des indicateurs dans gdalwarp (-te, -tr, etc.) pour essayer de reproduire la version R, puis comparer les valeurs des pixels et voir à quel point elles sont différentes.
J'ai constaté à plusieurs reprises que l'utilisation du -orderdrapeau («l'ordre du polynôme utilisé pour la déformation») gdalwarpmême sans utiliser les GCP produisait des résultats plus précis.
christoph

Réponses:

10

Belle question 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 regarder dans le code source pour les deux approches, mais je m'attends à ce que ce soit une exagération considérable.
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 en utilisant une reprojection du voisin le plus proche, elles correspondent comme prévu.

Mikkel Lydholm Rasmussen
la source
Merci pour cette astuce avec les méthodes de projection! Si je trouve du temps, j'examinerai plus en profondeur ceux-ci (Cependant, je connais mieux R que C).
EDi