Skip to content

Commit

Permalink
Change GeoBox constructor to accept shape
Browse files Browse the repository at this point in the history
`Geobox(width, height, ...)` interface is problematic for several
reasons:

1. not compatible with `ndarray`, where shape is used instead
2. no type-safety, easy to get order wrong
3. two parameters to configure single property

Replacing `width, height` with a single parameter `shape`. This could
be a simple `(int, int)` tuple for easier interop with `ndarray`. But
one can choose to use more typesafe `XY[int]`.

Example of the change at call site

```python
   # before change:
   GeoBox(512, 256, A, "espg:4326")

   # becomes
   ## if array order is "native" at call-site
   GeoBox((256, 512), A, "epsg:4326")
   ## if width/height orded is preferred
   GeoBox(ixy_(512, 256), A, "epsg:4326")
```
  • Loading branch information
Kirill888 committed Feb 4, 2022
1 parent 2aaa536 commit f8c94d9
Show file tree
Hide file tree
Showing 7 changed files with 95 additions and 84 deletions.
6 changes: 3 additions & 3 deletions odc/geo/_xr_interop.py
Expand Up @@ -280,9 +280,9 @@ def _locate_geo_info(src: XarrayObject) -> GeoState:
crs = _get_crs_from_attrs(src, sdims)

if transform is not None:
width = _xx.shape[0]
height = _yy.shape[0]
geobox = GeoBox(width, height, transform, crs)
nx = _xx.shape[0]
ny = _yy.shape[0]
geobox = GeoBox((ny, nx), transform, crs)

return GeoState(spatial_dims=sdims, transform=transform, crs=crs, geobox=geobox)

Expand Down
114 changes: 63 additions & 51 deletions odc/geo/geobox.py
Expand Up @@ -20,6 +20,7 @@
)
from .math import clamp, is_affine_st, is_almost_int
from .roi import align_up, polygon_path, roi_normalise, roi_shape
from .types import XY

# pylint: disable=invalid-name
MaybeInt = Optional[int]
Expand All @@ -45,18 +46,23 @@ class GeoBox:
Defines the location and resolution of a rectangular grid of data,
including it's :py:class:`~odc.geo.crs.CRS`.
:param shape: Shape in pixels ``(ny, nx)``
:param crs: Coordinate Reference System
:param affine: Affine transformation defining the location of the geobox
"""

def __init__(self, width: int, height: int, affine: Affine, crs: MaybeCRS):
def __init__(
self, shape: Union[Tuple[int, int], XY[int]], affine: Affine, crs: MaybeCRS
):
assert is_affine_st(
affine
), "Only axis-aligned geoboxes are currently supported"
self.width = width
self.height = height
if isinstance(shape, XY):
shape = shape.shape

self._shape = shape
self.affine = affine
self.extent = polygon_from_transform(width, height, affine, crs=crs)
self.extent = polygon_from_transform(shape[1], shape[0], affine, crs=crs)

@staticmethod
def from_geopolygon(
Expand Down Expand Up @@ -92,16 +98,16 @@ def from_geopolygon(
geopolygon = geopolygon.to_crs(crs)

bounding_box = geopolygon.boundingbox
offx, width = _align_pix(
offx, nx = _align_pix(
bounding_box.left, bounding_box.right, resolution[1], align[1]
)
offy, height = _align_pix(
offy, ny = _align_pix(
bounding_box.bottom, bounding_box.top, resolution[0], align[0]
)
affine = Affine.translation(offx, offy) * Affine.scale(
resolution[1], resolution[0]
)
return GeoBox(crs=crs, affine=affine, width=width, height=height)
return GeoBox((ny, nx), crs=crs, affine=affine)

def buffered(self, xbuff: float, ybuff: Optional[float] = None) -> "GeoBox":
"""
Expand All @@ -115,9 +121,10 @@ def buffered(self, xbuff: float, ybuff: Optional[float] = None) -> "GeoBox":
)
affine = self.affine * Affine.translation(-bx, -by)

ny, nx = (sz + 2 * b for sz, b in zip(self.shape, (by, bx)))

return GeoBox(
width=self.width + 2 * bx,
height=self.height + 2 * by,
(ny, nx),
affine=affine,
crs=self.crs,
)
Expand All @@ -137,11 +144,11 @@ def __getitem__(self, roi) -> "GeoBox":

roi = roi_normalise(roi, self.shape)
ty, tx = (s.start for s in roi)
h, w = roi_shape(roi)
ny, nx = roi_shape(roi)

affine = self.affine * Affine.translation(tx, ty)

return GeoBox(width=w, height=h, affine=affine, crs=self.crs)
return GeoBox(shape=(ny, nx), affine=affine, crs=self.crs)

def __or__(self, other) -> "GeoBox":
"""A geobox that encompasses both self and other."""
Expand All @@ -153,7 +160,7 @@ def __and__(self, other) -> "GeoBox":

def is_empty(self) -> bool:
"""Check if geobox is "empty"."""
return self.width == 0 or self.height == 0
return self._shape == (0, 0)

def __bool__(self) -> bool:
return not self.is_empty()
Expand All @@ -166,10 +173,20 @@ def transform(self) -> Affine:
"""Linear mapping from pixel space to CRS."""
return self.affine

@property
def width(self) -> int:
"""Width in pixels (nx)."""
return self._shape[1]

@property
def height(self) -> int:
"""Height in pixels (ny)."""
return self._shape[0]

@property
def shape(self) -> Tuple[int, int]:
"""Shape in pixels ``(height, width)``."""
return self.height, self.width
return self._shape

@property
def crs(self) -> Optional[CRS]:
Expand Down Expand Up @@ -206,9 +223,10 @@ def coordinates(self) -> Dict[str, Coordinate]:
"""
yres, xres = self.resolution
yoff, xoff = self.affine.yoff, self.affine.xoff
ny, nx = self.shape

xs = numpy.arange(self.width) * xres + (xoff + xres / 2)
ys = numpy.arange(self.height) * yres + (yoff + yres / 2)
xs = numpy.arange(nx) * xres + (xoff + xres / 2)
ys = numpy.arange(ny) * yres + (yoff + yres / 2)

units = self.crs.units if self.crs is not None else ("1", "1")

Expand All @@ -233,14 +251,16 @@ def __str__(self):
return f"GeoBox({self.geographic_extent})"

def __repr__(self):
return f"GeoBox({self.width}, {self.height}, {self.affine!r}, {self.crs})"
return (
f"GeoBox(({self._shape[0], self._shape[1]}), {self.affine!r}, {self.crs})"
)

def __eq__(self, other):
if not isinstance(other, GeoBox):
return False

return (
self.shape == other.shape
self._shape == other._shape
and self.transform == other.transform
and self.crs == other.crs
)
Expand Down Expand Up @@ -290,7 +310,8 @@ def bounding_box_in_pixel_domain(geobox: GeoBox, reference: GeoBox) -> BoundingB
raise ValueError("Incompatible grids")

tx, ty = round(c), round(f)
return BoundingBox(tx, ty, tx + geobox.width, ty + geobox.height)
ny, nx = geobox.shape
return BoundingBox(tx, ty, tx + nx, ty + ny)


def geobox_union_conservative(geoboxes: List[GeoBox]) -> GeoBox:
Expand All @@ -309,10 +330,7 @@ def geobox_union_conservative(geoboxes: List[GeoBox]) -> GeoBox:
)

affine = reference.affine * Affine.translation(*bbox[:2])

return GeoBox(
width=bbox.width, height=bbox.height, affine=affine, crs=reference.crs
)
return GeoBox(shape=bbox.shape, affine=affine, crs=reference.crs)


def geobox_intersection_conservative(geoboxes: List[GeoBox]) -> GeoBox:
Expand Down Expand Up @@ -342,9 +360,7 @@ def geobox_intersection_conservative(geoboxes: List[GeoBox]) -> GeoBox:

affine = reference.affine * Affine.translation(*bbox[:2])

return GeoBox(
width=bbox.width, height=bbox.height, affine=affine, crs=reference.crs
)
return GeoBox(shape=bbox.shape, affine=affine, crs=reference.crs)


def scaled_down_geobox(src_geobox: GeoBox, scaler: int) -> GeoBox:
Expand All @@ -363,13 +379,13 @@ def scaled_down_geobox(src_geobox: GeoBox, scaler: int) -> GeoBox:
"""
assert scaler > 1

H, W = (X // scaler + (1 if X % scaler else 0) for X in src_geobox.shape)
ny, nx = (X // scaler + (1 if X % scaler else 0) for X in src_geobox.shape)

# Since 0,0 is at the corner of a pixel, not center, there is no
# translation between pixel plane coords due to scaling
A = src_geobox.transform * Affine.scale(scaler, scaler)

return GeoBox(W, H, A, src_geobox.crs)
return GeoBox((ny, nx), A, src_geobox.crs)


def _round_to_res(value: float, res: float) -> int:
Expand All @@ -383,10 +399,10 @@ def flipy(gbox: GeoBox) -> GeoBox:
:returns: GeoBox covering the same region but with Y-axis flipped
"""
H, W = gbox.shape
A = Affine.translation(0, H) * Affine.scale(1, -1)
ny, _ = gbox.shape
A = Affine.translation(0, ny) * Affine.scale(1, -1)
A = gbox.affine * A
return GeoBox(W, H, A, gbox.crs)
return GeoBox(gbox.shape, A, gbox.crs)


def flipx(gbox: GeoBox) -> GeoBox:
Expand All @@ -395,10 +411,10 @@ def flipx(gbox: GeoBox) -> GeoBox:
:returns: GeoBox covering the same region but with X-axis flipped
"""
H, W = gbox.shape
A = Affine.translation(W, 0) * Affine.scale(-1, 1)
_, nx = gbox.shape
A = Affine.translation(nx, 0) * Affine.scale(-1, 1)
A = gbox.affine * A
return GeoBox(W, H, A, gbox.crs)
return GeoBox(gbox.shape, A, gbox.crs)


def translate_pix(gbox: GeoBox, tx: float, ty: float) -> GeoBox:
Expand All @@ -408,9 +424,8 @@ def translate_pix(gbox: GeoBox, tx: float, ty: float) -> GeoBox:
``(0,0)`` of the new GeoBox will be at the same location as pixel ``(tx, ty)`` in the original
GeoBox.
"""
H, W = gbox.shape
A = gbox.affine * Affine.translation(tx, ty)
return GeoBox(W, H, A, gbox.crs)
return GeoBox(gbox.shape, A, gbox.crs)


def pad(gbox: GeoBox, padx: int, pady: MaybeInt = None) -> GeoBox:
Expand All @@ -424,19 +439,20 @@ def pad(gbox: GeoBox, padx: int, pady: MaybeInt = None) -> GeoBox:

pady = padx if pady is None else pady

H, W = gbox.shape
ny, nx = gbox.shape
A = gbox.affine * Affine.translation(-padx, -pady)
return GeoBox(W + padx * 2, H + pady * 2, A, gbox.crs)
shape = (ny + pady * 2, nx + padx * 2)
return GeoBox(shape, A, gbox.crs)


def pad_wh(gbox: GeoBox, alignx: int = 16, aligny: MaybeInt = None) -> GeoBox:
"""
Expand GeoBox such that width and height are multiples of supplied number.
"""
aligny = alignx if aligny is None else aligny
H, W = gbox.shape
ny, nx = (align_up(sz, n) for sz, n in zip(gbox.shape, (aligny, alignx)))

return GeoBox(align_up(W, alignx), align_up(H, aligny), gbox.affine, gbox.crs)
return GeoBox((ny, nx), gbox.affine, gbox.crs)


def zoom_out(gbox: GeoBox, factor: float) -> GeoBox:
Expand All @@ -450,9 +466,9 @@ def zoom_out(gbox: GeoBox, factor: float) -> GeoBox:
GeoBox covering the same region but with different pixels (i.e. lower or higher resolution)
"""

H, W = (max(1, math.ceil(s / factor)) for s in gbox.shape)
ny, nx = (max(1, math.ceil(s / factor)) for s in gbox.shape)
A = gbox.affine * Affine.scale(factor, factor)
return GeoBox(W, H, A, gbox.crs)
return GeoBox((ny, nx), A, gbox.crs)


def zoom_to(gbox: GeoBox, shape: Tuple[int, int]) -> GeoBox:
Expand All @@ -462,12 +478,9 @@ def zoom_to(gbox: GeoBox, shape: Tuple[int, int]) -> GeoBox:
:returns:
GeoBox covering the same region but with different number of pixels and therefore resolution.
"""
H, W = gbox.shape
h, w = shape

sx, sy = W / float(w), H / float(h)
sy, sx = (N / float(n) for N, n in zip(gbox.shape, shape))
A = gbox.affine * Affine.scale(sx, sy)
return GeoBox(w, h, A, gbox.crs)
return GeoBox(shape, A, gbox.crs)


def rotate(gbox: GeoBox, deg: float) -> GeoBox:
Expand All @@ -482,10 +495,10 @@ def rotate(gbox: GeoBox, deg: float) -> GeoBox:
in that view arrow should point down (this is assuming usual case of inverted
y-axis)
"""
h, w = gbox.shape
c0 = gbox.transform * (w * 0.5, h * 0.5)
ny, nx = gbox.shape
c0 = gbox.transform * (nx * 0.5, ny * 0.5)
A = Affine.rotation(deg, c0) * gbox.transform
return GeoBox(w, h, A, gbox.crs)
return GeoBox(gbox.shape, A, gbox.crs)


def affine_transform_pix(gbox: GeoBox, transform: Affine) -> GeoBox:
Expand All @@ -502,9 +515,8 @@ def affine_transform_pix(gbox: GeoBox, transform: Affine) -> GeoBox:
``X_old_pix = transform * X_new_pix``
"""
H, W = gbox.shape
A = gbox.affine * transform
return GeoBox(W, H, A, gbox.crs)
return GeoBox(gbox.shape, A, gbox.crs)


class GeoboxTiles:
Expand Down
6 changes: 1 addition & 5 deletions odc/geo/gridspec.py
Expand Up @@ -123,11 +123,7 @@ def tile_geobox(self, tile_index: Tuple[int, int]) -> GeoBox:
"""
ry, rx = self.resolution
tx, ty = self._tile_txy(tile_index)
h, w = self._shape
geobox = GeoBox(
width=w, height=h, crs=self.crs, affine=Affine(rx, 0, tx, 0, ry, ty)
)
return geobox
return GeoBox(self._shape, crs=self.crs, affine=Affine(rx, 0, tx, 0, ry, ty))

def __getitem__(self, idx: Tuple[int, int]) -> GeoBox:
"""Lookup :py:class:`~odc.geo.geobox.GeoBox` of a given tile."""
Expand Down
9 changes: 5 additions & 4 deletions tests/test_gbox_ops.py
Expand Up @@ -8,6 +8,7 @@

from odc.geo import geobox as gbx
from odc.geo import geom as geometry
from odc.geo import ixy_
from odc.geo.geobox import GeoBox

# pylint: disable=pointless-statement,too-many-statements
Expand All @@ -16,7 +17,7 @@


def test_gbox_ops():
s = GeoBox(1000, 100, Affine(10, 0, 12340, 0, -10, 316770), epsg3857)
s = GeoBox(ixy_(1000, 100), Affine(10, 0, 12340, 0, -10, 316770), epsg3857)
assert s.shape == (100, 1000)

d = gbx.flipy(s)
Expand Down Expand Up @@ -147,7 +148,7 @@ def test_gbox_tiles():
A = Affine.identity()
H, W = (300, 200)
h, w = (10, 20)
gbox = GeoBox(W, H, A, epsg3857)
gbox = GeoBox(ixy_(W, H), A, epsg3857)
tt = gbx.GeoboxTiles(gbox, (h, w))
assert tt.shape == (300 / 10, 200 / 20)
assert tt.base is gbox
Expand All @@ -160,7 +161,7 @@ def test_gbox_tiles():

H, W = (11, 22)
h, w = (10, 9)
gbox = GeoBox(W, H, A, epsg3857)
gbox = GeoBox(ixy_(W, H), A, epsg3857)
tt = gbx.GeoboxTiles(gbox, (h, w))
assert tt.shape == (2, 3)
assert tt[1, 2] == gbox[10:11, 18:22]
Expand All @@ -182,7 +183,7 @@ def test_gbox_tiles():

(H, W) = (11, 22)
(h, w) = (10, 20)
tt = gbx.GeoboxTiles(GeoBox(W, H, A, epsg3857), (h, w))
tt = gbx.GeoboxTiles(GeoBox(ixy_(W, H), A, epsg3857), (h, w))
assert tt.chunk_shape((0, 0)) == (h, w)
assert tt.chunk_shape((0, 1)) == (h, 2)
assert tt.chunk_shape((1, 1)) == (1, 2)
Expand Down

0 comments on commit f8c94d9

Please sign in to comment.