startingvalue <- 10000
precision <- 0.01
rfactor <- 256*256 * precision
gfactor <- 256 * precision
r <- floor((startingvalue +dtm)*(1/precision) / 256 / 256)
g <- floor((startingvalue +dtm - r*rfactor)*(1/precision) / 256)
b <- floor((startingvalue +dtm - r*rfactor - g*gfactor)*(1/precision))
rgb_dem <- c(r,g,b)
writeRaster(rgb_dem, "dtm_rgb.tif"), datatype="INT1U", NAflag=NA)Orthophoto imagery
Original data
Original data were RGB-Orthoimages which were present as 93 uncompressed image tiles with corresponding overviews.
The tiles had following properties:
- Bands: 3
- Width: 15000 px
- Height: 15000 px
- Data type: Byte - Eight bit unsigned integer
Preprocessing
For further processing, the tiles were mosaiced into one virtual raster (vrt) with the dimensions 165000 x 180000 pixel.
GDAL command: gdalbuildvrt mosaic.vrt tiles\*.tif
Compression
The images were merged in mosaics (via the virtual raster file) and stored as single files with varying compression types and settings. We differentiate between images with internal overviews and those without overviews. Also it is important to distinguish between lossless compression and lossy compression.
With Overviews
For a smooth visualization experience it is necessary to have overviews1 of the image at different zoom levels. Cloud optimized GeoTIFFs (COG) include them by default.
| format | compression | lossless | size (GB) | size (%) | comment |
|---|---|---|---|---|---|
| GeoTIFF-tiles + overviews | none | yes | 72 | 100 | original |
| COG | LZW | yes | 68 | 94 | fast to write slow to read |
| COG | DEFLATE | yes | 55 | 76 | fast to read |
| COG | DEFLATE | yes | 45 | 63 | Predictor=YES |
| COG | ZSTD | yes | 55 | 76 | higher compression speed than Deflate |
| COG | LZMA | yes | 45 | 63 | very slow to write |
| COG | LERC | yes | 43 | 60 | MAX_Z_ERROR=0 |
| COG | LERC_ZSTD | yes | 43 | 60 | |
| COG | WEBP (100%) | yes | 34 | 47 | lossless WebP, very slow to write |
| COG | JPEG (75%) | no | 4 | 6 | YCbCr-Colorspace |
| COG | JPEG (75%) | no | 4 | 6 | YCbCr-Colorspace, NoData |
| COG | WEBP (75%) | no | 3 | 4 |
Lossless as well as lossy compression worked very well with the WEBP image format. The downside is that not all browsers and programs support the format (see here).
Without Overviews
If data should be stored just to archive it and for some analyzing workflows it might not be necessary to have overviews. To compare the size without overviews some regular mosaics were generated.
| format | compression | lossless | size (GB) | size (%) | comment |
|---|---|---|---|---|---|
| GeoTIFF-tiles | none | yes | 60 | 100 | original |
| GeoTIFF-mosaic | LZW | yes | 50 | 83 | |
| GeoTIFF-mosaic | JPEG (75%) | no | 8 | 13 | |
| GeoTIFF-mosaic | JPEG (75%) | no | 3 | 5 | YCbCr-Colorspace |
ERDAS format
Erdas Imagine HFA format consists of image files in .img format, overviews in .rrd format and additional .xml files. For size comparison a single ortho-image tile (12257 x 12257 px) was converted to Geotiff with internal overviews (COG).
| format | compression | lossless | size (MB) | size (%) | comment |
|---|---|---|---|---|---|
| HFA – Erdas Imagine .img + .rrd | none | yes | 576 + 192 | 100 | original |
| COG | DEFLATE | yes | 432 | 56 | Predictor=YES |
Elevation data
A completely different type of data is elevation data (e.g. digital terrain or canopy models). While it can also be stored as raster data its structure is completely different from photographic data such as above orthophoto imagery. Instead of three (or more) bands its data is usually located in a single band (just one value per location). But this single value can be much more differentiated (floating data instead of integer). For this reasons compression algorithms have a different performance on elevation data than they do for photographic data.
Original data
The image had following properties:
- Bands: 1
- Width: 4000
- Height: 4000
- Data type: Float32 - Thirty two bit floating point
Compression
One-band Elevation
| format | compression | lossless | size (MB) | size (%) | comment |
|---|---|---|---|---|---|
| GeoTIFF | none | yes | 63 | 100 | original |
| COG | LERC | yes | 83 | 135 | MAX_Z_ERROR=0 |
| COG | LZW | yes | 61 | 94 | fast to write slow to read |
| COG | DEFLATE | yes | 46 | 71 | fast to read |
| COG | LERC_DEFLATE | yes | 46 | 71 | |
| COG | ZSTD | yes | 44 | 68 | higher compression speed than Deflate |
| COG | LERC_ZSTD | yes | 44 | 68 | |
| COG | LZW | yes | 43 | 66 | Predictor=YES (predictor 3 for floating point) |
| COG | DEFLATE | yes | 35 | 54 | Predictor=YES (predictor 3 for floating point) |
| COG | ZSTD | yes | 35 | 54 | Predictor=FLOATING_POINT, LEVEL=1 |
| COG | DEFLATE | yes | 34 | 53 | Predictor=STANDARD |
| COG | ZSTD | yes | 34 | 53 | Predictor=FLOATING_POINT |
| COG | ZSTD | yes | 33 | 52 | Predictor=STANDARD, LEVEL=22 |
| COG | ZSTD | yes | 31 | 49 | Predictor=STANDARD |
| COG | LZMA | yes | 27 | 42 | very slow to write |
| COG | LERC_ZSTD | no | 12 | 18 | MAX_Z_ERROR=0.01 (precision of 1 cm) |
| COG | LERC_ZSTD | no | 8 | 12 | MAX_Z_ERROR=0.025 (precision of 2.5 cm) |
| COG | JPEG (75%) | no | only possible with 8-bit data | ||
| COG | JPEG (75%) | no | only possible with 8-bit data | ||
| COG | WEBP (75%) | no | only possible with 3 or 4 band data |
RGB-encoded
Another possibility to store the data is to encode the elevation in 3 channels. This has the advantage that the 3 bands can be of datatype byte (0-255) which reduces file size. This technique is often used for web mapping, where loaded tiles should be as small as possible to reduce loading times. Mapbox, Maptiler etc. use it with a precision factor of 0.1 which means that elevation is stored to the decimeter (e.g. 15.3 m). The formula to decode RGB-values to elevation is then height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1). A detailed overview about possibilities and limitations of this technique was written by Frédéric Rodrigo in this article.
Note
With this encoding formula there are certain limitations which data can be encoded. A starting value of -10000 for example means that the lowest elevation on earth (Mariana Trench), which is approximately 10984 m deep, can not be encoded in this format. Also evelation values with higher precision (e.g. in the centimeter range) can not precisely encoded with a factor of 0.1.
Encoding
Here we encoded the elevation with higher precision in R using the following script:
Decoding
It can be decoded with the formula elevation = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.01).
Bit-reduction
To reduce the file size even more we can reduce the bitdepth of the B-channel. Orignial 8-bit data encodes the elevation in 1cm steps. If we remove the last bit, data is stored in 128 steps, reducing the precision to 2cm steps, and if we go further removing the last 2 bits results in 4cm steps.
If data is visualized as Multidirectional Hillshade (small differences are expected to be visible earlier here than if height is visualized) the artificially introduced errors seem to be visible from approximately this scales (it is more clear in flat areas):
- 7bits: greater than 1:250
- 6bits: greater than 1:1.000
- 5bits: greater than 1:2.000
| format | compression | lossless | size (MB) | size (%) | comment |
|---|---|---|---|---|---|
| GeoTIFF | none | yes | 65 | 100 | original |
| GeoTIFF | none | (no) | 31 | 48 | RGB encoded |
| COG | LZMA | (yes) | 17 | 26 | RGB encoded as COG |
| COG | ZSTD | (yes) | 14 | 22 | RGB encoded as COG, Predictor=STANDARD |
| COG | WEBP | (yes) | 8 | 12 | RGB encoded as COG WEBP |
| GeoTIFF | WEBP | (yes) | 5 | 8 | RGB encoded as WEBP with 7bits for B-band |
| GeoTIFF | WEBP | (yes) | 4 | 6 | RGB encoded as WEBP with 6bits for B-band |
| GeoTIFF | WEBP | (yes) | 3 | 5 | RGB encoded as WEBP with 5bits for B-band |
Hillshade
If data is only used for visualization it might be an alternative to use a hillshade, this can also dramatically reduce file size.
| format | compression | lossless | size (MB) | size (%) | comment |
|---|---|---|---|---|---|
| GeoTIFF | none | yes | 65 | 100 | original |
| GeoTIFF | JPEG | no | 2 | 3 | hillshade JPEG 75% |
Conclusion
We have seen that a multitude of compression algorithms and settings exists which all have different outcomes, may it be in compression ratio, computation time or data precision. Thus, there is no single algorithm which fits all purposes, but there are barely any reasons to not use any compression. WEBP-compression showed impressive results for lossless and lossy compression of photographic imagery, the downside may be that it is not yet adopted by all programs. Otherwise, LERC-algorithm had good compression ratio for lossless compression and JPEG for lossy compression. Since lossy compression alters the original data slightly it may not be used when raw data is definitely needed (e.g. remote sensing analytics). For applications where vizualisation is the main goal (e.g. background imagery) lossy compression might be the better choice as the visual appearance remains the same but data compression is huge.
For elevation data the good compression ratios are achieved when the data is encoded in RGB-values and stored as type byte. It can be even enhanced if a byte-reduction is applied. However this needs dedicated encoding/decoding steps and might not be trivial to use for all purposes. If compression is not the highest goal there are good compression algorithms available for single-band (standard) elevation data.
But of course these tests may not be representitive for all geospatial raster data sources and purposes.
Additional Sources
Footnotes
Overviews:
lower-resolution, downsampled versions of the image, which are precomputed and stored with the image. Each overview level halves the height and width pixels, so each additional level reduces the resolution by a factor of 2. These overviews enable fast reads at a range of resolutions by minimizing the amount of bytes that need to be read.
source: Kyle Barron↩︎
Citation
For attribution, please cite this work as:
Wiesehahn, Jens. 2022. “GeoTiff Image Compression.” April
9, 2022. https://wiesehahn.github.io/posts/image_compression/.