0

When writing then reading a SpatRaster with certain file formats (ESRI .hdr labeling format, netCDF and ENVI, but not GeoTIFF or Erdas HFA), the raster origin, resolution, and coordinates change slightly. Enough for equality tests for these parameters to return the FALSE value. Here is an example, with GeoTiff then 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

Why do such shifts occur? Is there a way to prevent them, and if not, what can I do to correct them after reading?

blackgreen
  • 34,072
  • 23
  • 111
  • 129
Jean-Luc Dupouey
  • 341
  • 1
  • 6

1 Answers1

0

The shifts occur because an insufficient number of decimals is used to store a fraction, aggravated by a faulty approach to storing the extent or a raster. The FLT file uses "origin / resolution" approach to georeferencing

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"

As you can see, there is ULXMAP and ULYMAP (upper left x and y coordinates of the center of the cell), but the coordinates of the other corner are not stored. So the must be computed.

In this case the min (left) x-coordinate is

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

And the max (right) x-coordinate is

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

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

Where it should be

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

A very small difference, but the loss of precision can easily be avoided by storing the max x and y coordinates instead of the x and y resolution. Unfortunately that is not wat is done here, and even GDAL used this (I think flawed) approach internally.

I think the difference between filetypes you see is due to the numbers of decimals that is stored (double vs single precision).

If this were to cause problems (normally such small differences are ignored) you could fix this by assigning the correct extent

 ext(r3) <- c(-180, 180, -90, 90)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63