How to merge multiple GeoTIFFs into one single GeoTIFF and change the resolution

Question

I'm working on my Bachelor Thesis right now and I have a about 300gb of digital elevation database. This database has a 1m resolution and every single 1x1km square is stored into a single GeoTIFF. Now I need all of this data merged into a single grid.

I tried to work with a batch, but the only solution i could find is to first transform every GeoTIFF into a sgrd and then merge them all together and resample it. Is there a more efficient way? It takes very very long for all this data. My batch looked like this:

@echo off
setlocal enabledelayedexpansion

:: Bestimme das Verzeichnis, in dem das Skript liegt
set "SCRIPT_DIR=%~dp0"

:: Setze das SAGA Moduls PATH
set "SAGA_CMD_PATH=C:\Program Files\saga-6.4.0_x64\saga-6.4.0_x64\saga_cmd.exe"

:: Schleife durch alle TIFF-Dateien im Verzeichnis und konvertiere sie in SAGA-Grid-Dateien
for %%G in ("%SCRIPT_DIR%*.tif") do (
    echo Konvertiere %%G in %%~nG.sgrd
    "%SAGA_CMD_PATH%" io_gdal 0 -GRIDS "%%~nG.sgrd" -FILES "%%G"
    if %ERRORLEVEL% neq 0 (
        echo Fehler bei der Konvertierung von %%G
        exit /b 1
    )
)

:: Erzeuge eine Liste der SAGA-Grid-Dateien im Verzeichnis
set "GRIDS="
for %%H in ("%SCRIPT_DIR%*.sgrd") do (
    set "GRIDS=!GRIDS!;%%H"
)

:: Entferne das erste Semikolon
set "GRIDS=%GRIDS:~1%"

:: Überprüfen, ob SAGA-Grid-Dateien gefunden wurden
if "%GRIDS%"=="" (
    echo Keine SAGA-Grid-Dateien im Verzeichnis %SCRIPT_DIR% gefunden.
    exit /b 1
)

:: Setze die Ausgabedatei
set "OUTPUT_FILE=%SCRIPT_DIR%merged_output.sgrd"

:: Führe das Mosaicking aus
echo Führe das Mosaicking aus
"%SAGA_CMD_PATH%" grid_tools 3 -GRIDS "%GRIDS%" -MERGED "%OUTPUT_FILE%"

:: Überprüfen, ob das Mosaicking erfolgreich war
if %ERRORLEVEL% equ 0 (
    echo Mosaicking erfolgreich abgeschlossen. Ausgabedatei: %OUTPUT_FILE%
) else (
    echo Fehler beim Mosaicking.
    exit /b 1
)

:: Konvertiere die endgültige SAGA-Grid-Datei zurück in GeoTIFF
echo Konvertiere das Mosaik in GeoTIFF
"%SAGA_CMD_PATH%" io_gdal 2 -GRIDS "%OUTPUT_FILE%" -FILES "%SCRIPT_DIR%merged_output.tif" -OUT_CELL_SIZE 30.0

:: Aufräumen: Lösche die temporären SAGA-Grid-Dateien
del "%SCRIPT_DIR%*.sgrd"
del "%SCRIPT_DIR%*.mgrd"
del "%SCRIPT_DIR%*.sdat"
del "%SCRIPT_DIR%*.prj"

echo Fertig.

Does anyone has an Idea to do it more efficient? Maybe with Python?

Answer 1

I highly recommend becoming comfortable with the GDAL command line applications if you are going to be working with raster files at all regularly.

For your question you can use gdalbuildvrt to construct a virtual single tiff from your original files:

gdalbuildvrt merged.vrt *.tif 

Then to convert that raster and change the resolution you would want gdal_translate with something like:

gdal_translate -co COMPRESSED=LZW -co TILED=YES -tr 100 100 merged.vrt new_tiff.tif

Finally,I would recommend that you read Paul Ramsey's Compression for Dummies blog post to learn more about how to make imagery much smaller (and faster) to read and store.

Answer 2

Welcome to GIS community! My solution is based on rasterio which is a Python package for raster processing.

The first function called find_bill_files_path saves the paths of all .tiff files to a list. Its parameter is path_dir where the tiff images are saved locally.

The second function called merge_tiff merges all raster images into one raster. The two parameters are: (1) the list with the paths from all images, (2) the path of the merged image that will be created.

The third function called resampling_image resamples an input raster image to a new cell size and saves it to an output raster file. You may select the new resolution.

I have used .vrt (virtual rastr) instead of .tiff for the merged image. VRTs allow for efficient management of large datasets by describing the metadata of the actual raster files, reducing the need to store and process the entire dataset.

Be aware that resampling_image function is not necessary since I could have integrated resolution change into merge_tiff function.

import rasterio
import os
from rasterio.merge import merge
from rasterio.warp import calculate_default_transform, reproject, Resampling

def find_bill_files_path(path_dir):
    tiff_files = []
    
    for file in os.listdir(path_dir):
        if file.endswith('.tif') :
            # Construct the full file path
            input_raster_path = os.path.join(path_dir, file)
            bill_files.append(input_raster_path)
    
    return tiff_files


def merge_tiff(tiff_files, output_file):
    src_files_to_mosaic = []
    for fp in tiff_files:
        src = rasterio.open(fp)
        src_files_to_mosaic.append(src)
    
    # Merge function from Rasterio
    mosaic, out_trans = merge(src_files_to_mosaic)
    
    # Metadata of the output file
    out_meta = src_files_to_mosaic[0].meta.copy()
    out_meta.update({
        "driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": out_trans
    })
    
    # Save the merged TIFF
    with rasterio.open(output_file, 'w', **out_meta) as dest:
        dest.write(mosaic)
    
    # Close all the source files
    for src in src_files_to_mosaic:
        src.close()

def resampling_image(input_raster_path, output_raster_path, cell_size=6):
    with rasterio.open(input_raster_path) as src:
        # Define the desired new cell size
        new_cell_size = cell_size  # Desired cell size in the same units as the original raster
    
        # Compute the new transform and dimensions
        transform, width, height = calculate_default_transform(
            src.crs, src.crs, src.width, src.height, *src.bounds, resolution=new_cell_size)
    
        # Define the new metadata
        new_meta = src.meta.copy()
        new_meta.update({
            'transform': transform,
            'width': width,
            'height': height,
            'crs': src.crs
        })
    
        # Resample the raster
        with rasterio.open(output_raster_path, 'w', **new_meta) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=src.crs,
                    resampling=Resampling.bilinear)  # Change to Resampling.nearest for nearest neighbor

# Implementation
tiff_files_list = find_tiff_files_path(r"F:\DEM\all_images")
output_file_merged = r"F:\DEM\outputs\merged_image.vrt"
output_file_merged_tiff_resampled = r"F:\DEM\outputs\merged_image_resampled.tiff"
new_resolution = 10 

merge_tiff(tiff_files_list, output_file_merged)
resampling_image(output_file_merged, output_file_merged_tiff_resampled, new_resolution )

Answer 3

I edited the script from @IliasMachairas a little bit and finished with the following:

import rasterio
import os
from rasterio.merge import merge
from rasterio.warp import calculate_default_transform, reproject, Resampling

def merge_tiff(tiff_files, output_file):
    print("Merging TIFF files...")
    try:
        src_files_to_mosaic = []
        for fp in tiff_files:
            src = rasterio.open(fp)
            src_files_to_mosaic.append(src)
        
        # Merge function from Rasterio
        mosaic, out_trans = merge(src_files_to_mosaic)
        
        # Metadata of the output file
        out_meta = src_files_to_mosaic[0].meta.copy()
        out_meta.update({
            "driver": "GTiff",
            "height": mosaic.shape[1],
            "width": mosaic.shape[2],
            "transform": out_trans
        })
        
        # Save the merged TIFF
        with rasterio.open(output_file, 'w', **out_meta) as dest:
            dest.write(mosaic)
        
        # Close all the source files
        for src in src_files_to_mosaic:
            src.close()
        print(f"Merged file saved to {output_file}")
    except Exception as e:
        print(f"Error during merging: {e}")

def resample_image(input_raster_path, output_raster_path, cell_size=30):
    print("Resampling image...")
    try:
        with rasterio.open(input_raster_path) as src:
            # Compute the new transform and dimensions
            transform, width, height = calculate_default_transform(
                src.crs, src.crs, src.width, src.height, *src.bounds, resolution=cell_size)
        
            # Define the new metadata
            new_meta = src.meta.copy()
            new_meta.update({
                'transform': transform,
                'width': width,
                'height': height,
                'crs': src.crs
            })
        
            # Resample the raster
            with rasterio.open(output_raster_path, 'w', **new_meta) as dst:
                for i in range(1, src.count + 1):
                    reproject(
                        source=rasterio.band(src, i),
                        destination=rasterio.band(dst, i),
                        src_transform=src.transform,
                        src_crs=src.crs,
                        dst_transform=transform,
                        dst_crs=src.crs,
                        resampling=Resampling.bilinear)  # Change to Resampling.nearest for nearest neighbor
            print(f"Resampled file saved to {output_raster_path}")
    except Exception as e:
        print(f"Error during resampling: {e}")

def process_geotiffs(input_folder):
    tiff_files = [os.path.join(input_folder, file) for file in os.listdir(input_folder) if file.endswith('.tif')]
    if not tiff_files:
        print("No TIFF files found in the input folder.")
        return
    
    output_folder = os.path.join(input_folder, 'outputs')
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)
    
    merged_output_file = os.path.join(output_folder, 'merged_image.tif')
    resampled_output_file = os.path.join(output_folder, 'merged_image_resampled.tif')
    
    merge_tiff(tiff_files, merged_output_file)
    resample_image(merged_output_file, resampled_output_file)

if __name__ == "__main__":
    print("Starting script...")
    input_folder = os.path.dirname(os.path.abspath(__file__))
    process_geotiffs(input_folder)
    print("Script completed.")

Now you can put the skript in the same folder with the GeoTIFF's and execute it. It will create a subfolder outputs automatically and store the results inside the folder. Make sure to install rasterio in your Python Environment before you execute it.