Fastest Possible Use of GDAL to Merge, Reproject & Convert

Question

I have a ~10GB, Australia-wide, single band dataset consisting of 3,351 GeoTIFFs across 8 projections (Map Grid of Australia [MGA], zones 49 to 56). The data has a 2m resolution but only covers urban areas; the images are mostly NODATA.

I'm converting them into a single GDA94 lat/long Cloud Optimised GeoTIFF (COG) to serve the data from AWS S3 and to simplify visualisation & analysis of the data at a national scale.

Whilst merging each MGA zone's data & reprojecting them to GDA94 was fast enough (around 1-2 hours), merging the resulting 8 images and the conversion to a single COG took 8-12 hours.

In my Python script I've tested many combinations of gdal.BuildVRT, gdal.Warp, gdal.Translate & gdal_merge with a subset of the images but can't get it to process any faster (noting it's possible to make it slower with a poor combination).

The fastest process I've found is:

1. Merge & reproject each MGA zone's set of images to GDA94 using gdal.Warp in a single step
2. Create a VRT of the 8 GDA94 images
3. Use gdal.Translate to merge the 8 images in the VRT and output as a single COG

Is there a faster way of doing this?:

# mosaic and transform to GDA94 lat/long for each MGA zone (aka UTM South zones on GDA94 datum)
for zone in mga_zones:
    files_to_mosaic = glob.glob(os.path.join(input_dict["input_path"], f"*{input_dict['glob_pattern']}Z{zone}*.tif"))
    num_images = len(files_to_mosaic)

    # store the interim images in memory
    interim_file = f"/vsimem/temp_Z{zone}_{input_dict['name']}.tif"

    # merge and convert to GDA94
    gdal.Warp(interim_file, files_to_mosaic, options="-multi -wm 80% -t_srs EPSG:4283 -co BIGTIFF=YES -co COMPRESS=DEFLATE -co NUM_THREADS=ALL_CPUS -overwrite")
    warped_files_to_mosaic.append(interim_file)

# mosaic all merged files and output as a single Cloud Optimised GeoTIFF (COG) for all of AU
vrt_file = f"temp_au_{input_dict['name']}.vrt"
vrt = gdal.BuildVRT(vrt_file, warped_files_to_mosaic)

gdal.Translate(input_dict["output_file"], vrt, format="COG", options="-co BIGTIFF=YES -co COMPRESS=DEFLATE -co NUM_THREADS=ALL_CPUS")
vrt = None
os.remove(vrt_file)

It's step 3 that takes 8-12 hours on an AWS EC2 instance with 64 cores and 512GB RAM (with Amazon Linux, i.e. Fedora). It's running on Python 3.9 with GDAL 3.3.2 in a Conda environment.

Observations:

• Using gdal.BuildVRT followed by gdal.Warp for step 1 caused a 2-3x slowdown (?)
• Using gdal_merge.py & gdal.Translate for steps 2 & 3 is slower than gdal.BuildVRT with gdal.Translate
• Setting gdal.SetCacheMax to 480GB didn't appear to have any effect on performance (noting I ran a subset of images for this testing)
• Using processes with all CPUs gives a modest speedup. Python processes mostly run at 100-200% CPU; indicating limited parallelisation
• RAM is potentially underutilised; presumedly because this is a CPU heavy process. The Python process usually sits around 10-20% memory usage
• I also tested Rasterio to see the difference between C & "pure" Python and it was significantly slower with Rasterio; somewhat expected but was running out of ideas...
• Running the equivalent process in a Bash script made no noticable difference

Answer 1

On further testing: discovered that gdal.Warp supports reprojecting, mosaicing & converting images with different projections in one go.

Simplifying the code to one line sped up the process by over 10x.

This takes 1 hour 15 mins to output a 22GB image from the 3,351 input files across 8 projections:

gdal.Warp(output_file, input_files, format="COG", 
      options="-overwrite -multi -wm 80% -t_srs EPSG:4283 -co TILED=YES -co BIGTIFF=YES -co COMPRESS=DEFLATE -co NUM_THREADS=ALL_CPUS")