Create RGBA COG with GDAL > 2.3

create_rgb_bigtiff_cog.sh

#!/bin/bash

set -Eeuo pipefail

BLOCKSIZE=256
OVERVIEW_BLOCKSIZE=256
MIN_OVERVIEW_SIZE=256
KEEP_TEMP=0
MOSAIC_VRT="mosaic.vrt"
MOSAIC_RGB_VRT="mosaicrgb.vrt"
TIFF_FILES=file.list

usage() {
  cat >&2 <<EOF
USAGE: $0 [-k] [-t tile_size] <SOURCE_DIR> <OUTPUT_TIFF>

     -h            display this help and exit
     -t TILE_SIZE  internal tiles size for the source data and overviews
     -m SIZE       Maximum width or height of the smallest overview level. Defaults to 256.
     -k            keep temp files from the processing.
EOF
  exit 1
}

function cleanup {
  if [ "$KEEP_TEMP" == 0 ]; then
    rm -f $MOSAIC_VRT
    rm -f $MOSAIC_RGB_VRT
    rm -f $MOSAIC_RGB_VRT.*
    rm -f $TIFF_FILES
  fi
}
trap cleanup EXIT

function version {
  echo "$@" | awk -F. '{ printf("%d%03d%03d%03d\n", $1,$2,$3,$4); }';
}

function to_abs_path {
    local target="$1"
    if [ "$target" == "." ]; then
        echo "$(pwd)"
    elif [ "$target" == ".." ]; then
        echo "$(dirname "$(pwd)")"
    else
        echo "$(cd "$(dirname "$1")"; pwd)/$(basename "$1")"
    fi
}

# Note: for JPEG compression, the above method produce cloud optimized files only if using GDAL 2.2 
# (or a dev version >= r36879). For older versions, the IFD of the overviews will be written towards
# the end of the file. A recent version of GDAL (2.2 or dev version >= r37257) built against internal
# libtiff (or libtiff >= 4.0.8) will also help reducing the amount of bytes read for JPEG compressed
# files with YCbCr subsampling. Lastly GDAL 2.3 is required for ZSTD compression, which is much faster
# than deflate for the steps before the final gdal_translate that compresses to JPEG
# See https://trac.osgeo.org/gdal/wiki/CloudOptimizedGeoTIFF#HowtogenerateitwithGDAL
GDAL_VERSION=$(gdal_translate --version | cut -c 6- | rev | cut -c 22- | rev)
if [ $(version $GDAL_VERSION) -lt $(version "2.3.0") ]; then
    echo "ERROR: Need to have at least GDAL version 2.3.0 installed, only have $GDAL_VERSION" >&2
    exit 1
fi

if ! (gdal_translate --format GTiff || true) | grep --quiet "<Value>ZSTD"; then
  echo "GDAL ZSTD compressions not found" >&2
  exit 1
fi
if ! (gdal_translate --format GTiff || true) | grep --quiet "LIBTIFF=INTERNAL"; then
  echo "GDAL Internal not LIBTIFF used" >&2
  exit 1
fi

while getopts hkt:m: opt; do
  case $opt in
    h)
      usage
      ;;
    k) KEEP_TEMP=1
      ;;
    m) MIN_OVERVIEW_SIZE=$OPTARG
      ;;
    t)
      BLOCKSIZE=$OPTARG
      OVERVIEW_BLOCKSIZE=$OPTARG
      ;;
    *)
      usage
      ;;
  esac
done
shift $((OPTIND -1))

echo "Internal Block Size: $BLOCKSIZE"
echo "Internal Overview Block Size: $OVERVIEW_BLOCKSIZE"
echo "Min overview height size: $MIN_OVERVIEW_SIZE"

if [ "$#" -ne 2 ]; then
    usage
fi

SOURCE_DIR=$(to_abs_path $1)
OUTPUT_TIFF=$(to_abs_path $2)

if [ ! -d  "$SOURCE_DIR" ]; then
    echo "ERROR: The source directory $SOURCE_DIR does not exist or is not a directory" >&2
    usage
fi
if [ -z "$OUTPUT_TIFF" ]; then
    echo "ERROR: define the OUTPUT_TIFF parameter" >&2
    usage
fi

# Create a mosaic of the tiles, with transparency mask
OUTPATH_PATH=$(dirname "${OUTPUT_TIFF}")
cd "$OUTPATH_PATH"

find "$SOURCE_DIR" -maxdepth 1 -name "*.tif" > $TIFF_FILES
gdalbuildvrt -addalpha -input_file_list $TIFF_FILES $MOSAIC_VRT
gdal_translate -of VRT -b 1 -b 2 -b 3 $MOSAIC_VRT $MOSAIC_RGB_VRT

# Create external overviews by generating each overview level in its 
# own file by cascading calls to gdaladdo. Only in the final gdal_translate
# COG generation stage we will  use --config GDAL_TIFF_INTERNAL_MASK YES. 
# This is a workaround to achieve better performance of very large rasters. 
# See https://trac.osgeo.org/gdal/ticket/5067#comment:2
#
# Note: Also it seems that external overview building doesn't automatically 
# generate overviews of the associate mask. So need to generate external mask
# at those intermediate stages  (that is do no set GDAL_TIFF_INTERNAL_MASK),
# and run gdaladdo on the .msk  file, generating a .msk.ovr, and then run 
# it on the .msk.ovr, etc. See this on how to generate the msk file from 
# a rgba vrt file. See: 
# https://lists.osgeo.org/pipermail/gdal-dev/2014-February/038107.html
gdal_translate \
  -of "GTiff" \
  -b 4 \
  -mo "INTERNAL_MASK_FLAGS_1=2" \
  -mo "INTERNAL_MASK_FLAGS_2=2" \
  -mo "INTERNAL_MASK_FLAGS_3=2" \
  -co COMPRESS=ZSTD \
  -co INTERLEAVE=BAND \
  -co NBITS=1 \
  --config GDAL_NUM_THREADS ALL_CPUS \
  $MOSAIC_VRT $MOSAIC_RGB_VRT.msk

# TODO autoselect levels like gdaladdo 2.3+ does
OVERVIEW_FILE=$MOSAIC_RGB_VRT
OVERVIEW_MSK_FILE=$MOSAIC_RGB_VRT.msk
for LEVEL in 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192
do
  echo "Building overview level: $LEVEL"
  gdaladdo \
    --config COMPRESS_OVERVIEW ZSTD \
    --config GDAL_TIFF_INTERNAL_MASK NO \
    --config GDAL_NUM_THREADS ALL_CPUS \
    -ro \
    -r average \
    $OVERVIEW_FILE 2
  # For the first iteration the first gdaladdo builds a mask, so skip the specific
  # mask building process to save a little time. For the subsequent iterations
  # the link to the mask will be lost due to the naming convention and the we
  # need to run the gdaladdo mask process
  if [ "$LEVEL" -ne 2 ]
  then
    gdaladdo \
      --config COMPRESS_OVERVIEW ZSTD \
      --config GDAL_TIFF_INTERNAL_MASK NO \
      --config GDAL_NUM_THREADS ALL_CPUS \
      -ro \
      -r average \
      $OVERVIEW_MSK_FILE 2
  fi
  OVERVIEW_FILE=${OVERVIEW_FILE}.ovr
  OVERVIEW_MSK_FILE=${OVERVIEW_MSK_FILE}.ovr
done

# Create the Cloud Optimised Geotiff
# Note Increasing the block size can reduce the size of the IFD. But larger 
# blocks can cause more bytes to be pulled for random access if the
# compression rate is not high.
gdal_translate \
  -of GTiff \
  -co BIGTIFF=YES \
  -co TILED=YES \
  -co BLOCKXSIZE=$BLOCKSIZE \
  -co BLOCKYSIZE=$BLOCKSIZE \
  -co COMPRESS=JPEG \
  -co JPEG_QUALITY=85 \
  -co PHOTOMETRIC=YCBCR \
  -co COPY_SRC_OVERVIEWS=YES \
  -co ALPHA=YES \
  --config GDAL_TIFF_OVR_BLOCKSIZE $OVERVIEW_BLOCKSIZE \
  --config GDAL_TIFF_INTERNAL_MASK YES \
  $MOSAIC_RGB_VRT "$OUTPUT_TIFF"

# Finally validate the COG GeoTiff
validate_cloud_optimized_geotiff.py "$OUTPUT_TIFF"

# vim: set ts=4 sts=4 sw=4 et: