Diferença entre gdalwarp e projectRaster

9

Estou tentando projetar uma varredura. Em R existe a projectRaster()função para isso (abaixo de um exemplo totalmente reprodutível):

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

O que funciona bem. No entanto, é bastante lento.

A fim de aumentar a velocidade que eu gostaria de usar gdalwarp(com um SSD, o custo de leitura e gravação de / para disco / R não é muito alto).

No entanto, não consigo reproduzir os resultados do 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, no entanto, os resultados são 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)

O que estou perdendo / fazendo errado?

Existem outras alternativas (mais rápidas) para projectRaster()?

raster gdal r coordinate-system gdalwarp EDi
fonte
Ninguém? I fornecido um exemplo totalmente reproduzível (deve trabalhar com Linux ou Mac) ...
EDi
O que você está esperando? As duas opções usam o mesmo proj.4?
Espero que ambos os métodos produzam a mesma rasterização reprojetada, a mesma extensão e o mesmo valor em (-100, 50). No entanto, eles aparentemente não o fazem :(
EDi
11
Os dois programas estão criando grades diferentes para serem distorcidas. Mesmo que a amostragem bilinear seja exatamente a mesma, os pontos que estão sendo interpolados estão em lugares diferentes e você terá respostas diferentes. As origens e os tamanhos de pixel são diferentes. Você pode definir alguns sinalizadores em gdalwarp (-te, -tr, etc.) para tentar reproduzir a versão R e, em seguida, comparar os valores de pixel e ver quão diferentes eles são.
Eu descobri em várias ocasiões que o uso da -orderflag (a "ordem do polinômio usado para deformar") gdalwarpmesmo sem usar GCPs produzia resultados mais precisos.
christoph

Respostas:

10

Pergunta agradável e reproduzível. Pessoalmente, eu esperaria que o motivo da diferença esteja nas implementações da reprojeção bilinear. Obviamente, você pode procurar no código fonte as duas abordagens, mas eu esperaria que isso fosse um grande exagero.
Parece que a implementação do R introduz "erros" / "alterações" maiores que a versão bruta do GDAL (pelo menos em minhas versões e testes - o projectRaster introduz alterações em torno de + -0,01 enquanto o GDAL fornece valores em torno de + -0,002).

Se você comparar as duas abordagens usando uma reprojeção de vizinho mais próximo, elas corresponderão conforme o esperado.

Mikkel Lydholm Rasmussen
fonte
Obrigado por esta dica com os métodos de projeção! Se eu encontrar tempo, examinarei mais profundamente essas questões (no entanto, estou mais familiarizado com R do que com C).
Edi 21/07