opendatacube / odc-geo

GeoBox and geometry utilities extracted from datacube-core
https://odc-geo.readthedocs.io/en/latest/
Apache License 2.0
77 stars 12 forks source link

Differences in shape using from_geopolygon with odc.geo.GeoBox and datacube.utils.geometry.GeoBox #130

Open snowman2 opened 6 months ago

snowman2 commented 6 months ago

With the datacube 1.9 release coming soon, I wanted to make sure that some minor differences are expected between the behavior of the GeoBox classes. This issue will also serve as a reference for others who encounter this.

Library versions:

datacube 1.9.0rc1 & 1.8.17
odc-geo 0.4.2

Base setup:

from odc.geo.geom import Geometry
from odc.geo.geobox import GeoBox as OdcGeoBox
from odc.geo import res_
from datacube.utils.geometry import GeoBox as DcGeoBox

geom = Geometry(
    geom={
        "type": "Polygon",
        "coordinates": (
            (
                (-93.3043568, 41.8665694),
                (-93.3043349, 41.8701339),
                (-93.3008453, 41.8701405),
                (-93.2985946, 41.8690252),
                (-93.2979561, 41.8685513),
                (-93.2985067, 41.8665745),
                (-93.3043568, 41.8665694),
            ),
        ),
    },
    crs="EPSG:4326",
)
dst_crs="EPSG:6933"

odc.geo:

odc_geo_box = OdcGeoBox.from_geopolygon(
    geopolygon=geom,
    resolution=res_(3),
    crs=dst_crs,
)
odc_geo_box.shape.yx
(115, 207)
odc_geo_box = OdcGeoBox.from_geopolygon(
    geopolygon=geom,
    resolution=res_(3),
    crs=dst_crs,
    tight=True
)
odc_geo_box.shape.yx
(114, 206)

datacube.utils.geometry:

dc_geo_box = DcGeoBox.from_geopolygon(
    geopolygon=geom,
    resolution=[-3, 3],
    crs=dst_crs,
)
dc_geo_box.shape
(114, 207)

I am going to dig into the history

snowman2 commented 6 months ago

Difference from #38 released in v0.1.2. Specifically this commit 16714c29 according to git bisect:

commit 16714c29aa88402eea93aae0565929688fe27338
Author: Kirill Kouzoubov <kirill888@gmail.com>
Date:   Fri Apr 22 20:59:55 2022 +1000

    Replace alignment with anchor

    `align=` was too confusing and hard to document, concept of anchor I
    feel is more reasonable. You can anchor pixel edges to 0,0 or pixel
    center to 0,0, or you can leave it floating. Or you can pick any
    location within a pixel and say that 0,0 should coincide with that.

    center == (0.5, 0.5)
    edge == (0, 0)

    Also replacing _align_pix methods, that were hard to decipher, with
    snap_grid method in odc.geo.math
snowman2 commented 6 months ago
from odc.geo.math import snap_grid

bbox = geom.to_crs(dst_crs).boundingbox
tol = 0.01
offx, nx = snap_grid(bbox.left, bbox.right, rx, None, tol=tol)
offy, ny = snap_grid(bbox.bottom, bbox.top, ry, None, tol=tol)
(ny, nx), (offx, offy)
((114, 206), (-9002590.31883444, 4888352.245956851))
tol = 0.01
offx, nx = snap_grid(bbox.left, bbox.right, rx, 0, tol=tol)
offy, ny = snap_grid(bbox.bottom, bbox.top, ry, 0, tol=tol)
(ny, nx), (offx, offy)
((115, 207), (-9002592.0, 4888353.0))
from datacube.utils.geometry._base import _align_pix

offx, nx = _align_pix(bbox.left, bbox.right, rx, 0)
offy, ny = _align_pix(bbox.bottom, bbox.top, ry, 0)
(ny, nx), (offx, offy)
((114, 207), (-9002592.0, 4888353.0))
SpacemanPaul commented 6 months ago

There are some minor differences between odc-geo GeoBox and datacube.utils GeoBox, yes.

1.9 will use the odc-geo GeoBox internally, but still provide the old datacube.utils class with a deprecation warning.

Kirill888 commented 6 months ago

I'm not surprised that difference exists. I guess a more important question is which one is "more correct", it's good to have such definition, ideally backed by rigorous tests.

SpacemanPaul commented 6 months ago

IMO the way odc-geo provides clear, consistent for its behaviour in potentially ambiguous areas - AND provides convenient tools for interoperating with systems that made different choices those areas - is the strongest argument for migrating to odc-geo.

Core's current behaviour is unclear, inconsistent and poorly documented.