GeoTiff Image Compression

Introduction

Geospatial raster data can quickly balloon to enormous sizes, creating challenges for storage and distribution. But with dozens of compression algorithms available—each with different trade-offs between file size, speed, and data fidelity—choosing the right approach isn’t straightforward. We tested various compression formats available for GeoTIFF images via GDAL on two common data types: RGB orthoimagery and digital elevation models. This post compares traditional methods with newer approaches (WEBP, LERC, RGB-encoded elevation) to help you select appropriate compression for your use case.

Orthophoto imagery

Original data

As test data we used five ortho-image tiles distributed in heterogeneous landscapes of Lower Saxony, Germany. They are available as Cloud Optimized Geotiff from the State Office for Geoinformation and Surveying of Lower Saxony (LGLN) via s3 storage. Each of them covers an area of 2x2 km with a resolution of 20 cm (4000x4000 px) and comes with 4 spectral bands (Red, Green, Blue and Infrared channel) uncompressed in data format Byte (eight bit unsigned integer).

# GDAL commands via OSGeo4W Shell

>gdalinfo /vsicurl/https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/324905842/2024-09-05/dop20rgbi_32_490_5842_2_ni_2024-09-05.tif
Driver: GTiff/GeoTIFF
Files: /vsicurl/https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/324905842/2024-09-05/dop20rgbi_32_490_5842_2_ni_2024-09-05.tif
Size is 10000, 10000
Coordinate System is:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            MEMBER["European Terrestrial Reference Frame 2020"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6┬░E and 12┬░E: Austria; Denmark - onshore and offshore; Germany - onshore and offshore; Italy - onshore and offshore; Norway including Svalbard - onshore and offshore; Spain - offshore."],
        BBOX[36.53,6,84.01,12.01]],
    USAGE[
        SCOPE["Pan-European conformal mapping at scales larger than 1:500,000."],
        AREA["Europe between 6┬░E and 12┬░E and approximately 36┬░30'N to 84┬░N."],
        BBOX[36.53,6,84.01,12.01]],
    ID["EPSG",25832]]
Data axis to CRS axis mapping: 1,2
Origin = (490000.000000000000000,5844000.000000000000000)
Pixel Size = (0.200000000000000,-0.200000000000000)
Metadata:
  OVR_RESAMPLING_ALG=BILINEAR
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  490000.000, 5844000.000) (  8d51' 6.70"E, 52d44'44.76"N)
Lower Left  (  490000.000, 5842000.000) (  8d51' 6.92"E, 52d43'40.03"N)
Upper Right (  492000.000, 5844000.000) (  8d52'53.36"E, 52d44'44.88"N)
Lower Right (  492000.000, 5842000.000) (  8d52'53.53"E, 52d43'40.15"N)
Center      (  491000.000, 5843000.000) (  8d52' 0.13"E, 52d44'12.46"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  NoData Value=0
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  NoData Value=0
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  NoData Value=0
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
Band 4 Block=512x512 Type=Byte, ColorInterp=Undefined
  NoData Value=0
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313

Compression

The images were converted and stored on disk with various compression types and settings via GDAL. It is important to distinguish between lossless compression and lossy compression. While lossless compression ensures that tha underlying image data stays entirely unchanged this is not the case with lossy compression. Lossy compression can nevertheless ensure, that images stay visually unchanged (only minor changes which are not perceived by humans).

# pak::pak("USDAForestService/gdalraster")
library(gdalraster)
library(purrr)
library(fs)

# Set gdal configurations (for reproducability?)
# Multi-threaded (comment GDAL_NUM_THREADS to use default single threaded
set_config_option("GDAL_NUM_THREADS", "32")
set_config_option("GDAL_CACHEMAX", "8000")
set_config_option("OVERVIEWS", "IGNORE_EXISTING")

# get orthoimages from LGLN open data
urls <- c(
  "https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/324905842/2024-09-05/dop20rgbi_32_490_5842_2_ni_2024-09-05.tif",
  "https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/326085740/2022-05-09/dop20rgbi_32_608_5740_2_ni_2022-05-09.tif",
  "https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/326085732/2022-05-09/dop20rgbi_32_608_5732_2_ni_2022-05-09.tif",
  "https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/326125846/2024-09-21/dop20rgbi_32_612_5846_2_ni_2024-09-21.tif",
  "https://dop20-rgbi.s3.eu-de.cloud-object-storage.appdomain.cloud/326045852/2024-09-21/dop20rgbi_32_604_5852_2_ni_2024-09-21.tif"
)

dsns <- paste0("/vsicurl/", urls) # prefix for virtual file source



# Create a function to measure performance for each datasource and compression
raster_compression <- function(dsns, options) {
  results <- list()

  for (datasourcename in dsns) {
    for (option in options) {
      output_file <- file_temp(ext = ".tif")

      write_ds <- function() {
        createCopy(
          src_filename = datasourcename,
          dst_filename = output_file,
          format = "COG",
          options = option$setting
        )
      }

      # Get write time
      writetime <- system.time(write_ds())["elapsed"]

      # Get file size in MB
      data_size_uncompressed <- vsi_stat(datasourcename, "size") / (1024^2)
      data_size_compressed <- vsi_stat(output_file, "size") / (1024^2)

      # Get read time
      img <- new(GDALRaster, output_file)
      readtime <- system.time(read_ds(img))["elapsed"]

      # Store results
      results[[length(results) + 1]] <- list(
        datasource = basename(datasourcename),
        settings = option$naming,
        uncompressed_mb = data_size_uncompressed,
        compressed_mb = data_size_compressed,
        writetime_sec = writetime,
        readtime_sec = readtime
      )
    }
  }

  # Convert list to data frame
  results_df <- do.call(rbind, lapply(results, data.frame))
  return(results_df)
}



# create list of options
options <- list(
  # lossless
  list(setting = c("COMPRESS=NONE")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=YES")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=NO")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=1")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=22")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=22")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100")),
  list(setting = c("COMPRESS=LERC", "MAX_Z_ERROR=0")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=6")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=22")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=1")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=5")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=9")),
  # lossy overviews
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "OVERVIEW_QUALITY=75")),
  # lossy
  list(setting = c("COMPRESS=WEBP", "QUALITY=75"))
  # only available in certain gdal builds
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=NO", "JXL_EFFORT=1")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=NO", "JXL_EFFORT=5")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=NO", "JXL_EFFORT=9")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=NO", "JXL_EFFORT=5", "JXL_DISTANCE=2"))
)

for (i in 1:length(options)) {
  options[[i]]$naming <- paste(options[[i]]$setting, collapse = ", ")
}


# apply function on list
benchmark_results <- raster_compression(dsns, options)

Comparison of Compression Performances
Size and processing time (mean ± standard deviation) depending on compression settings
Settings1	Size		Time		Speed
	Absolute (MB)	Relative (%)	Write (s)	Read (s)	Write (MB/s)	Read (MB/s)
COMPRESS=NONE	539 ± 0	100 ± 0	44.7 ± 5.7	5.1 ± 0.3	12.2 ± 1.4	105.6 ± 5.8
COMPRESS=LZW, PREDICTOR=NO	512 ± 24	95 ± 4	41.4 ± 2.7	4.9 ± 0.2	13.1 ± 0.9	104.7 ± 8.7
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=1	429 ± 6	80 ± 1	46.4 ± 6.9	4.7 ± 0.2	11.8 ± 1.8	91.9 ± 5.4
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=9	408 ± 20	76 ± 4	47.9 ± 5.8	4.8 ± 0.2	11.4 ± 1.4	85.7 ± 8.0
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=1	407 ± 18	76 ± 3	45.5 ± 1.9	4.9 ± 0.2	11.9 ± 0.5	84.0 ± 7.5
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=6	402 ± 20	75 ± 4	48.2 ± 6.8	4.8 ± 0.2	11.4 ± 1.6	83.8 ± 7.6
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=9	390 ± 22	72 ± 4	49.3 ± 5.0	4.8 ± 0.2	11.0 ± 1.0	82.3 ± 8.3
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=22	378 ± 25	70 ± 5	71.5 ± 2.6	4.7 ± 0.3	7.5 ± 0.3	80.2 ± 10.0
COMPRESS=LZW, PREDICTOR=YES	368 ± 26	68 ± 5	42.6 ± 4.8	4.9 ± 0.4	12.8 ± 1.4	76.3 ± 11.0
COMPRESS=LERC, MAX_Z_ERROR=0	335 ± 26	62 ± 5	46.0 ± 11.4	4.9 ± 0.2	12.2 ± 2.5	68.8 ± 8.7
COMPRESS=LZMA, LEVEL=1	334 ± 24	62 ± 4	48.2 ± 15.1	5.3 ± 0.2	11.9 ± 2.8	63.8 ± 6.7
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=1	333 ± 27	62 ± 5	40.8 ± 2.8	4.8 ± 0.3	13.3 ± 0.9	70.0 ± 9.6
COMPRESS=LERC_ZSTD, LEVEL=1	332 ± 27	62 ± 5	48.1 ± 5.4	4.8 ± 0.2	11.3 ± 1.2	69.1 ± 8.8
COMPRESS=LERC_ZSTD, LEVEL=9	332 ± 27	62 ± 5	49.2 ± 8.5	4.9 ± 0.1	11.2 ± 1.9	68.1 ± 7.0
COMPRESS=LERC_ZSTD, LEVEL=22	332 ± 27	62 ± 5	47.3 ± 5.7	5.0 ± 0.2	11.5 ± 1.4	66.1 ± 8.1
COMPRESS=LERC_DEFLATE, LEVEL=1	332 ± 27	62 ± 5	48.6 ± 7.1	4.8 ± 0.3	11.3 ± 1.5	68.9 ± 9.6
COMPRESS=LERC_DEFLATE, LEVEL=6	331 ± 27	62 ± 5	45.6 ± 5.7	4.9 ± 0.3	12.0 ± 1.6	68.4 ± 9.7
COMPRESS=LERC_DEFLATE, LEVEL=9	331 ± 27	61 ± 5	46.1 ± 4.4	5.0 ± 0.2	11.8 ± 1.1	66.8 ± 8.3
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=1	321 ± 20	60 ± 4	45.7 ± 7.2	4.9 ± 0.2	12.0 ± 1.9	65.8 ± 6.2
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=6	320 ± 21	59 ± 4	43.5 ± 9.3	4.9 ± 0.2	12.8 ± 2.3	66.1 ± 7.0
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=9	318 ± 24	59 ± 4	46.3 ± 3.6	4.9 ± 0.2	11.7 ± 0.9	64.5 ± 7.9
COMPRESS=LZMA, LEVEL=9	316 ± 23	59 ± 4	44.5 ± 3.2	5.2 ± 0.2	12.2 ± 0.9	61.3 ± 6.9
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=9	313 ± 19	58 ± 4	46.9 ± 10.3	4.8 ± 0.2	11.9 ± 2.2	65.8 ± 7.1
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=22	301 ± 15	56 ± 3	73.9 ± 2.1	4.7 ± 0.2	7.3 ± 0.2	64.3 ± 6.5
COMPRESS=WEBP, QUALITY=100	224 ± 17	42 ± 3	45.0 ± 7.2	4.8 ± 0.2	12.2 ± 1.7	47.2 ± 5.9
COMPRESS=WEBP, QUALITY=100, OVERVIEW_QUALITY=75	212 ± 15	39 ± 3	50.2 ± 5.9	4.8 ± 0.2	10.9 ± 1.3	43.9 ± 5.1
COMPRESS=WEBP, QUALITY=75	89 ± 12	16 ± 2	45.8 ± 8.2	4.8 ± 0.2	12.0 ± 2.0	18.4 ± 3.1
1 WEBP compression with QUALITY < 100 implies lossy compression, all other tested algorithms are lossless compression
Source: Original images were five uncompressed Orthophotos (2x2km) from LGLN. Original images were stored as uncompressed COGs (size: 539 ± 0 MB).

Elevation data

A 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

As elevation test data we used five terrain model tiles distributed in heterogeneous landscapes of Lower Saxony, Germany. They are also available as Cloud Optimized Geotiff from the State Office for Geoinformation and Surveying of Lower Saxony (LGLN) via s3 storage. Each of them covers an area of 1x1 km with a resolution of 1 m (1000x1000 px) and comes with a singe band (elecation) in data format Float32 (thirty two bit floating point), compressed with LZW.

# GDAL commands via OSGeo4W Shell


>gdalinfo /vsicurl/https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325735712/2016-04-04/dgm1_32_573_5712_1_ni_2016.tif
Driver: GTiff/GeoTIFF
Files: /vsicurl/https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325735712/2016-04-04/dgm1_32_573_5712_1_ni_2016.tif
Size is 1000, 1000
Coordinate System is:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            MEMBER["European Terrestrial Reference Frame 2020"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6┬░E and 12┬░E: Austria; Denmark - onshore and offshore; Germany - onshore and offshore; Italy - onshore and offshore; Norway including Svalbard - onshore and offshore; Spain - offshore."],
        BBOX[36.53,6,84.01,12.01]],
    USAGE[
        SCOPE["Pan-European conformal mapping at scales larger than 1:500,000."],
        AREA["Europe between 6┬░E and 12┬░E and approximately 36┬░30'N to 84┬░N."],
        BBOX[36.53,6,84.01,12.01]],
    ID["EPSG",25832]]
Data axis to CRS axis mapping: 1,2
Origin = (573000.000000000000000,5713000.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  TIFFTAG_SOFTWARE=LAStools is a product of rapidlasso GmbH, Germany
  TIFFTAG_ARTIST=created by LAStools (c) martin.isenburg@rapidlasso.com
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  573000.000, 5713000.000) ( 10d 3'11.36"E, 51d33'48.06"N)
Lower Left  (  573000.000, 5712000.000) ( 10d 3'10.61"E, 51d33'15.69"N)
Upper Right (  574000.000, 5713000.000) ( 10d 4' 3.28"E, 51d33'47.59"N)
Lower Right (  574000.000, 5712000.000) ( 10d 4' 2.52"E, 51d33'15.22"N)
Center      (  573500.000, 5712500.000) ( 10d 3'36.94"E, 51d33'31.64"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  NoData Value=-9999
  Overviews: 500x500

Compression

One-band Elevation

# pak::pak("USDAForestService/gdalraster")
library(gdalraster)
library(purrr)
library(fs)

# Sset gdal configurations (for reproducability?)
set_config_option("GDAL_NUM_THREADS", "16")
set_config_option("GDAL_CACHEMAX", "4000")
set_config_option("OVERVIEWS", "IGNORE_EXISTING")

urls <- c(
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325735712/2016-04-04/dgm1_32_573_5712_1_ni_2016.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326025743/2018-04-09/dgm1_32_602_5743_1_ni_2018.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326065787/2019-02-27/dgm1_32_606_5787_1_ni_2019.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326065850/2020-04-07/dgm1_32_606_5850_1_ni_2020.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325045930/2017-02-15/dgm1_32_504_5930_1_ni_2017.tif")

dsns <- paste0("/vsicurl/", urls)  # prefix for virtual file source



# Create a function to measure performance for each datasource and compression
raster_compression <- function(dsns, options) {
  results <- list()
  
  for (datasourcename in dsns) {
    for (option in options) {

      output_file <- file_temp(ext = ".tif")
      
      write_ds <- function(){
        createCopy(
          src_filename = datasourcename,
          dst_filename = output_file,
          format = "COG",
          options = option$setting
        )
      } 
      
      # Get write time
      writetime <- system.time(write_ds())["elapsed"]
      
      # Get file size in MB
      data_size_uncompressed <- vsi_stat(datasourcename, "size") / (1024^2) 
      data_size_compressed <- vsi_stat(output_file, "size") / (1024^2) 
      
      # Get read time
      img <- new(GDALRaster, output_file)
      readtime <- system.time(read_ds(img))["elapsed"]
      
      # Store results
      results[[length(results) + 1]] <- list(
        datasource = basename(datasourcename),
        settings = option$naming,
        uncompressed_mb = data_size_uncompressed,
        compressed_mb = data_size_compressed,
        writetime_sec = writetime,
        readtime_sec = readtime
      )
    }
  }
  
  # Convert list to data frame
  results_df <- do.call(rbind, lapply(results, data.frame))
  return(results_df)
}


# create list of options
options <- list(
  # lossless
  list(setting = c("COMPRESS=NONE")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=YES")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=NO")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=1")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=22")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=22")),
  list(setting = c("COMPRESS=LERC", "MAX_Z_ERROR=0")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=6")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=22"))
)

for (i in 1:length(options)) {
  options[[i]]$naming <- paste(options[[i]]$setting, collapse = ", ")
}


# apply function on list
benchmark_results <- raster_compression(dsns, options)

Comparison of Compression Performances
Size and processing time (mean ± standard deviation) depending on compression settings
Settings1	Absolute Size (MB)	Relative Size (%)	Runtime Write (s)	Runtime Read (s)
COMPRESS=NONE	5.00 ± 0.00	126.9 ± 9.6	1.61 ± 0.70	0.03 ± 0.01
COMPRESS=LERC, MAX_Z_ERROR=0	4.78 ± 0.00	121.3 ± 9.1	0.09 ± 0.01	0.04 ± 0.01
COMPRESS=LZW, PREDICTOR=NO	3.96 ± 0.31	100.0 ± 0.0	0.14 ± 0.01	0.04 ± 0.02
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=1	3.72 ± 0.18	94.3 ± 6.0	0.07 ± 0.01	0.03 ± 0.01
COMPRESS=LERC_ZSTD, LEVEL=1	3.72 ± 0.18	94.1 ± 6.1	0.10 ± 0.02	0.04 ± 0.01
COMPRESS=LZW, PREDICTOR=YES	3.37 ± 0.30	85.6 ± 11.8	0.12 ± 0.01	0.04 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=1	3.13 ± 0.35	78.8 ± 2.9	0.12 ± 0.02	0.04 ± 0.01
COMPRESS=LERC_DEFLATE, LEVEL=1	3.11 ± 0.35	78.5 ± 3.0	0.12 ± 0.03	0.04 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=6	3.02 ± 0.32	76.2 ± 2.4	0.15 ± 0.03	0.04 ± 0.01
COMPRESS=LERC_DEFLATE, LEVEL=6	3.00 ± 0.32	75.7 ± 2.4	0.17 ± 0.02	0.04 ± 0.01
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=9	2.89 ± 0.23	72.9 ± 1.8	0.17 ± 0.02	0.03 ± 0.01
COMPRESS=LERC_ZSTD, LEVEL=9	2.89 ± 0.23	72.9 ± 1.8	0.16 ± 0.03	0.05 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=9	2.86 ± 0.30	72.1 ± 2.6	0.22 ± 0.02	0.05 ± 0.01
COMPRESS=LERC_DEFLATE, LEVEL=9	2.84 ± 0.30	71.7 ± 2.6	0.24 ± 0.02	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=22	2.72 ± 0.26	68.7 ± 2.0	1.34 ± 0.06	0.04 ± 0.01
COMPRESS=LERC_ZSTD, LEVEL=22	2.71 ± 0.26	68.5 ± 2.1	0.50 ± 0.03	0.03 ± 0.02
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=1	2.60 ± 0.24	66.0 ± 8.9	0.11 ± 0.01	0.04 ± 0.01
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=1	2.56 ± 0.25	65.1 ± 9.1	0.09 ± 0.01	0.04 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=6	2.56 ± 0.24	65.0 ± 8.9	0.14 ± 0.03	0.04 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=9	2.53 ± 0.24	64.2 ± 8.9	0.33 ± 0.04	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=9	2.51 ± 0.23	63.7 ± 8.5	0.14 ± 0.02	0.04 ± 0.02
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=22	2.48 ± 0.22	63.0 ± 8.5	2.18 ± 0.10	0.04 ± 0.01
COMPRESS=LZMA, LEVEL=1	2.20 ± 0.09	55.9 ± 3.7	0.40 ± 0.01	0.06 ± 0.01
COMPRESS=LZMA, LEVEL=9	1.98 ± 0.13	50.0 ± 1.8	0.84 ± 0.02	0.06 ± 0.01
1 All tested algorithms are lossless compression, WEBP which performes best with Orthoimages is not possible with 1-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 storage and loading times. There are two different formats for this, varying in their encoding formula: Mapbox-encoded and Terrarium-encoded.

Bit-reduction

To reduce the file size even more we can reduce the bitdepth of the B-channel in RGB encoded elevation data. Orignial 8-bit data with Mapbox-encoding encodes the elevation in 10 cm steps. If we remove the last bit, data is stored in 128 steps, reducing the precision to 20 cm steps, and if we go further removing the last 2 bits results in 40 cm steps. Terrarium-encoding results in a precision of ~4 mm (1/256) with original 8-bit data, reducing bitdepth reduces precision by a factor of 2 here as well.

Mapbox-encoding

Mapbox encoding is usually the default in most applications which use RGB-encoded elevation. Mapbox, Maplibe, 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). RGB channels each store 8-bits, which allows for 16,777,216 unique values. 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.

# pak::pak("USDAForestService/gdalraster")
library(gdalraster)
library(purrr)
library(fs)

# Sset gdal configurations (for reproducability?)
set_config_option("GDAL_NUM_THREADS", "16")
set_config_option("GDAL_CACHEMAX", "4000")
set_config_option("OVERVIEWS", "IGNORE_EXISTING")

urls <- c(
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325735712/2016-04-04/dgm1_32_573_5712_1_ni_2016.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326025743/2018-04-09/dgm1_32_602_5743_1_ni_2018.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326065787/2019-02-27/dgm1_32_606_5787_1_ni_2019.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326065850/2020-04-07/dgm1_32_606_5850_1_ni_2020.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325045930/2017-02-15/dgm1_32_504_5930_1_ni_2017.tif")

dsns <- paste0("/vsicurl/", urls)  # prefix for virtual file source



# create list of options for RGB encoding
options <- list(
  # lossless
  list(setting = c("COMPRESS=NONE")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=YES")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=NO")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=1")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=22")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=22")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100")),
  list(setting = c("COMPRESS=LERC", "MAX_Z_ERROR=0")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=6")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=22")),
  # only available with certain gdal builds
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=1")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=5")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=9")),
  # lossy
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "DISCARD_LSB=0,0,1")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "DISCARD_LSB=0,0,2")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "DISCARD_LSB=0,0,3"))
)

for (i in 1:length(options)) {
  options[[i]]$naming <- paste(options[[i]]$setting, collapse = ", ")
}



#### benchmark mapbox encoding

# Create a function to measure performance for each datasource and compression
raster_compression <- function(dsns, options) {
  results <- list()
  
  for (datasourcename in dsns) {
    
    
    # convert DTM to RGB encoded (Mapbox)
    mapbox_encoded <- "to_terrainrgb <- function(dtm) {
        startingvalue <- 10000
        precision <- 0.1
        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))
        return(c(r,g,b))
    }
    
    to_terrainrgb(ELEV)"
    
    out_file <- tempfile(fileext = ".tif")
    
    
    dtm_mapbox <- calc(expr = mapbox_encoded,
                       rasterfiles = datasourcename,
                       var.names = "ELEV",
                       dstfile = out_file,
                       dtName = "Byte",
                       out_band = 1:3,
                       nodata_value = 0)
    
    
    
    for (option in options) {
      
      output_file <- file_temp(ext = ".tif")
      
      write_ds <- function(){
        createCopy(
          src_filename = dtm_mapbox,
          dst_filename = output_file,
          format = "COG",
          options = option$setting
        )
      } 
      
      # Get write time
      writetime <- system.time(write_ds())["elapsed"]
      
      # Get file size in MB
      data_size_uncompressed <- vsi_stat(datasourcename, "size") / (1024^2) 
      data_size_compressed <- vsi_stat(output_file, "size") / (1024^2) 
      
      # Get read time
      img <- new(GDALRaster, output_file)
      readtime <- system.time(read_ds(img))["elapsed"]
      
      # Store results
      results[[length(results) + 1]] <- list(
        datasource = basename(datasourcename),
        settings = option$naming,
        uncompressed_mb = data_size_uncompressed,
        compressed_mb = data_size_compressed,
        writetime_sec = writetime,
        readtime_sec = readtime
      )}
  }
  
  # Convert list to data frame
  results_df <- do.call(rbind, lapply(results, data.frame))
  return(results_df)
}

# apply function on list
benchmark_results_mapbox <- raster_compression(dsns, options)

Comparison of Compression Performances
Size and processing time (mean ± standard deviation) depending on compression settings
Settings1	Absolute Size (MB)	Relative Size (%)	Runtime Write (s)	Runtime Read (s)
COMPRESS=NONE	3.75 ± 0.00	95.2 ± 7.2	0.15 ± 0.02	0.06 ± 0.01
COMPRESS=LZW, PREDICTOR=NO	1.13 ± 0.54	27.8 ± 11.0	0.15 ± 0.02	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=1	0.74 ± 0.26	18.5 ± 4.9	0.15 ± 0.03	0.06 ± 0.02
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=1	0.69 ± 0.26	17.2 ± 5.1	0.15 ± 0.01	0.05 ± 0.02
COMPRESS=LERC, MAX_Z_ERROR=0	0.66 ± 0.16	16.5 ± 2.8	0.18 ± 0.04	0.05 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=6	0.60 ± 0.26	14.9 ± 5.3	0.21 ± 0.04	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=9	0.58 ± 0.22	14.5 ± 4.3	0.21 ± 0.04	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=1	0.57 ± 0.13	14.3 ± 2.1	0.16 ± 0.01	0.06 ± 0.00
COMPRESS=LZMA, LEVEL=1	0.57 ± 0.22	14.2 ± 4.4	0.25 ± 0.02	0.08 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=9	0.54 ± 0.26	13.4 ± 5.2	0.49 ± 0.08	0.05 ± 0.02
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=1	0.53 ± 0.12	13.2 ± 2.0	0.14 ± 0.02	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=22	0.49 ± 0.22	12.1 ± 4.4	1.58 ± 0.15	0.05 ± 0.02
COMPRESS=LZMA, LEVEL=9	0.48 ± 0.22	11.9 ± 4.3	0.78 ± 0.08	0.06 ± 0.01
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=9	0.47 ± 0.11	11.9 ± 1.9	0.26 ± 0.06	0.05 ± 0.02
COMPRESS=LZW, PREDICTOR=YES	0.47 ± 0.15	11.8 ± 2.8	0.18 ± 0.04	0.05 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=6	0.46 ± 0.12	11.5 ± 2.1	0.21 ± 0.01	0.05 ± 0.02
COMPRESS=LERC_ZSTD, LEVEL=1	0.39 ± 0.11	9.8 ± 1.9	0.20 ± 0.03	0.07 ± 0.02
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=9	0.39 ± 0.11	9.8 ± 2.0	0.64 ± 0.06	0.05 ± 0.01
COMPRESS=LERC_DEFLATE, LEVEL=1	0.38 ± 0.11	9.5 ± 1.9	0.19 ± 0.02	0.06 ± 0.03
COMPRESS=LERC_ZSTD, LEVEL=9	0.38 ± 0.11	9.4 ± 2.0	0.18 ± 0.02	0.06 ± 0.02
COMPRESS=LERC_DEFLATE, LEVEL=6	0.37 ± 0.11	9.3 ± 2.0	0.19 ± 0.02	0.06 ± 0.01
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=22	0.37 ± 0.11	9.2 ± 2.0	1.79 ± 0.10	0.05 ± 0.01
COMPRESS=LERC_DEFLATE, LEVEL=9	0.37 ± 0.11	9.1 ± 2.0	0.23 ± 0.03	0.06 ± 0.02
COMPRESS=LERC_ZSTD, LEVEL=22	0.36 ± 0.11	9.1 ± 2.0	0.27 ± 0.03	0.06 ± 0.02
COMPRESS=WEBP, QUALITY=100	0.23 ± 0.06	5.9 ± 1.0	0.44 ± 0.12	0.06 ± 0.02
COMPRESS=WEBP, QUALITY=100, DISCARD_LSB=0,0,1	0.23 ± 0.06	5.9 ± 1.0	0.44 ± 0.10	0.06 ± 0.02
COMPRESS=WEBP, QUALITY=100, DISCARD_LSB=0,0,2	0.23 ± 0.06	5.9 ± 1.0	0.44 ± 0.10	0.06 ± 0.02
COMPRESS=WEBP, QUALITY=100, DISCARD_LSB=0,0,3	0.23 ± 0.06	5.9 ± 1.0	0.47 ± 0.12	0.04 ± 0.01
1 All tested algorithms are lossless compression, except for DISCARD_LSB where the last 1/2/3 bits of the blue band are discarted.
Source: Original images were converted to RGB-encoded images first (Mapbox-encoding)

Terrarium-encoding

The Terrarium-encoding uses a different formula to store elevation in three bands, it can be decoded as:

height = (R * 256 + G + B / 256) - 32768

This allows for a precision of 0.00390625 m (1 / 256).

# pak::pak("USDAForestService/gdalraster")
library(gdalraster)
library(purrr)
library(fs)

# Sset gdal configurations (for reproducability?)
set_config_option("GDAL_NUM_THREADS", "16")
set_config_option("GDAL_CACHEMAX", "4000")
set_config_option("OVERVIEWS", "IGNORE_EXISTING")

urls <- c(
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325735712/2016-04-04/dgm1_32_573_5712_1_ni_2016.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326025743/2018-04-09/dgm1_32_602_5743_1_ni_2018.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326065787/2019-02-27/dgm1_32_606_5787_1_ni_2019.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/326065850/2020-04-07/dgm1_32_606_5850_1_ni_2020.tif",
  "https://dgm.s3.eu-de.cloud-object-storage.appdomain.cloud/325045930/2017-02-15/dgm1_32_504_5930_1_ni_2017.tif")

dsns <- paste0("/vsicurl/", urls)  # prefix for virtual file source



# create list of options for RGB encoding
options <- list(
  # lossless
  list(setting = c("COMPRESS=NONE")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=YES")),
  list(setting = c("COMPRESS=LZW", "PREDICTOR=NO")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=6")),
  list(setting = c("COMPRESS=DEFLATE", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=1")),
  list(setting = c("COMPRESS=LZMA", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=YES", "LEVEL=22")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=1")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=9")),
  list(setting = c("COMPRESS=ZSTD", "PREDICTOR=NO", "LEVEL=22")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100")),
  list(setting = c("COMPRESS=LERC", "MAX_Z_ERROR=0")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=6")),
  list(setting = c("COMPRESS=LERC_DEFLATE", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=1")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=9")),
  list(setting = c("COMPRESS=LERC_ZSTD", "LEVEL=22")),
  # only available with certain gdal builds
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=1")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=5")),
  # list(setting = c("COMPRESS=JXL", "JXL_LOSSLESS=YES", "JXL_EFFORT=9")),
  # lossy
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "DISCARD_LSB=0,0,1")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "DISCARD_LSB=0,0,2")),
  list(setting = c("COMPRESS=WEBP", "QUALITY=100", "DISCARD_LSB=0,0,3"))
)

for (i in 1:length(options)) {
  options[[i]]$naming <- paste(options[[i]]$setting, collapse = ", ")
}



#### benchmark terrarium encoding

# Create a function to measure performance for each datasource and compression
raster_compression <- function(dsns, options) {
  results <- list()
  
  for (datasourcename in dsns) {
    
    
    # convert DTM to RGB encoded (Terrarium)
    terrarium_encoded <- "to_terrainrgb <- function(dtm) {
        v <- dtm + 32768
        r <- floor(v/256)
        g <- floor(v %% 256)
        b <- floor((v - floor(v)) * 256)
        return(c(r,g,b))
        }
        to_terrainrgb(ELEV)"
    
    out_file <- tempfile(fileext = ".tif")
    
    dtm_terrarium <- calc(expr = terrarium_encoded,
                          rasterfiles = datasourcename,
                          var.names = "ELEV",
                          dstfile = out_file,
                          dtName = "Byte",
                          out_band = 1:3,
                          nodata_value = 0)
    
    
    
    for (option in options) {
      
      output_file <- file_temp(ext = ".tif")
      
      write_ds <- function(){
        createCopy(
          src_filename = dtm_terrarium,
          dst_filename = output_file,
          format = "COG",
          options = option$setting
        )
      } 
      
      # Get write time
      writetime <- system.time(write_ds())["elapsed"]
      
      # Get file size in MB
      data_size_uncompressed <- vsi_stat(datasourcename, "size") / (1024^2) 
      data_size_compressed <- vsi_stat(output_file, "size") / (1024^2) 
      
      # Get read time
      img <- new(GDALRaster, output_file)
      readtime <- system.time(read_ds(img))["elapsed"]
      
      # Store results
      results[[length(results) + 1]] <- list(
        datasource = basename(datasourcename),
        settings = option$naming,
        uncompressed_mb = data_size_uncompressed,
        compressed_mb = data_size_compressed,
        writetime_sec = writetime,
        readtime_sec = readtime
      )
    }
  }
  
  # Convert list to data frame
  results_df <- do.call(rbind, lapply(results, data.frame))
  return(results_df)
}

# apply function on list
benchmark_results_terrarium <- raster_compression(dsns, options)

Comparison of Compression Performances
Size and processing time (mean ± standard deviation) depending on compression settings
Settings1	Absolute Size (MB)	Relative Size (%)	Runtime Write (s)	Runtime Read (s)
COMPRESS=NONE	3.75 ± 0.00	95.2 ± 7.2	0.15 ± 0.02	0.04 ± 0.02
COMPRESS=LZW, PREDICTOR=NO	2.66 ± 0.40	66.9 ± 4.8	0.18 ± 0.02	0.05 ± 0.01
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=1	2.21 ± 0.37	55.4 ± 5.0	0.17 ± 0.03	0.06 ± 0.01
COMPRESS=LZMA, LEVEL=1	1.97 ± 0.43	49.3 ± 6.9	0.38 ± 0.04	0.08 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=1	1.90 ± 0.29	47.8 ± 3.7	0.16 ± 0.02	0.06 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=6	1.86 ± 0.31	46.7 ± 4.2	0.18 ± 0.02	0.06 ± 0.02
COMPRESS=DEFLATE, PREDICTOR=NO, LEVEL=9	1.81 ± 0.32	45.5 ± 4.5	0.29 ± 0.04	0.05 ± 0.01
COMPRESS=LZMA, LEVEL=9	1.80 ± 0.37	45.2 ± 5.6	0.70 ± 0.03	0.07 ± 0.03
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=9	1.71 ± 0.24	43.0 ± 2.7	0.21 ± 0.02	0.06 ± 0.02
COMPRESS=ZSTD, PREDICTOR=NO, LEVEL=22	1.62 ± 0.26	40.8 ± 3.4	1.33 ± 0.03	0.05 ± 0.02
COMPRESS=LZW, PREDICTOR=YES	1.61 ± 0.25	40.5 ± 3.4	0.20 ± 0.03	0.06 ± 0.02
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=1	1.43 ± 0.20	36.1 ± 2.4	0.17 ± 0.01	0.04 ± 0.01
COMPRESS=LERC, MAX_Z_ERROR=0	1.41 ± 0.17	35.6 ± 1.8	0.18 ± 0.03	0.06 ± 0.02
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=1	1.41 ± 0.19	35.4 ± 2.3	0.19 ± 0.03	0.07 ± 0.03
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=6	1.40 ± 0.19	35.1 ± 2.3	0.22 ± 0.02	0.05 ± 0.02
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=9	1.36 ± 0.17	34.2 ± 2.0	0.23 ± 0.01	0.05 ± 0.03
COMPRESS=LERC_ZSTD, LEVEL=1	1.33 ± 0.16	33.6 ± 2.2	0.18 ± 0.01	0.06 ± 0.01
COMPRESS=DEFLATE, PREDICTOR=YES, LEVEL=9	1.33 ± 0.19	33.4 ± 2.4	0.37 ± 0.03	0.06 ± 0.02
COMPRESS=LERC_DEFLATE, LEVEL=1	1.28 ± 0.15	32.4 ± 2.0	0.20 ± 0.01	0.05 ± 0.02
COMPRESS=LERC_DEFLATE, LEVEL=6	1.27 ± 0.15	32.1 ± 1.9	0.19 ± 0.01	0.06 ± 0.01
COMPRESS=LERC_ZSTD, LEVEL=9	1.27 ± 0.15	32.1 ± 1.9	0.20 ± 0.02	0.06 ± 0.01
COMPRESS=LERC_DEFLATE, LEVEL=9	1.27 ± 0.14	32.0 ± 1.9	0.24 ± 0.02	0.08 ± 0.03
COMPRESS=LERC_ZSTD, LEVEL=22	1.25 ± 0.14	31.6 ± 1.8	0.28 ± 0.03	0.07 ± 0.01
COMPRESS=ZSTD, PREDICTOR=YES, LEVEL=22	1.23 ± 0.17	30.9 ± 2.2	1.46 ± 0.03	0.05 ± 0.01
COMPRESS=WEBP, QUALITY=100	0.85 ± 0.11	21.3 ± 1.4	0.97 ± 0.14	0.06 ± 0.01
COMPRESS=WEBP, QUALITY=100, DISCARD_LSB=0,0,1	0.85 ± 0.11	21.3 ± 1.4	1.00 ± 0.12	0.05 ± 0.01
COMPRESS=WEBP, QUALITY=100, DISCARD_LSB=0,0,2	0.85 ± 0.11	21.3 ± 1.4	1.00 ± 0.16	0.05 ± 0.02
COMPRESS=WEBP, QUALITY=100, DISCARD_LSB=0,0,3	0.85 ± 0.11	21.3 ± 1.4	0.97 ± 0.17	0.05 ± 0.01
1 All tested algorithms are lossless compression, except for DISCARD_LSB where the last 1/2/3 bits of the blue band are discarted.
Source: Original images were converted to RGB-encoded images first (Terrarium-encoding)

Hillshade

If data is only used for visualization it might be an alternative to use a hillshade, this can also dramatically reduce file size, as it is only nessecary to store e.g. 256 greyscale values. However, this limits the data usability dramatically to certain use-cases which is probably not wanted.

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. 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. To get a decent image quality at high compression rates it might be a viable solution to vary compression ratio (or algorithm) by zoom level. We can have lower compression with higher quality at high zoom levels (or original data) and high compression rates at lower zoom levels (overviews) as this is artificially resampled data anyway.

For elevation data good compression ratios are achieved when the data is encoded in RGB-values and stored as type byte. 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

• GeoTiff Compression for Dummies
• Der Weg zum optimalen GeoTIFF
• Guide to GeoTIFF compression and optimization with GDAL
• TIFF compression options by Keith Jenkins
• Optimization of RGB DEM tiles for dynamic hill shading with Mapbox GL or MapLibre GL
• Elevation compression and bulk distribution
• GeoTiff compression benchmarking
• Best Practices for Cloud-Optimized GeoTIFFs (COGs) for the use in Swiss Geospatial Infastructures