Rounding issue with R, raster, and lat-lon-data

145 views Asked by At

I am importing and exporting ASTER digital elevation models in R using the raster package. These data have a geographic coordinate system, so that the extent of the raster and its cellsize are in degrees.

I find that there must be some kind of arithmetic underflow when importing and exporting one and the same dataset. For example, I import a toy ASTER scene ASTGTMV003_N48W122_dem.tif into R (I uploaded the data set here) and want to find out its resolution:

# example scene downloaded from EarthExplorer (https://earthexplorer.usgs.gov/)
r <- raster("ASTGTMV003_N48W122_dem.tif") 
res(r)

And R tells me that the resolution (x, y) is: 0.0002777778 0.0002777778 QGIS, however, suggests that the same dataset has much more digits: 0.0002777777777777779944,-0.0002777777777777779944.

So here comes the twist: I now export this dataset as a GeoTIFF and import it again and keep all parameters unchanged. R would tell me that nothing has changed:

writeRaster(r, "test.tif") # export GTiff
r2 <- raster("test.tif")   # import again
res(r) == res(r2)          # which returns TRUE TRUE

Yet in GIS, I find that my new test.tif has a resolution of 0.0002777777778394882884,-0.0002777777778394882884, which is different from the original ASTER scene. While R seems to be forgiving in this regard, GIS tells me that the extents and the cellsizes are NOT the same. I reproduced this example with other ASTER and SRTM scenes and the problem persists. Yet the digits are always different from the 13th digit onward.

I speculate that the problem is tied to the uneven number of rows and columns (3601, 3601). At least, using ALOS3D 30m data with 3600 rows and columns, I couldn't reproduce this issue.

Did anyone face a similar issue with extents and resolutions when exporting and importing SRTM or ASTER data?

0

There are 0 answers