Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

generalize buffer option #519

Merged
merged 4 commits into from
Sep 1, 2022
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions rio_tiler/errors.py
Original file line number Diff line number Diff line change
@@ -13,8 +13,8 @@ class TileOutsideBounds(RioTilerError):
"""Z-X-Y Tile is outside image bounds."""


class IncorrectTileBuffer(RioTilerError):
"""Tile buffer is a float but not half of an integer"""
class InvalidBufferSize(RioTilerError):
"`buffer` must be a multiple of `0.5` (e.g: 0.5, 1, 1.5, ...)."


class PointOutsideBounds(RioTilerError):
41 changes: 16 additions & 25 deletions rio_tiler/io/cogeo.py
Original file line number Diff line number Diff line change
@@ -7,7 +7,7 @@
import attr
import numpy
import rasterio
from morecantile import BoundingBox, Tile, TileMatrixSet
from morecantile import Tile, TileMatrixSet
from rasterio import transform
from rasterio.crs import CRS
from rasterio.enums import Resampling
@@ -19,12 +19,7 @@

from .. import reader
from ..constants import WEB_MERCATOR_TMS, WGS84_CRS
from ..errors import (
ExpressionMixingWarning,
IncorrectTileBuffer,
NoOverviewWarning,
TileOutsideBounds,
)
from ..errors import ExpressionMixingWarning, NoOverviewWarning, TileOutsideBounds
from ..expression import apply_expression, get_expression_blocks, parse_expression
from ..models import BandStatistics, ImageData, Info
from ..types import BBox, DataMaskType, Indexes, NoData, NumType
@@ -312,6 +307,7 @@ def tile(
indexes: Optional[Indexes] = None,
expression: Optional[str] = None,
tile_buffer: Optional[NumType] = None,
buffer: Optional[float] = None,
**kwargs: Any,
) -> ImageData:
"""Read a Web Map tile from a COG.
@@ -323,7 +319,8 @@ def tile(
tilesize (int, optional): Output image size. Defaults to `256`.
indexes (int or sequence of int, optional): Band indexes.
expression (str, optional): rio-tiler expression (e.g. b1/b2+b3).
tile_buffer (int or float, optional): Buffer on each side of the given tile. It must be a multiple of `0.5`. Output **tilesize** will be expanded to `tilesize + 2 * tile_buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
tile_buffer (int or float, optional): Buffer on each side of the given tile. It must be a multiple of `0.5`. Output **tilesize** will be expanded to `tilesize + 2 * tile_buffer` (e.g 0.5 = 257x257, 1.0 = 258x258). DEPRECATED
buffer (float, optional): Buffer on each side of the given tile. It must be a multiple of `0.5`. Output **tilesize** will be expanded to `tilesize + 2 * tile_buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
kwargs (optional): Options to forward to the `COGReader.part` method.
Returns:
@@ -336,25 +333,12 @@ def tile(
)

tile_bounds = self.tms.xy_bounds(Tile(x=tile_x, y=tile_y, z=tile_z))
if tile_buffer is not None:
if tile_buffer % 0.5:
raise IncorrectTileBuffer(
"`tile_buffer` must be a multiple of `0.5` (e.g: 0.5, 1, 1.5, ...)."
)

x_res = (tile_bounds.right - tile_bounds.left) / tilesize
y_res = (tile_bounds.top - tile_bounds.bottom) / tilesize

# Buffered Tile Bounds
tile_bounds = BoundingBox(
tile_bounds.left - x_res * tile_buffer,
tile_bounds.bottom - y_res * tile_buffer,
tile_bounds.right + x_res * tile_buffer,
tile_bounds.top + y_res * tile_buffer,
if tile_buffer:
warnings.warn(
"`tile_buffer` is deprecated, use `buffer`.", DeprecationWarning
)

# Buffered Tile Size
tilesize += int(tile_buffer * 2)
buffer = tile_buffer

return self.part(
tile_bounds,
@@ -365,6 +349,7 @@ def tile(
max_size=None,
indexes=indexes,
expression=expression,
buffer=buffer,
**kwargs,
)

@@ -378,6 +363,7 @@ def part(
max_size: Optional[int] = None,
height: Optional[int] = None,
width: Optional[int] = None,
buffer: Optional[float] = None,
**kwargs: Any,
) -> ImageData:
"""Read part of a COG.
@@ -391,6 +377,7 @@ def part(
max_size (int, optional): Limit the size of the longest dimension of the dataset read, respecting bounds X/Y aspect ratio.
height (int, optional): Output height of the array.
width (int, optional): Output width of the array.
buffer (float, optional): Buffer on each side of the given aoi. It must be a multiple of `0.5`. Output **image size** will be expanded to `output imagesize + 2 * buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
kwargs (optional): Options to forward to the `rio_tiler.reader.part` function.
Returns:
@@ -423,6 +410,7 @@ def part(
bounds_crs=bounds_crs,
dst_crs=dst_crs,
indexes=indexes,
buffer=buffer,
**kwargs,
)

@@ -562,6 +550,7 @@ def feature(
max_size: Optional[int] = None,
height: Optional[int] = None,
width: Optional[int] = None,
buffer: Optional[NumType] = None,
**kwargs: Any,
) -> ImageData:
"""Read part of a COG defined by a geojson feature.
@@ -575,6 +564,7 @@ def feature(
max_size (int, optional): Limit the size of the longest dimension of the dataset read, respecting bounds X/Y aspect ratio.
height (int, optional): Output height of the array.
width (int, optional): Output width of the array.
buffer (int or float, optional): Buffer on each side of the given aoi. It must be a multiple of `0.5`. Output **image size** will be expanded to `output imagesize + 2 * buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
kwargs (optional): Options to forward to the `COGReader.part` method.
Returns:
@@ -601,6 +591,7 @@ def feature(
width=width,
height=height,
vrt_options=vrt_options,
buffer=buffer,
**kwargs,
)

43 changes: 40 additions & 3 deletions rio_tiler/reader.py
Original file line number Diff line number Diff line change
@@ -15,7 +15,12 @@
from rasterio.warp import transform_bounds

from .constants import WGS84_CRS
from .errors import AlphaBandWarning, PointOutsideBounds, TileOutsideBounds
from .errors import (
AlphaBandWarning,
InvalidBufferSize,
PointOutsideBounds,
TileOutsideBounds,
)
from .models import ImageData
from .types import BBox, DataMaskType, Indexes, NoData
from .utils import _requested_tile_aligned_with_internal_tile as is_aligned
@@ -136,6 +141,7 @@ def part(
minimum_overlap: Optional[float] = None,
vrt_options: Optional[Dict] = None,
max_size: Optional[int] = None,
buffer: Optional[float] = None,
**kwargs: Any,
) -> ImageData:
"""Read part of a dataset.
@@ -166,6 +172,11 @@ def part(
UserWarning,
)

if buffer and buffer % 0.5:
raise InvalidBufferSize(
"`buffer` must be a multiple of `0.5` (e.g: 0.5, 1, 1.5, ...)."
)

if bounds_crs:
bounds = transform_bounds(bounds_crs, dst_crs, *bounds, densify_pts=21)

@@ -192,8 +203,6 @@ def part(
src_dst, bounds, height, width, dst_crs=dst_crs
)

window = windows.Window(col_off=0, row_off=0, width=vrt_width, height=vrt_height)

if max_size and not (width and height):
if max(vrt_width, vrt_height) > max_size:
ratio = vrt_height / vrt_width
@@ -206,6 +215,34 @@ def part(

out_height = height or vrt_height
out_width = width or vrt_width

if buffer:
height = height or vrt_height
width = width or vrt_width

# Get output resolution
x_res = (bounds[2] - bounds[0]) / out_width
y_res = (bounds[3] - bounds[1]) / out_height

# apply buffer to bounds
bounds = (
bounds[0] - x_res * buffer,
bounds[1] - y_res * buffer,
bounds[2] + x_res * buffer,
bounds[3] + y_res * buffer,
)

# new output size
out_height += int(buffer * 2)
out_width += int(buffer * 2)

# re-calculate the transform given the new bounds, height and width
vrt_transform, vrt_width, vrt_height = get_vrt_transform(
src_dst, bounds, out_height, out_width, dst_crs=dst_crs
)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This approach will be a bit slower because of the call to get_vrt_transform but it will produce better result

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I'm not sure there is a way around this because we need to run get_vrt_transform first before doing these buffer calculations 🤷


window = windows.Window(col_off=0, row_off=0, width=vrt_width, height=vrt_height)

if padding > 0 and not is_aligned(src_dst, bounds, out_height, out_width, dst_crs):
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just a FYI:
padding is different than buffer because it adds pixel around the VRT while buffer adds pixel around the output image

vrt_transform = vrt_transform * Affine.translation(-padding, -padding)
orig_vrt_height = vrt_height
14 changes: 7 additions & 7 deletions tests/test_io_cogeo.py
Original file line number Diff line number Diff line change
@@ -20,7 +20,7 @@
from rio_tiler.errors import (
AlphaBandWarning,
ExpressionMixingWarning,
IncorrectTileBuffer,
InvalidBufferSize,
NoOverviewWarning,
TileOutsideBounds,
)
@@ -199,7 +199,7 @@ def test_tile_valid_default():
# We are using a file that is aligned with the grid so no resampling should be involved
with COGReader(COG_WEB) as cog:
img = cog.tile(147, 182, 9)
img_buffer = cog.tile(147, 182, 9, tile_buffer=10)
img_buffer = cog.tile(147, 182, 9, buffer=10)
assert img_buffer.width == 276
assert img_buffer.height == 276
assert not img.bounds == img_buffer.bounds
@@ -214,26 +214,26 @@ def test_tile_invalid_bounds():


def test_tile_with_incorrect_float_buffer():
with pytest.raises(IncorrectTileBuffer):
with pytest.raises(InvalidBufferSize):
with COGReader(COGEO) as cog:
cog.tile(43, 24, 7, tile_buffer=0.8)
cog.tile(43, 24, 7, buffer=0.8)


def test_tile_with_int_buffer():
with COGReader(COGEO) as cog:
data, mask = cog.tile(43, 24, 7, tile_buffer=1)
data, mask = cog.tile(43, 24, 7, buffer=1)
assert data.shape == (1, 258, 258)
assert mask.all()

with COGReader(COGEO) as cog:
data, mask = cog.tile(43, 24, 7, tile_buffer=0)
data, mask = cog.tile(43, 24, 7, buffer=0)
assert data.shape == (1, 256, 256)
assert mask.all()


def test_tile_with_correct_float_buffer():
with COGReader(COGEO) as cog:
data, mask = cog.tile(43, 24, 7, tile_buffer=0.5)
data, mask = cog.tile(43, 24, 7, buffer=0.5)
assert data.shape == (1, 257, 257)
assert mask.all()

55 changes: 55 additions & 0 deletions tests/test_reader.py
Original file line number Diff line number Diff line change
@@ -387,3 +387,58 @@ def test_point():
# Test with COG + Alpha Band
assert not reader.point(src_dst, [-104.77519499, 38.95367054])[0]
assert reader.point(src_dst, [-104.77519499, 38.95367054], masked=False)[0] == 0


def test_part_with_buffer():
"""Make sure buffer works as expected."""
bounds = [
-6574807.42497772,
12210356.646387195,
-6261721.357121638,
12523442.714243278,
]
# Read part at full resolution
with rasterio.open(COG) as src_dst:
img_no_buffer = reader.part(src_dst, bounds, dst_crs=constants.WEB_MERCATOR_CRS)

x_size = img_no_buffer.width
y_size = img_no_buffer.height

x_res = (bounds[2] - bounds[0]) / x_size
y_res = (bounds[3] - bounds[1]) / y_size

nx = x_size + 4
ny = y_size + 4

# apply a 2 pixel buffer
bounds_with_buffer = (
bounds[0] - x_res * 2,
bounds[1] - y_res * 2,
bounds[2] + x_res * 2,
bounds[3] + y_res * 2,
)
with rasterio.open(COG) as src_dst:
img = reader.part(
src_dst,
bounds_with_buffer,
height=ny,
width=nx,
dst_crs=constants.WEB_MERCATOR_CRS,
)
assert img.width == nx
assert img.height == ny

with rasterio.open(COG) as src_dst:
imgb = reader.part(
src_dst, bounds, buffer=2, dst_crs=constants.WEB_MERCATOR_CRS
)
assert imgb.width == nx
assert imgb.height == ny

assert numpy.array_equal(img.data, imgb.data)
assert img.bounds == imgb.bounds

# No resampling is involved. Because we read the full resolution data
# all arrays should be equal
numpy.array_equal(img_no_buffer.data, imgb.data[:, 2:-2, 2:-2])
numpy.array_equal(img_no_buffer.data, img.data[:, 2:-2, 2:-2])
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here we make sure that reading with buffer (manually created and using buffer=) give the same result as non-buffer read