Irrelevant error raised during gdal.Translate

druidance commented 9 months ago

Expected behavior and actual behavior.

Expected behavior: For each folder:

all images are merged into mosaic;
the resulting mosaic is cropped with a polygon.

Actual behavior: When the loop goes for the second folder and dives into the merging algorythm, an error is raised: RuntimeError: Unknown option name '-dstnodata'.

Here's what I noticed during testing:

this error isn't raised if using CLI-style for providing Warp options:

options=f'-of {GEOTIFF_FORMAT} '
            f'-cutline {polygon_path} '
            f'-crop_to_cutline '
            f'-dstnodata {NO_DATA_VALUE} '
            f'-co BIGTIFF=YES'

the error isn't related to gdal.Translate where it's raised, but to gdal.Warp
the error isn't raised where it should be - when calling gdal.Warp
it seems the error shouldn't be raised at all, since CLI-style options work without any issue
Steps to reproduce the problem.

Prereq:
it is sufficient to have at least 2 folders each containing >= 2 images;
test images are taken from the open source 2019 NOAA NGS DSS Infrared 8 Bit Imagery: Savannah, GA but one can use any data for reproducing the issue (in this case, don't forget to provide shape files with polygons).

Specific images used:

folder src_1:
- 482500e3555000n.tif
- 482500e3557500n.tif
folder src_2:
- 482500e3560000n.tif
- 482500e3562500n.tif

import glob
import logging
import sys

from logging import getLogger
from pathlib import Path
from typing import TypeAlias

from osgeo import gdal

PATH: TypeAlias = str

GEOTIFF_FORMAT = 'GTiff'
VRT_FILE_FORMAT = '.vrt'

GEOTIFF_CREATION_OPTIONS = ['BIGTIFF=YES']

NO_DATA_VALUE = 0

logger = getLogger('RasterMerger')
logging.basicConfig(
        stream=sys.stdout,
        level=logging.INFO,
        format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',
        datefmt='%Y-%m-%d %H:%M:%S'
    )

# Turn GDAL exception ON
gdal.UseExceptions()
# Global GDAL settings to speed up all processes
gdal.SetConfigOption('GDAL_NUM_THREADS', 'ALL_CPUS')

def create_mosaic(
        src_images: list[PATH],
        mosaic_path: PATH) -> None:
    try:
        logger.info(f'Trying to create mosaic {mosaic_path} from images.')
        # Create temporary path for virtual dataset (VRT)
        vrt_path = str(Path(mosaic_path).with_suffix(VRT_FILE_FORMAT))

        # Create 'virtual dataset' of all tiles meant for merging
        # More info - https://gdal.org/programs/gdalbuildvrt.html
        gdal.BuildVRT(
            vrt_path,
            src_images
        )
        # Transform virtual dataset into GeoTIFF mosaic
        gdal.Translate(
            mosaic_path,
            vrt_path,
            format=GEOTIFF_FORMAT,
            options=GEOTIFF_CREATION_OPTIONS
        )
        logger.info(f'Successfully created mosaic {mosaic_path} from images.')

    except RuntimeError as create_err:
        logger.error(create_err)
        logger.info(f'Could not create mosaic {mosaic_path} from images.')
        raise

def crop_raster(
        polygon_path: PATH,
        src_raster_path: PATH,
        dst_raster_path: PATH) -> None:
    try:
        logger.info(f'Trying to crop {src_raster_path}...')
        gdal.Warp(
            dst_raster_path,
            src_raster_path,
            cutlineDSName=polygon_path,
            format=GEOTIFF_FORMAT,
            cropToCutline=True,
            dstNodata=NO_DATA_VALUE,
            options=GEOTIFF_CREATION_OPTIONS
        )
        logger.info(f'Successfully cropped {src_raster_path} with a polygon.')
    except RuntimeError as err:
        logger.error(err)
        logger.info(f'Could not crop raster {src_raster_path} with polygon.')
        raise

def main():
    top_folder_path = 'YOUR_TOP_FOLDER_PATH'
    src_1 = glob.glob(f'{top_folder_path}/src_1/*')
    src_2 = glob.glob(f'{top_folder_path}/src_2/*')

    for num, image_pack in enumerate((src_1, src_2)):
        mosaic_path = f'{top_folder_path}/{num}_mosaic.geotiff'
        cropped_image_path = f'{top_folder_path}/{num}_cropped.geotiff'
        polygon_path = f'{top_folder_path}/{num}_polygon.shp'

        create_mosaic(
            image_pack,
            mosaic_path
        )
        crop_raster(
            polygon_path,
            mosaic_path,
            cropped_image_path
        )

if __name__ == '__main__':
    main()

Operating system

Ubuntu 23.10 64 bit / Debian 12.4 (in Docker)

GDAL version and provenance

GDAL 3.7.1 / 3.6.2 (apt)

Python version

3.11.6

polygons_shp.zip

edit: typo

rouault commented 9 months ago

This is both a usage issue from yours and a GDAL bug:

To fix the issue on your side, you have to change options=GEOTIFF_CREATION_OPTIONS to creationOptions=GEOTIFF_CREATION_OPTIONS in both gdal.Translate() and gdal.Warp() invokations
GDAL fix in https://github.com/OSGeo/gdal/pull/9260

druidance commented 9 months ago

This is both a usage issue from yours and a GDAL bug:

To fix the issue on your side, you have to change options=GEOTIFF_CREATION_OPTIONS to creationOptions=GEOTIFF_CREATION_OPTIONS in both gdal.Translate() and gdal.Warp() invokations

GDAL fix in Python bindings: gdal.Translate()/gdal.Warp()/etc.: make sure not to … #9260

@rouault, thank you so much for providing a solution to my issue! I tested it out, and the updated code now works as expectated. You're superfast and superaccurate as always :+1:

OSGeo / gdal