buildvrt: add "common" value to "-resolution" argument

[dbaston]

Feature description

I am using gdalbuildvrt to mosaic several tiles from the Copernicus 30m global DEM. The horizontal resolution of this dataset varies by tile. When using gdalbuildvrt to mosaic adjacent tiles with 3-second and 5-second horizontal resolutions, I need to use a horizontal resolution of 1 second, but none of the automatic options available with -resolution can return this value.

I am thinking of implementing -resolution=common, which would attempt to calculate a resolution that can be multiplied by an integer to arrive at both of the input resolutions. An example implementation in Python:

import pytest
import fractions
import math

def common_res(a, b):
    frac_a = fractions.Fraction(a).limit_denominator(10000)
    frac_b = fractions.Fraction(b).limit_denominator(10000)

    common_denom = math.lcm(frac_a.denominator, frac_b.denominator)

    num_a = frac_a.numerator * int(common_denom / frac_a.denominator)
    num_b = frac_b.numerator * int(common_denom / frac_b.denominator)

    common_num = math.gcd(num_a, num_b)

    return common_num / common_denom

@pytest.mark.parametrize("a,b,expected", [
    (5/3600, 3/3600, 1/3600),
    (5, 3, 1),
    (5/3600, 2.5/3600, 2.5/3600),
    (1/10, 1, 1/10),
    (1/10, 1/3, 1/30),
    (1/17, 1/13, 1/221), 
    (1/17, 1/3600, 1/61200),
    ])
def test_common_res(a, b, expected):
    assert common_res(a, b) == expected

The key piece here is the Fraction.limit_denominator method. This could be adapted from CPython or other implementations such as this one.

If this seems like an appropriate addition, I can go ahead and implement it. We would probably want some guardrails to prevent disaggregation by more than a factor of N. This also fails to consider rasters whose resolutions may be "compatible" but whose origin points make them not so.

Additional context

No response

[rouault]

sounds like a good idea

[jratike80]

Agreed, good idea. The GDIWG elevation and orthophoto products also change the horizontal resolution at some latitudes.

[rouault]

should also be mapped to new gdal CLI:

$ gdal raster buildvrt --help
Usage: gdal raster buildvrt [OPTIONS] <INPUTS> <OUTPUT>
[...]
  --resolution <xres>,<yres>|average|highest|lowest>  Target resolution (in destination CRS units)

[dbaston]

Suppose I have two 0.5-degree rasters, one with origin at longitude -180 and one with origin at longitude -179.75. Should --resolution=common result in a resolution of 0.5, 0.25, or an error?

Also, is --resolution=gcd a better name?

[rouault]

Should --resolution=common result in a resolution of 0.5, 0.25, or an error?

I'd say 0.5, but we should perhaps require that users specify a --resampling=nearest|bilinear switch then.

Also, is --resolution=gcd a better name?

No strong preference. I'd tend to find "common" slightly more intuitive. It is vague, which gives us slightly more room for being creative if we'd want to change the heuristics.

[jratike80]

"common" only slightly more intuitive than "gcd"? Not in my world.
The documentation could mension that the resulting resolution may by much smaller than what user may guess (0.75 and 1.0 -> 0.25) but gives a benefit that nearest neighbor resampling hits always just one source pixel and pixel values do not get altered. If that is the main benefit as I think.