Search code examples
pythongdal

Merge overlapping rasters using GDAL

I have around 211 rasters, one for each area of the world. The gdalinfo for one of them is below. They are all the same except for the area of concern. Values are always between 1-6.

enter image description here

I would like to merge them into one giant raster for the world. I have managed to successfully do this by first building a VRT, and then writing the VRT to a single file:

gdalbuildvrt.exe -b 1 -q -input_file_list my_files.txt global_file.vrt

gdal_translate.exe -q -co PREDICTOR=2 -co COMPRESS=LZW -of GTiff -co BIGTIFF=YES -co TILED=YES -co NUM_THREADS=ALL_CPUS global_file.vrt global_file.tif

The result was around 15GB in size.

However my issue is that there are often areas where each country overlap. In those cases I need to take the max raster/pixel value. But gdal_translate does not do that, it just takes the last value that was written.

I read about the PixelFunction, (https://gis.stackexchange.com/questions/350233/how-can-i-modify-a-vrtrasterband-sub-class-etc-from-python) and tried to implement it, but I kept getting memory issues.

Does anyone have any ideas on a memory-safe way to combine a large list of rasters and to take the maximum value where they overlap? If the best way is the PixelFunction then let me know and I will provide more details about the errors.

Thanks, James

Solution

  • Try reducing your output block size, use BLOCKXSIZE=128 and BLOCKYSIZE=128:

    gdal_translate.exe -q -co PREDICTOR=2 -co COMPRESS=LZW -of GTiff -co BIGTIFF=YES -co BLOCKXSIZE=128 -co BLOCKYSIZE=128 -co TILED=YES -co NUM_THREADS=ALL_CPUS global_file.vrt global_file.tif

    This will require significantly less memory.

    If this doesn't work, your next option would be to convert all of your input files to 128x128 tiled geotiffs. Alas, in GDAL <3.7 you can't change the block size of the VRT driver, it is something that is coming in 3.7.