Skip to content

Commit

Permalink
Round gridspec tilesize to fix floating point precision issue
Browse files Browse the repository at this point in the history
  • Loading branch information
alexgleith committed Oct 1, 2024
1 parent 84650c4 commit ff8f323
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 3 deletions.
7 changes: 5 additions & 2 deletions odc/geo/gridspec.py
Original file line number Diff line number Diff line change
Expand Up @@ -66,9 +66,12 @@ def __init__(
self.crs = norm_crs_or_error(crs)
self._shape = tile_shape
self.resolution = resolution

# This is an arbitrary rounding of 12 decimal places
# I don't know how to make it more robust
self.tile_size = xy_(
tile_shape.x * abs(resolution.x),
tile_shape.y * abs(resolution.y),
round(tile_shape.x * abs(resolution.x), 12),
round(tile_shape.y * abs(resolution.y), 12),
)
self.origin = origin

Expand Down
32 changes: 31 additions & 1 deletion tests/test_gridspec.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
import pytest
from pytest import approx

from odc.geo import CRS, res_, resyx_, xy_, yx_
from odc.geo import CRS, res_, resyx_, xy_, yx_, XY
from odc.geo.geom import polygon
from odc.geo.gridspec import GridSpec
from odc.geo.testutils import SAMPLE_WKT_WITHOUT_AUTHORITY
Expand All @@ -18,6 +18,36 @@
# pylint: disable=comparison-with-itself,unnecessary-comprehension


def test_gridspec_small():
print("Starting test for GridSpec")
WGS84GRID30 = GridSpec(
"EPSG:4326", tile_shape=(5000, 5000), resolution=0.0003, origin=XY(-180, -90)
)

assert WGS84GRID30.tile_shape == (5000, 5000)
assert WGS84GRID30.tile_size == XY(1.5, 1.5)

# Tile is at (-180 + 50*1.5) and (-90 + 50*1.5)
tile = (50, 50)
geobox = WGS84GRID30.tile_geobox(tile)
affine = geobox.affine

# Affine should be like this: (0.0003, 0, -105.0, 0, -0.0003, -13.5)
assert affine.a == 0.0003
assert affine.c == -105.0 # -180 + 50 * 1.5 * 0.0003
assert affine.f == -13.50 # -90 + 50 * 1.5 * 0.0003

# Tile is at (-180 + 200*1.5) and (-90 + 75*1.5)
tile = (200, 75)
geobox = WGS84GRID30.tile_geobox(tile)
affine = geobox.affine

# Affine should be like this: (0.0003, 0, 120.0, 0, -0.0003, 24.0)
assert affine.a == 0.0003
assert affine.c == 120.0 # -180 + 200 * 1.5 * 0.0003
assert affine.f == 24.0 # -90 + 75 * 1.5 * 0.0003


def test_gridspec():
gs = GridSpec(
crs=CRS("EPSG:4326"),
Expand Down

0 comments on commit ff8f323

Please sign in to comment.