read_box fails if reading with the bounding box of the tiff

[jarmniku]

At least with as_crs=3857 I cannot correctly read_box() using the bounding box returned by tif_bBox_converted. This problem happens also when shrinking another dimension but keeping the other to match with the original bouding box.

The geotiff I am playing with is attached. (W5131F.tif, provided by Maanmittauslaitos 2023/05, Finland)

Code:

geotif = GeoTiff("W5131F.tif", as_crs=3857, band=0)
arr = geotif.read_box(geotif.tif_bBox_converted)
print(f"{arr.shape=}")

This should print something like 3000 by 3000 as the shape, but instead the output is:

arr.shape=(73, 3000)

I debugged for a while and noted that i_min and j_min in get_int_box() get negative values. Limiting them to be equal or greater than zero seems to fix something. Not sure if I can survive with that.

W5131F.zip

[KipCrossing]

Thanks for picking that up. If you've already debugged, please open a PR and I'll merge asap! Thanks!

[jarmniku]

I am not currently sure if my fix is the correct way to fix the issue. I am writing a solution that can automatically build up a single elevation map for a given rectangle area, using provided small geotiff tiles (6 km x 6 km). This required a bit of code logic so that the overlapping regions between the given area and the tiles are found. After the mentioned fix, the next problem seems to be that there are gaps between the tiles. It might be due to some kind of coordinate misalignment, see the image below. In this case, nine geotiff tiles are used to construct the shown elevation map.

[st-korn]

Yes, I got a similar problem!

Only 4 of 9 read_box work correctly:

I wrote a small program to test this. You can use it with any geotiff file:

from geotiff import GeoTiff

def read_square(x,y):
    arr = gtf.read_box(((x-square_size//2, y+square_size//2), (x+square_size//2, y-square_size//2)))
    print(arr.shape)

gtf = GeoTiff("ASTGTMV003_N43E039_dem.tif", as_crs=3857)
box = gtf.tif_bBox_converted
square_size = min(box[1][0]-box[0][0], box[0][1]-box[1][1]) // 3
print(box, square_size)
# ((4341444.67989728, 5465463.676723338), (4452795.092771329, 5311950.706664268)) 37116.0

# Center
read_square((box[1][0]+box[0][0])//2, (box[0][1]+box[1][1])//2)
# (869, 1199) - OK

# Right
read_square(box[1][0], (box[0][1]+box[1][1])//2)
#(869, 600) - OK

# Right-bottom
read_square(box[1][0], box[1][1])
# (438, 600) - OK

# Bottom
read_square((box[1][0]+box[0][0])//2, box[1][1])
# (438, 1199) - OK

# Left-bottom
read_square(box[0][0], box[1][1])
# (438, 0) - ISSUE

# Left
read_square(box[0][0], (box[0][1]+box[1][1])//2)
#(869, 0) - ISSUE

# Left-top
read_square(box[0][0], box[0][1])
# (0, 0) - ISSUE

# Top
read_square((box[1][0]+box[0][0])//2, box[0][1])
# (0, 1199) - ISSUE

I will try to fix it and do pull request...