0

在写入然后读取具有某些文件格式(ESRI .hdr 标签格式、netCDF 和 ENVI,但不是 GeoTIFF 或 Erdas HFA)的 SpatRaster 时,光栅原点、分辨率和坐标略有变化。这些参数的相等性测试足以返回 FALSE 值。这是一个例子,GeoTiff 然后是 flt。

Nrow=45; Ncol=108
r1 <- rast(nrow=Nrow,ncol=Ncol,vals=1:(Nrow*Ncol),ext=c(-180,180,-60,90))
origin(r1)
# [1] 0 0

writeRaster(r1,"test.tif",overwrite=TRUE)

r2 <- rast("test.tif")
origin(r2)  
# origin is OK, same as r1
# [1] 0 0
 
res(r1)==res(r2)
# [1] TRUE TRUE

writeRaster(r1,"test.flt",overwrite=TRUE)
r3 <- rast("test.flt")
origin(r3)  
# origin has been shifted
# [1] 1.705303e-13 5.684342e-14

res(r1)==res(r3)   
# the lag is large enough for r3 to be considered as a different raster from r1
# [1] FALSE FALSE

这里有三个问题:

. 为什么会发生这样的转变?

. 有没有办法阻止它们?

. 如果没有,阅读后纠正它们的最佳方法是什么?

谢谢。

4

1 回答 1

0

发生偏移的原因是用于存储分数的小数位数不足,存储范围或栅格的错误方法加剧了这种情况。该FLT文件使用“原点/分辨率”方法进行地理配准

readLines("test.hdr")
 [1] "BYTEORDER      I"                 "LAYOUT         BIL"               "NROWS          45"               
 [4] "NCOLS          108"               "NBANDS         1"                 "NBITS          32"               
 [7] "BANDROWBYTES   432"               "TOTALROWBYTES  432"               "PIXELTYPE      FLOAT"            
[10] "ULXMAP         -178.333333333333" "ULYMAP         88.3333333333333"  "XDIM           3.33333333333333" 
[13] "YDIM           3.33333333333333"  "NODATA         1.#QNAN"    

可以看到,有ULXMAPand ULYMAP(单元格中心的左上x和y坐标),但是没有存储另一个角的坐标。所以必须计算。

在这种情况下,最小(左)x 坐标是

XDIM <- 3.33333333333333
xmin <- -178.333333333333 - 0.5 * XDIM
print(xmin, 16)
# -179.9999999999997

并且最大(右)x坐标是

NCOLS <- 108
xmax <- xmin + NCOLS * XDIM
print(xmax, 16)
#180

yres = (xmax - xmin) / NCOLS
print(yres, 17)
[1] 3.3333333333333304

它应该在哪里

print(360/NCOLS, 17)
[1] 3.3333333333333335

一个非常小的差异,但可以通过存储最大 x 和 y 坐标而不是 x 和 y 分辨率来轻松避免精度损失。不幸的是,这不是在这里完成的,甚至 GDAL 在内部也使用了这种(我认为有缺陷的)方法。

我认为您看到的文件类型之间的差异是由于存储的小数位数(双精度与单精度)。

如果这会导致问题(通常会忽略如此小的差异),您可以通过分配正确的范围来解决此问题

 ext(r3) <- c(-180, 180, -90, 90)
于 2021-04-10T05:11:02.797 回答