Skip to content

Instantly share code, notes, and snippets.

@palmerj
Last active April 6, 2022 20:42
Show Gist options
  • Save palmerj/2e4a46fbcf0c97212e6e77fced22e885 to your computer and use it in GitHub Desktop.
Save palmerj/2e4a46fbcf0c97212e6e77fced22e885 to your computer and use it in GitHub Desktop.
Create RGBA COG with GDAL > 2.3
#!/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:
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment