Floating Point Considerations

Loophole 1: x/y coordinates precision loss

This example shows how many digits are available in both 32 and 64 bits floating point representation:

# Single-precision_floating-point_format has 24 bits of precision (~7 digits)
>>> x = np.array([12345.678901234567890123456789], dtype='float32')
>>> for i in range(5):
        print(repr(x[0]))
        x = np.nextafter(x, np.inf)
12345.678711
12345.679688
12345.680664
12345.681641
12345.682617

# Double-precision_floating-point_format has 53 bits of precision (~16 digits)
>>> x = np.array([12345.678901234567890123456789], dtype='float64')
>>> for i in range(5):
        print(repr(x[0]))
        x = np.nextafter(x, np.inf)
12345.678901234567
12345.678901234569
12345.678901234571
12345.678901234573
12345.678901234574

This example show a problematic pixel iteration using float32:

# Real data
>>> raster_top_left_y = 6826596.88515215
>>> resolution_y = -0.1442

>>> for i in range(10):
        x = np.array([raster_top_left_y + resolution_y * i], dtype='float64')[0]
        print(repr(x))
6826596.885152
6826596.740952
6826596.596752
6826596.452552
6826596.308352
6826596.164152
6826596.019952
6826595.875752
6826595.731552
6826595.587352

>>> for i in range(10):
        x = np.array([raster_top_left_y + resolution_y * i], dtype='float32')[0]
        print(repr(x))
6826597.000000
6826596.500000
6826596.500000
6826596.500000
6826596.500000
6826596.000000
6826596.000000
6826596.000000
6826595.500000
6826595.500000

Conclusion: Always use double-precision when manipulating x/y coordinates

Loophole 2: Conversion from spatial to raster

Ideally, the pixel corresponding to a spatial coordinate should be the one containing the point. This means that when converting from spatial to raster, the coordinates should be floored.

The following example highlights a loophole:

# Real data
>>> fp = Footprint(gt=(739710.21737267415, 0.14422515974562847, 0.0, 6826855.2830756251, 0.0, -0.14422515974561684), rsize=(3082, 2996))

>>> px = (10., 10.)
>>> px
(10.0, 10.0)

>>> px = fp.affine * px
>>> px
(739711.6596242716, 6826853.840824028)

>>> px = ~fp.affine * px
>>> px
(9.999999999068677, 10.0)
# Directly flooring this value would result in an error of 1 pixel in x. Instead we should round before flooring.

>>> px = tuple(np.around(px, 5)) # Round at 5 digits
>>> px
(10.0, 10.0)

>>> px = tuple(np.floor(px))
>>> px
(10.0, 10.0)

Conclusion: Always use Footprint.spatial_to_raster for your conversions