Suitable raster format for elevation data with sub-meter accuracy?

Question

I've got some elevation data represented in ascii gridded xyz-files which I want to convert to a conventional raster format. What's a bit special is that the elevation values have sub-meter accuracy, as they have one decimal (for example: '101.7')

As a rule I prefer to use open source software and stick to conventional, open, widely supported formats.

So the first thing I tried was to convert the data using gdal to geotiff, but specifying the data-type as Float32/64 distorted the elevation values somewhat. When I check the pixel-value in qgis it reports it as 101.70001220703125.

The command I use to convert an ascii file into raster is:

gdal_translate -ot Float32 -a_srs "EPSG:32633" -co "TILED=YES" -co "COMPRESS=LZW" -co "PREDICTOR=2"  -co "ZLEVEL=9" file.txt file.tif

Questions:
What would be suitable raster formats for sub-meter elevation data?
Is there any other software I should consider using instead of gdal?
Is this distortion of decimal precision perhaps inherently unavoidable? (Some doubts have been raised in the comments below though).

Edit:
Some more information about the original text format: I cannot disseminate this data and provide a real example, but the original format is rather trivial. The text files contain UTM coordinates for a specific zone, 33N for example. The spacing between the points are regular, with a spacing of 5 meters between them. Just as an example two adjacent XYZ values would be:
680000 5420000 101.7
680005 5420000 101.3

Actually, it should be possible to use just the above two lines and invoke gdal_translate as described above, it should be sufficient to discern how gdal handles the numerical precision.

Answer 1

Try the IEEE 754 Converter http://www.h-schmidt.net/FloatConverter/IEEE754.html. Read the note:

Rounding errors: Not every decimal number can be expressed exactly as a floating point number. This can be seen when entering "0.1" and examining its binary representation which is either slightly smaller or larger, depending on the last bit.

You can convert decimal presentation of 101.7 into binary presentation 01000010110010110110011001100110 but it is a little bit less than 101.7. Turn the last binary digit into 1 and you get 01000010110010110110011001100111 but now it is a little bit more than 101.7. In the 32-bit binary world there is nothing between those two values. You must just tolerate it or to use the method suggested by Mikkel Lydholm Rasmussen and multiply by ten and use integer data type.

EDIT:

@mkennedy correctly pointed out that the observed error was far too big compared with the rounding error calculated with IEEE754 converter:

Observed: 101.70001220703125
IEE754:   101.70000457763672

Thus the observed value added 0.0122 millimeters to the DEM height value while the rounding error from 32 bit adds only 0.0046 millimeters. This is theoretically interesting while perhaps the unexplained error of 0.0076 mm is not meaningful for most practical work.

I did some more debugging and I wrote this test XYZ file "float_error.xyz":

ncols         1
nrows         2
xllcorner     680000
yllcorner     5420000
cellsize      5.0
NODATA_value  -9999
101.7
101.3

Gdalinfo from this file is:

gdalinfo float_error.xyz -stats
Driver: AAIGrid/Arc/Info ASCII Grid
Files: float_error.xyz
Size is 1, 2
Coordinate System is `'
Origin = (680000.000000000000000,5420010.000000000000000)
Pixel Size = (5.000000000000000,-5.000000000000000)
Corner Coordinates:
Upper Left  (  680000.000, 5420010.000)
Lower Left  (  680000.000, 5420000.000)
Upper Right (  680005.000, 5420010.000)
Lower Right (  680005.000, 5420000.000)
Center      (  680002.500, 5420005.000)
Band 1 Block=1x1 Type=Float32, ColorInterp=Undefined
  Minimum=101.300, Maximum=101.700, Mean=101.500, StdDev=0.200
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=101.69999694824
    STATISTICS_MEAN=101.5
    STATISTICS_MINIMUM=101.30000305176
    STATISTICS_STDDEV=0.19999694824219

Next I converted XYZ into GeoTIFF as:

gdal_translate -of GTiff -ot Float32 float_error.xyz float_error.tif
Input file size is 1, 2
0...10...20...30...40...50...60...70...80...90...100 - done.

Finally new gdalinfo for finding out the conversion error:

gdalinfo float_error.tif
Driver: GTiff/GeoTIFF
Files: float_error.tif
Size is 1, 2
Coordinate System is `'
Origin = (680000.000000000000000,5420010.000000000000000)
Pixel Size = (5.000000000000000,-5.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  680000.000, 5420010.000)
Lower Left  (  680000.000, 5420000.000)
Upper Right (  680005.000, 5420010.000)
Lower Right (  680005.000, 5420000.000)
Center      (  680002.500, 5420005.000)
Band 1 Block=1x2 Type=Float32, ColorInterp=Gray
  Min=101.300 Max=101.700
  Minimum=101.300, Maximum=101.700, Mean=101.500, StdDev=0.200
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=101.69999694824
    STATISTICS_MEAN=101.5
    STATISTICS_MINIMUM=101.30000305176
    STATISTICS_STDDEV=0.19999694824219

My GDAL 2.0 seems to produce slightly different result, and it is an exact match with the result that the IEEE754 converter gives for binary value 01000010110010110110011001100110

Conclusion 1: GDAL 2.0 seems to make the conversion from ASCII Grid decimal values into 32 floating points as accurately as possible.

Conclusion 2: The Vincent's value 101.70001220703125 is still unexplained. The GDAL version may be different. I suppose that Click with QGIS method is accurate and does not involve any interpolation.

Answer 2

Projection matters more than format. If you are not getting enough precision, switch from a UTM projection to something centered more around you area. And if you try float64, you'll do better, but my guess it that you are then shaving yaqs for your use case. If you really want to preserve exact numbers, scale elevation and use ints.