Różnica między gdalwarp a projectRaster

9

Próbuję wyświetlić raster. W R jest do tego projectRaster()funkcja (poniżej w pełni odtwarzalnego przykładu):

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

Które działa dobrze. Jest jednak dość wolny.

Aby zwiększyć prędkość, gdalwarpzdecydowałem się zamiast tego użyć (z dyskiem SSD koszt odczytu i zapisu z / na dysk / R nie jest bardzo wysoki).

Nie mogę jednak odtworzyć wyników projectRaster()użycia 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)

Wydaje się, że działa, ale wyniki są inne:

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

Czego mi brakuje / robię źle?

Czy istnieją inne (szybsze) alternatywy projectRaster()?

raster gdal r coordinate-system gdalwarp EDi
źródło
Nikt?
Podałem
Czego oczekujesz Czy obie opcje używają tego samego proj.4?
Oczekuję, że obie metody dadzą ten sam ponownie rzutowany raster, ten sam zakres i tę samą wartość przy (-100, 50). Jednak najwyraźniej nie :(
EDi
1
Oba programy tworzą różne siatki, na których można się wypaczać. Nawet jeśli próbkowanie dwuliniowe było dokładnie takie samo, interpolowane punkty znajdują się w różnych miejscach i można uzyskać różne odpowiedzi. Początki i rozmiary pikseli są różne. Możesz ustawić kilka flag w gdalwarp (-te, -tr itp.), Aby spróbować odtworzyć wersję R, a następnie porównać wartości pikseli i zobaczyć, jak się różnią.
Wielokrotnie stwierdziłem, że użycie -orderflagi („kolejność wielomianu używanego do wypaczania”) gdalwarpnawet bez użycia GCP dało bardziej dokładne wyniki.
christoph

Odpowiedzi:

10

Ładne i powtarzalne pytanie. Osobiście oczekiwałbym, że przyczyną tej różnicy są implementacje dwuliniowej reprojektcji. Możesz oczywiście zajrzeć do kodu źródłowego dla tych dwóch podejść, ale spodziewałbym się, że będzie to ogromna przesada.
Wygląda na to, że implementacja R wprowadza większe „błędy” / „zmiany” niż surowa wersja GDAL (przynajmniej w moich wersjach i testach - projectRaster wprowadza zmiany około + -0,01, podczas gdy GDAL podaje wartości około + -0,002).

Jeśli porównasz oba podejścia przy użyciu reprojection najbliższego sąsiada, pasują one zgodnie z oczekiwaniami.

Mikkel Lydholm Rasmussen
źródło
Dzięki za tę wskazówkę dotyczącą metod projekcji! Jeśli znajdę czas, przyjrzę się im głębiej (jednak bardziej znam R niż C).
EDi