gdal.BuildVRT incorrectly changes the resolution (with resolution="average")

The gdal.Dataset returned from gdal.BuildVRT changes both the x- and y-resolution based on the number of inputs provided, even if all inputs have the exact same resolution.

Some things I have observed so far:

It happens only with the default resolution="average", both "highest" and "lowest" avoid it.
It happens both for the Python interface and the command-line utility (Library\bin\gdalbuildvrt.exe)
It's independent of stack or mosaic (separate = True/False)
It already happens on the gdal.Dataset, before writing the VRT (but it does also end up the in VRT).
It happens both for the geotiff and netcdf4 formats as the input (either on-disk or in-mem)

Based on the example below, using the standard sum function in Python seems to replicate it exactly:
assert gt_ds[1] == sum([reference_gt[1]]*n)/n

Expected behavior and actual behavior.

The resolution returned from gdal.BuildVRT should not change (for identical input resolutions).

Steps to reproduce the problem.

The example below has two different geotransforms (one commented out) for which it both occurs. The geotransform of ~3000 meters has much larger differences compared to the ~3km one.

Both are real uses cases, being the geotransform for the original grid of Meteosat Second Generation, The example code below doesn't write any real data (didn't seem to matter), but an actual input could be downloaded from the url below:

https://msgcpp-adaguc.knmi.nl/adaguc-server?dataset=msgrt&service=wcs&request=getcoverage&coverage=lwe_precipitation_rate&FORMAT=NetCDF4&time=current

(which was also the source of the -0.1 nodata from #7486 btw)

from osgeo import gdal

def create_test_file(fn="/vsimem/tmp_file123.tif"):

    # https://msgcpp.knmi.nl/data-access.html
    gt = (-5570.248248450553, 3.0004031511048237, 0.0, 5570.248248450553, 0.0, -3.0004031511048237)

    # # https://github.com/pytroll/satpy/blob/main/satpy/etc/areas.yaml
    # gt = (-5570248.686685662, 3000.4032785810186, 0.0, 5570248.686685662, 0.0, -3000.4032785810186)

    drv = gdal.GetDriverByName("GTiff")

    ds = drv.Create(fn, 1, 1, 1, gdal.GDT_Float32)
    ds.SetGeoTransform(gt)
    ds = None

    return fn

fn = create_test_file()

ds = gdal.Open(fn)
reference_gt = ds.GetGeoTransform()
print(reference_gt)
ds = None

for n in [1, 13, 15, 18]:

    fn_tmp = f"/vsimem/tmp_vrt_{n}"
    ds = gdal.BuildVRT(fn_tmp, [fn]*n, separate=False)
    gt_ds = ds.GetGeoTransform() # from the Dataset
    ds = None

    print(f"{n=:<4d}:", gt_ds)

    gdal.Unlink(fn_tmp)

Which for me outputs:

(-5570.248248450553, 3.0004031511048237, 0.0, 5570.248248450553, 0.0, -3.0004031511048237)
n=1   : (-5570.248248450553, 3.0004031511048237, 0.0, 5570.248248450553, 0.0, -3.0004031511048237)
n=13  : (-5570.248248450553, 3.000403151104824,  0.0, 5570.248248450553, 0.0, -3.000403151104824)
n=15  : (-5570.248248450553, 3.0004031511048246, 0.0, 5570.248248450553, 0.0, -3.0004031511048246)
n=18  : (-5570.248248450553, 3.000403151104825,  0.0, 5570.248248450553, 0.0, -3.000403151104825)

I checked that example for up to 1000 inputs, which results in the following deviations from the input:

Using the ~3km (KNMI) geotransform:

gdal_vrt_input_amount_difference_knmi

Using the ~3000m (Satpy) geotransform, note the different y-limits:

gdal_vrt_input_amount_difference_satpy

As mentioned above sum([n_values])/n replicates it, so it could be related to a summation taking place when calculating the average resolution. An example:

import numpy as np
import math

n = 73
values = [3000.4032785810186] * n

sum(values) / n              # 3000.403278581014 (wrong)
np.asarray(values).mean()    # 3000.403278581019 (wrong)
np.asarray(values).sum() / n # 3000.403278581019 (wrong)
math.fsum(values) / n        # 3000.4032785810186 (correct)

IIRC Numpy uses pairwise summation internally. https://en.wikipedia.org/wiki/Kahan_summation_algorithm

Operating system

Windows 10 Pro (64bit) 22H2

GDAL version and provenance

Windows 11 64bit
python                 3.10.9          h4de0772_0_cpython   conda-forge
gdal                   3.6.2           py310h644bc08_9      conda-forge
libgdal                3.6.2           h221c0cb_9           conda-forge
numpy                  1.23.5          py310h4a8f9c9_0      conda-forge

There is discussion about fast and almost-as-accurate-as-possible mean "fmean" in https://bugs.python.org/issue35904. Out of curiosity, does the difference in the 12th decimal of the resolution make practical issues for you?

In terms of results it wouldn't matter to me, in a case like this I might for example use the coordinates to calculate solar angles, for which such an accuracy isn't necessary at all.

It's more an issue (annoyance) in terms of workflow for me. I usually either religiously check that all spatial properties match before combining different sources, or throw an Exception. Or when I expect differences I might automatically pass it through gdal.Wap for example, which also isn't great for a case like this. And inheriting output properties from inputs like this causes it to potentially percolate to downstream outputs.

Now that I'm aware, I can override the default "average" for cases where I know upfront everything should be the same anyway, which is the bulk of my use cases. It took a while before I realized where it originated, which is part of posting it here, it could confirm others noticing similar results.

I appreciate the rabbit hole that's floating point precision, so if it's difficult to avoid in this case I can live with it. Thanks for the link, Raymond Hettinger has done some great work (and talks) on these topics.

I'm guessing this doesn't happen in Python (but C/C++), so it won't be as simple replacing a sum with fsum. 😅

OSGeo / gdal