Estoy tratando de proyectar una trama. En R existe la projectRaster()
función para esto (debajo de un ejemplo completamente reproducible):
# 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')
Que funciona bien. Sin embargo, es bastante lento.
Para aumentar la velocidad, pensé usarlo gdalwarp
en su lugar (con un SSD, el costo de leer y escribir desde / hacia el disco / R no es muy alto).
Sin embargo, no puedo reproducir los resultados del projectRaster()
uso 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)
Parece funcionar, sin embargo, los resultados son diferentes:
# 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é me estoy perdiendo / haciendo mal?
¿Hay otras alternativas (más rápidas) para projectRaster()
?
-order
bandera (el "orden del polinominal usado para la deformación")gdalwarp
incluso sin usar GCP produjo resultados más precisos.Respuestas:
Pregunta agradable y reproducible. Personalmente, esperaría que la razón de la diferencia esté en las implementaciones de la reproyección bilineal. Obviamente, puede buscar en el código fuente los dos enfoques, pero espero que sea una gran exageración.
Parece que la implementación de R introduce "errores" / "cambios" más grandes que la versión sin procesar de GDAL (al menos en mis versiones y pruebas - projectRaster introduce cambios alrededor de + -0.01 mientras que GDAL da valores alrededor de + -0.002).
Si compara ambos enfoques utilizando una reproyección de vecino más cercano, coinciden como se esperaba.
fuente