Unterschied zwischen gdalwarp und projectRaster

9

Ich versuche ein Raster zu projizieren. In R gibt es dazu die projectRaster()Funktion (unter einem vollständig reproduzierbaren Beispiel):

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

Welches funktioniert gut. Es ist jedoch ziemlich langsam.

Um die Geschwindigkeit zu erhöhen, verwende ich gdalwarpstattdessen (mit einer SSD sind die Kosten für das Lesen und Schreiben von / auf Festplatte / R nicht sehr hoch).

Ich kann jedoch die Ergebnisse der projectRaster()Verwendung nicht reproduzieren 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)

Es scheint zu funktionieren, aber die Ergebnisse sind unterschiedlich:

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

Was vermisse ich / mache ich falsch?

Gibt es andere (schnellere) Alternativen zu projectRaster()?

raster gdal r coordinate-system gdalwarp EDi
quelle
Niemand? Ich lieferte ein vollständig reproduzierbares Beispiel (sollte mit Linux oder Mac
funktionieren
Was erwarten Sie? Verwenden beide Optionen dasselbe proj.4?
Ich erwarte, dass beide Methoden das gleiche neu projizierte Raster, das gleiche Ausmaß und den gleichen Wert bei (-100, 50) ergeben. Allerdings offenbar nicht :(
EDi
1
Die beiden Programme erstellen unterschiedliche Gitter, auf die man sich verziehen kann. Selbst wenn die bilineare Abtastung genau gleich wäre, befinden sich die zu interpolierenden Punkte an verschiedenen Stellen, und Sie hätten unterschiedliche Antworten. Die Herkunft und die Pixelgröße sind unterschiedlich. Sie können einige Flags in gdalwarp (-te, -tr usw.) setzen, um zu versuchen, die R-Version zu reproduzieren, und dann die Pixelwerte vergleichen und sehen, wie unterschiedlich sie sind.
Ich habe mehrfach festgestellt, dass die Verwendung des -orderFlags (die "Reihenfolge des zum Verziehen verwendeten Polynoms") gdalwarpauch ohne Verwendung von GCPs genauere Ergebnisse liefert.
Christoph

Antworten:

10

Schöne und reproduzierbare Frage. Persönlich würde ich erwarten, dass der Grund für den Unterschied in den Implementierungen der bilinearen Reprojektion liegt. Sie können natürlich den Quellcode für die beiden Ansätze untersuchen, aber ich würde erwarten, dass dies ein großer Overkill ist.
Es scheint, dass die R-Implementierung größere "Fehler" / "Änderungen" als die unformatierte GDAL-Version einführt (zumindest in meinen Versionen und Tests - projectRaster führt Änderungen um + -0,01 ein, während GDAL Werte um + -0,002 liefert).

Wenn Sie beide Ansätze mit einer Neuprojektion des nächsten Nachbarn vergleichen, stimmen sie wie erwartet überein.

Mikkel Lydholm Rasmussen
quelle
Danke für diesen Hinweis mit den Projektionsmethoden! Wenn ich Zeit finde, werde ich mir diese genauer ansehen (ich kenne R jedoch besser als C).
EDi