Inconsistent Output Dimensions from gdal.warp() with Identical outputBounds

What is the bug?

I am encountering an issue with gdal.warp() where the output dimensions are inconsistent when mosaicing two sets of images, despite using the same outputBounds and source projection (EPSG:32613) for all images, and converting them to EPSG:4326. The first set of images produces an output with one set of dimensions, while another set results in a different set of dimensions.

Steps to reproduce the issue

Define the area of interest (AOI) as a polygon with specific coordinates.
Obtain two set of images from here to create mosaic.
Use gdal.Warp with outputBounds set to the bounds of the AOI, converting dataset to EPSG:4326 on both set of datasets.

Here's the code snippet used for creating the mosaic:

from osgeo import gdal
from shapely.geometry import shape

aoi_geometry = {
                "type": "Polygon",
                "coordinates": [
                    [
                        [-106.40856573808874,35.51543260935861],[-106.40856573808874,35.10139620198477],[-105.92962613828232,35.10139620198477],[-105.92962613828232,35.51543260935861],[-106.40856573808874,35.51543260935861],
                    ]
                ]
            }

bbox = shape(aoi_geometry).bounds
sorted_mosaic_day_wise = ['/Users/Downloads/2024_03_03/T13SDU_20240303T174231_B02_10m.tif',
                          '/Users/Downloads/2024_03_03/T13SDV_20240303T174231_B02_10m.tif']

gdal.Warp(
    '/Users/Downloads/20240303_mosaiced.tif',
    sorted_mosaic_day_wise,
    options=gdal.WarpOptions(
        multithread=True,
        dstSRS="EPSG:4326",
        outputBounds=bbox,
        resampleAlg="lanczos",
        dstNodata=0,
        srcNodata=0,
        creationOptions=[
            "COMPRESS=LZW",
            "NUM_THREADS=ALL_CPUS" ],))

Expected Behavior The output mosaic should have consistent dimensions across different sets of input images when using the same outputBounds and projection. Example:- Width - 4793 Height - 4143

Actual Behavior The dimensions of the output mosaic are inconsistent between different sets of input images, even though the same outputBounds and projection conversion (EPSG:32613 to EPSG:4326) are used. The first two images result in one set of output dimensions:Width-4793 Height-4143 while the next two images produce a different set of dimensions Width-4792 Height-4143.

Versions and provenance

GDAL version: 3.7.3 Python version: 3.9.7 shapely: 2.0.1

Additional context

No response

After some digging and conversation with @rouault :

we know the target extent but not the resolution and we need to reproject.
But, starting from source rasters with different extents when we reproject the resolution can become different, so there is no way we can possibly know a common resolution starting from rasters with different extent.

I mean an "exact" solution does not exist.

a simplified version of @rouault idea to get a consistent resolution is:

take the central point C of the target extent
reproject it to the source CRS and have P1
compute P2, with P2.x = P1.x + src_gt[1] and P2.y = P1.y . Reproject P2 to the target CRS. That should give you the target horizontal resolution
do something similar with the y direction, and that should give you the target vertical resolution

A possible workaround is to specify the target resolution.

OSGeo / gdal