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")

# 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

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.


  • 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
  • I've seen the same before, the handling of VRT is generic and having a more specific understanding of the data can allow you to cut some corners. I have achieved huge speed benefits by taking control of the operation in C++ (up to the disc/network read-write speed), using python or C#/VB.net you could make the process faster with GDAL bindings and taking charge of direct raster I/O. Commented Oct 27, 2021 at 7:10
  • 1
    Create the 8 GDA94 images as tiled -co TILED=YES and help GDAL to use memory efficiently. For writing a 256x256 sized tile into COG GDAL must read at least 256 whole rows from a striped TIFF and they are perhaps tens of thousands of pixels wide. GDAL usually does not need or gain from huge cachemax. If GDAL can read a small chunks of pixel data, process, and write to disk, the memory requirement is not big at all even the images may be.
    – user30184
    Commented Oct 27, 2021 at 8:04

1 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")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.