A script I found helpful when debugging some cog issues.

generate_tiff.py

import os

import numpy as np
from osgeo import gdal, osr


def create_gradient_data(width, height, band_index, dtype=gdal.GDT_Byte):
    """Create structured gradient data for debugging rendering issues."""
    if dtype == gdal.GDT_Byte:
        max_val = 255
        np_dtype = np.uint8
    elif dtype == gdal.GDT_UInt16:
        max_val = 65535
        np_dtype = np.uint16
    else:
        max_val = 255
        np_dtype = np.uint8

    # Create coordinate grids
    x = np.linspace(0, 1, width)
    y = np.linspace(0, 1, height)
    X, Y = np.meshgrid(x, y)

    if band_index == 0:  # Red band - horizontal gradient
        data = (X * max_val).astype(np_dtype)
    elif band_index == 1:  # Green band - vertical gradient
        data = (Y * max_val).astype(np_dtype)
    elif band_index == 2:  # Blue band - diagonal gradient
        data = ((X + Y) / 2 * max_val).astype(np_dtype)
    else:  # Additional bands - radial gradient from center
        center_x, center_y = 0.5, 0.5
        distance = np.sqrt((X - center_x)**2 + (Y - center_y)**2)
        max_distance = np.sqrt(0.5)  # Maximum distance from center to corner
        data = ((1 - distance / max_distance) * max_val).astype(np_dtype)

    return data

def create_toy_cog_web_mercator(filename, width, height, num_bands, dtype=gdal.GDT_Byte):
    """Create a Cloud Optimized GeoTIFF (COG) file in EPSG:3857 (Web Mercator)."""
    temp_filename = filename.replace(".tif", "_temp_3857.tif")

    web_mercator_wkt = 'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3857"]]'

    # Define bounds directly in Web Mercator meters (EPSG:3857)
    # These coordinates are around Cuba area but in meters
    web_mercator_ulx = -9200000  # Upper left X (easting)
    web_mercator_uly = 2650000   # Upper left Y (northing)

    # Create a square 100km x 100km area
    size_meters = 100000  # 100km
    web_mercator_lrx = web_mercator_ulx + size_meters
    web_mercator_lry = web_mercator_uly - size_meters

    # Calculate pixel sizes in meters
    pixel_size_x = (web_mercator_lrx - web_mercator_ulx) / width
    pixel_size_y = (web_mercator_lry - web_mercator_uly) / height  # Will be negative

    # Step 1: Create a basic tiled, compressed GeoTIFF
    creation_options_step1 = [
        "TILED=YES",
        "COMPRESS=DEFLATE",
        "PREDICTOR=2",
        "BIGTIFF=IF_NEEDED",
        "BLOCKXSIZE=256",
        "BLOCKYSIZE=256",
    ]

    driver_gtiff = gdal.GetDriverByName("GTiff")
    if driver_gtiff is None:
        print("Error: GTiff driver not found.")
        return

    dataset = driver_gtiff.Create(
        temp_filename, width, height, num_bands, dtype, options=creation_options_step1
    )
    if dataset is None:
        print(f"Error: Could not create temporary dataset '{temp_filename}'.")
        return

    # Set GeoTransform using Web Mercator coordinates
    dataset.SetGeoTransform(
        (web_mercator_ulx, pixel_size_x, 0, web_mercator_uly, 0, pixel_size_y)
    )

    # Set projection to Web Mercator (EPSG:3857)
    dataset.SetProjection(web_mercator_wkt)

    for i in range(num_bands):
        band = dataset.GetRasterBand(i + 1)
        data = create_gradient_data(width, height, i, dtype)
        band.WriteArray(data)
        band.SetNoDataValue(0)

    # Build overviews
    dataset.BuildOverviews("AVERAGE", [2, 4, 8, 16])

    dataset = None  # Close the dataset

    # Step 2: Translate to COG format
    print(
        f"Optimizing '{temp_filename}' to COG as '{filename}' using gdal.Translate(format=\"COG\")..."
    )
    try:
        gdal.Translate(
            filename,
            temp_filename,
            format="COG",
            creationOptions=[
                "COMPRESS=DEFLATE",
                "PREDICTOR=2",
                "BLOCKSIZE=512"
            ]
        )
        print(f"Successfully created COG: {filename}")
    except Exception as e:
        print(f"Error during COG translation. Falling back to regular GeoTIFF: {e}")
        os.rename(temp_filename, filename)
    finally:
        if os.path.exists(temp_filename) and temp_filename != filename:
            os.remove(temp_filename)

if __name__ == "__main__":
    output_dir = "test_cogs_3857"
    os.makedirs(output_dir, exist_ok=True)

    print("Generating COG files in EPSG:3857 (Web Mercator)...")
    create_toy_cog_web_mercator(
        os.path.join(output_dir, "cuba_toy_cog_3band_3857.tif"), 256, 256, 3
    )
    create_toy_cog_web_mercator(
        os.path.join(output_dir, "cuba_toy_cog_1band_uint16_3857.tif"),
        512,
        512,
        1,
        dtype=gdal.GDT_UInt16,
    )
    print(
        f"Check files in '{output_dir}' with 'gdalinfo --json <file>' to confirm projection and COG compliance."
    )