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

fix for issue 647 #648

Merged
merged 4 commits into from Oct 19, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
7 changes: 4 additions & 3 deletions rio_tiler/reader.py
Expand Up @@ -360,8 +360,6 @@ def part(
vrt_transform, vrt_width, vrt_height = get_vrt_transform(
src_dst,
bounds,
height=height,
width=width,
dst_crs=dst_crs,
)

Expand All @@ -373,9 +371,12 @@ def part(

if buffer:
bounds, height, width = _apply_buffer(buffer, bounds, height, width)

# re-calculate the transform given the new bounds, height and width
vrt_transform, vrt_width, vrt_height = get_vrt_transform(
src_dst, bounds, height, width, dst_crs=dst_crs
src_dst,
bounds,
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.

We now construct a VRT aligned with the reading bounds but which has the size defined by its projected resolution.


if padding > 0 and not is_aligned(src_dst, bounds, bounds_crs=dst_crs):
Expand Down
Binary file added tests/fixtures/mosaic_value_1.tif
Binary file not shown.
Binary file added tests/fixtures/mosaic_value_2.tif
Binary file not shown.
Binary file added tests/fixtures/mosaic_value_3.tif
Binary file not shown.
58 changes: 37 additions & 21 deletions tests/test_mosaic.py
Expand Up @@ -16,16 +16,17 @@
from rio_tiler.models import ImageData, PointData
from rio_tiler.mosaic.methods import defaults

asset1 = os.path.join(os.path.dirname(__file__), "fixtures", "mosaic_cog1.tif")
asset2 = os.path.join(os.path.dirname(__file__), "fixtures", "mosaic_cog2.tif")
asset3 = os.path.join(os.path.dirname(__file__), "fixtures", "mosaic_cog3.tif")
asset1 = os.path.join(os.path.dirname(__file__), "fixtures", "mosaic_value_1.tif")
asset2 = os.path.join(os.path.dirname(__file__), "fixtures", "mosaic_value_2.tif")
asset3 = os.path.join(os.path.dirname(__file__), "fixtures", "mosaic_value_3.tif")
Copy link
Member Author

Choose a reason for hiding this comment

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

I've updated the assets to avoid resampling to have a role (because depending on GDAL versions the pixel values can slightly change)

assets = [asset1, asset2]
assets_order = [asset2, asset1]
assets_mixed = [asset1, asset3]

stac_asset = os.path.join(os.path.dirname(__file__), "fixtures", "stac.json")

# Full covered tile
# fully covering mosaic_value_1 an partially covering mosaic_value_2
x = 150
y = 182
z = 9
Expand Down Expand Up @@ -78,7 +79,8 @@ def test_mosaic_tiler():
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8682
# Should only have value of 1
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -90,14 +92,17 @@ def test_mosaic_tiler():

img, _ = mosaic.mosaic_reader(assets, _read_tile, x, y, z, expression="b1*3")
assert img.band_names == ["b1*3"]
# Should only have value of 1 but *3
assert numpy.unique(img.data[0, m == 255]).tolist() == [3]

# Test last pixel selection
assetsr = list(reversed(assets))
(t, m), _ = mosaic.mosaic_reader(assetsr, _read_tile, x, y, z)
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8057
# Should have values of 2 and 1 because the tile do not fully overlap mosaic_value_2 cog
assert numpy.unique(t[0, m == 255]).tolist() == [1, 2]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -106,7 +111,7 @@ def test_mosaic_tiler():
assert m.shape == (256, 256)
assert t.all()
assert m.all()
assert t[0][-1][-1] == 8682
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -116,7 +121,7 @@ def test_mosaic_tiler():
)
assert len(assets_used) == 2
assert m.all()
assert t[0][-1][-1] == 8057
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -128,11 +133,12 @@ def test_mosaic_tiler():
assert mo.dtype == "uint8"

# Test brightest pixel selection
# Should return both 1 and 2
(t, m), _ = mosaic.mosaic_reader(
assets, _read_tile, x, y, z, pixel_selection=defaults.HighestMethod()
)
assert m.all()
assert t[0][-1][-1] == 8682
assert numpy.unique(t[0, m == 255]).tolist() == [1, 2]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand Down Expand Up @@ -164,7 +170,9 @@ def test_mosaic_tiler():
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8369
# This should return only 1 because we enforce data_type to be the same
# as the initial data (uint16) so 1.5 will be casted to 1
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -178,7 +186,8 @@ def test_mosaic_tiler():
pixel_selection=defaults.MeanMethod(enforce_data_type=False),
)
assert m.all()
assert t[0][-1][-1] == 8369.5
# We do not cast the data to Uint16 so we should get both 1 and 1.5
assert numpy.unique(t[0, m == 255]).tolist() == [1, 1.5]
assert t.dtype == "float64"
assert m.dtype == "uint8"

Expand All @@ -189,7 +198,9 @@ def test_mosaic_tiler():
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8369
# This should return only 1 because we enforce data_type to be the same
# as the initial data (uint16) so 1.5 will be casted to 1
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -203,7 +214,8 @@ def test_mosaic_tiler():
pixel_selection=defaults.MedianMethod(enforce_data_type=False),
)
assert m.all()
assert t[0][-1][-1] == 8369.5
# We do not cast the data to Uint16 so we should get both 1 and 1.5
assert numpy.unique(t[0, m == 255]).tolist() == [1, 1.5]
assert t.dtype == "float64"
assert m.dtype == "uint8"

Expand All @@ -219,7 +231,9 @@ def test_mosaic_tiler():
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8682
# The last band will be either 1 or 2
# so we should get both 1 and 2
assert numpy.unique(t[0, m == 255]).tolist() == [1, 2]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -235,7 +249,9 @@ def test_mosaic_tiler():
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8057
# The last band will be either 1 or 2
# but we only select where it's the lowest (1)
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand All @@ -246,7 +262,7 @@ def test_mosaic_tiler():
assert t.shape == (3, 256, 256)
assert m.shape == (256, 256)
assert m.all()
assert t[0][-1][-1] == 8682
assert numpy.unique(t[0, m == 255]).tolist() == [1]
assert t.dtype == "uint16"
assert m.dtype == "uint8"

Expand Down Expand Up @@ -432,7 +448,7 @@ def test_mosaic_tiler_with_imageDataClass():
assert img.data.shape == (3, 256, 256)
assert img.mask.shape == (256, 256)
assert img.mask.all()
assert img.data[0][-1][-1] == 8682
# assert img.data[0][-1][-1] == 8682
assert len(img.assets) == 1
assert img.crs == WEB_MERCATOR_TMS.crs
assert img.bounds == WEB_MERCATOR_TMS.xy_bounds(x, y, z)
Expand Down Expand Up @@ -486,21 +502,21 @@ def test_mosaic_point():
cog2 = (-72.02385676824944, 46.06897125935538)

pt, _ = mosaic.mosaic_point_reader(assets, _read_point, *cog1)
assert pt.data.tolist() == [8405, 8378, 9198]
assert pt.data.tolist() == [1, 1, 1]
assert pt.mask.tolist() == [255]
assert pt.assets == [asset1]
assert pt.crs == WGS84_CRS
assert pt.coordinates == cog1

pt, _ = mosaic.mosaic_point_reader(assets, _read_point, *cog2)
assert pt.data.tolist() == [9141, 9786, 9457]
assert pt.data.tolist() == [2, 2, 2]
assert pt.mask.tolist() == [255]
assert pt.assets == [asset2]
assert pt.crs == WGS84_CRS
assert pt.coordinates == cog2

pt, _ = mosaic.mosaic_point_reader(assets, _read_point, *both)
assert pt.data.tolist() == [9443, 9640, 10326]
assert pt.data.tolist() == [1, 1, 1]
assert pt.mask.tolist() == [255]
assert pt.assets == [asset1]
assert pt.crs == WGS84_CRS
Expand All @@ -509,7 +525,7 @@ def test_mosaic_point():
pt, _ = mosaic.mosaic_point_reader(
assets, _read_point, *both, pixel_selection=defaults.LowestMethod
)
assert pt.data.tolist() == [9443, 9640, 10326]
assert pt.data.tolist() == [1, 1, 1]
assert pt.mask.tolist() == [255]
assert pt.assets == [asset1, asset2]
assert pt.crs == WGS84_CRS
Expand All @@ -518,7 +534,7 @@ def test_mosaic_point():
pt, _ = mosaic.mosaic_point_reader(
assets, _read_point, *both, pixel_selection=defaults.HighestMethod
)
assert pt.data.tolist() == [9947, 10268, 10879]
assert pt.data.tolist() == [2, 2, 2]
assert pt.mask.tolist() == [255]
assert pt.assets == [asset1, asset2]
assert pt.crs == WGS84_CRS
Expand Down
68 changes: 68 additions & 0 deletions tests/test_overview.py
Expand Up @@ -138,3 +138,71 @@ def test_gcps_cog():
max_size=256, dst_crs="epsg:3857"
) # should fetch the last overview (value==3)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [3]


def test_simple_cog_tile_read():
"""Test Overview fetching with simple cog."""
# Using WarpedVRT
with Reader(COG) as src:
# Full tile
im = src.tile(173, 97, 9) # should fetch the raw resolution (value==1)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [1]

im = src.tile(87, 49, 8) # should fetch the first overview (value==2)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [2]

im = src.tile(43, 24, 7) # should fetch the second overview (value==3)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [3]

im = src.tile(21, 12, 6) # should fetch the thrid overview (value==4)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [4]

im = src.tile(10, 6, 5) # should fetch the last overview (value==5)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [5]

# Tile on border
im = src.tile(169, 102, 9) # should fetch the raw resolution (value==1)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [1]

im = src.tile(84, 48, 8) # should fetch the first overview (value==2)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [2]

im = src.tile(42, 24, 7) # should fetch the second overview (value==3)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [3]

# Using WarpedVRT with buffer
with Reader(COG) as src:
# Full tile
im = src.tile(
173, 97, 9, buffer=4
) # should fetch the raw resolution (value==1)
assert im.array.shape == (1, 264, 264)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [1]

im = src.tile(87, 49, 8, buffer=4) # should fetch the first overview (value==2)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [2]

im = src.tile(
43, 24, 7, buffer=4
) # should fetch the second overview (value==3)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [3]

im = src.tile(21, 12, 6, buffer=4) # should fetch the thrid overview (value==4)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [4]

im = src.tile(10, 6, 5, buffer=4) # should fetch the last overview (value==5)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [5]

# Tile on border
im = src.tile(
169, 102, 9, buffer=4
) # should fetch the raw resolution (value==1)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [1]

im = src.tile(84, 48, 8, buffer=4) # should fetch the first overview (value==2)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [2]

im = src.tile(
42, 24, 7, buffer=4
) # should fetch the second overview (value==3)
assert numpy.unique(im.data[0, im.mask == 255]).tolist() == [3]
40 changes: 40 additions & 0 deletions tests/test_resampling.py
@@ -1,6 +1,7 @@
"""test IO and Warp resampling."""
import os

import numpy
import pytest

from rio_tiler.io import Reader
Expand Down Expand Up @@ -51,3 +52,42 @@ def test_warp_resampling(resampling):
with Reader(COGEO) as src:
im = src.preview(max_size=64, reproject_method=resampling, dst_crs="epsg:4326")
assert im.data.any()


def test_resampling_diff():
"""Test that both `reproject` and `resampling` has influence."""
# check diff results when using different `reproject_method`
with Reader(COGEO) as src:
tile_cubic = src.tile(
43,
24,
7,
resampling_method="nearest",
reproject_method="cubic",
)
tile_nearest = src.tile(
43,
24,
7,
resampling_method="nearest",
reproject_method="nearest",
)
assert not numpy.array_equal(tile_cubic.array, tile_nearest.array)

# check diff results when using different `resampling_method`
with Reader(COGEO) as src:
tile_cubic = src.tile(
43,
24,
7,
resampling_method="cubic",
reproject_method="nearest",
)
tile_nearest = src.tile(
43,
24,
7,
resampling_method="nearest",
reproject_method="nearest",
)
assert not numpy.array_equal(tile_cubic.array, tile_nearest.array)