Perbedaan antara gdalwarp dan projectRaster

9

Saya mencoba memproyeksikan Raster. Di R ada projectRaster()fungsi untuk ini (di bawah contoh yang sepenuhnya direproduksi):

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

Itu bekerja dengan baik. Namun itu cukup lambat.

Untuk menambah kecepatan, saya menggunakan gdalwarp(dengan SSD, biaya membaca dan menulis dari / ke disk / R tidak terlalu tinggi).

Namun, saya tidak dapat mereproduksi hasil projectRaster()menggunakan 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)

Tampaknya berfungsi, namun hasilnya berbeda:

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

Apa yang saya lewatkan / lakukan salah?

Apakah ada alternatif lain (lebih cepat) projectRaster()?

raster gdal r coordinate-system gdalwarp EDI
sumber
Tidak ada Saya memberikan contoh yang sepenuhnya direproduksi (harus bekerja dengan Linux atau Mac) ...
EDi
Apa yang kamu harapkan Apakah kedua opsi menggunakan proj.4 yang sama?
Saya berharap bahwa kedua metode menghasilkan raster yang diproyeksikan ulang yang sama, tingkat yang sama dan nilai yang sama di (-100, 50). Namun, mereka tampaknya tidak :(
EDi
1
Kedua program ini menciptakan grid yang berbeda untuk digunakan. Sekalipun pengambilan sampel bilinear persis sama, poin yang sedang diinterpolasi ada di tempat yang berbeda, dan Anda akan memiliki jawaban yang berbeda. Asal dan ukuran piksel berbeda. Anda dapat mengatur beberapa flag di gdalwarp (-te, -tr, dll.) Untuk mencoba dan mereproduksi versi R, lalu membandingkan nilai piksel dan melihat betapa berbedanya mereka.
Saya menemukan pada beberapa kesempatan bahwa menggunakan -orderbendera ("urutan polinominal yang digunakan untuk melengkungkan") gdalwarpbahkan tanpa menggunakan GCP menghasilkan hasil yang lebih akurat.
christoph

Jawaban:

10

Pertanyaan yang bagus dan dapat direproduksi. Secara pribadi, saya berharap bahwa alasan perbedaannya adalah dalam implementasi proyeksi bilinear. Anda jelas dapat melihat ke dalam kode sumber untuk dua pendekatan, tetapi saya berharap itu menjadi kerja keras yang sangat besar.
Tampaknya implementasi R memperkenalkan "kesalahan" / "perubahan" yang lebih besar daripada versi GDAL mentah (setidaknya dalam versi & tes saya - projectRaster memperkenalkan perubahan sekitar + -0,01 sementara GDAL memberikan nilai sekitar + -0,002).

Jika Anda membandingkan kedua pendekatan menggunakan proyeksi tetangga terdekat, mereka cocok dengan yang diharapkan.

Mikkel Lydholm Rasmussen
sumber
Terima kasih atas petunjuk ini dengan metode proyeksi! Jika saya menemukan waktu saya akan melihat lebih dalam ke dalamnya (Namun, saya lebih akrab dengan R kemudian dengan C).
EDI