transform_bounds does not appear to support bbox crossing the antimeridian

[yonda-yonda]

example (EPSG:32601)

Upper Left  (365474, 7696071) 
Lower Left  (365474, 7591951) 
Upper Right (477094, 7696071)
Lower Right (477094, 7591951)

rasterio. warp.transform_bounds return flattened bounds.

transform_bounds('EPSG:32601', 'EPSG:4326', 365474, 7591951, 477094, 7696071)
-> (-179.90689695793938, 68.40816857051695, 179.96969316611757, 69.37306968939305)

correct bounds
(179.5819541972332, 68.40816857051695, -177.55854062377162, 69.37306968939305)

It seems to need checking transformed points is crossing the antimeridian.
For example, like this.

def is_clockwise(xs, ys):
    if len(xs) != len(ys):
        raise TransformError("xs and ys arrays must be the same length") 
    
    area = 0
    if len(xs) > 2:
        if xs[0] == xs[-1] and ys[0] == ys[-1]:
            length = len(xs) - 1
        else:
            length = len(xs)
        
        for i in range(length - 1):
            area = area + xs[i] * ys[i+1] - ys[i] * xs[i+1]
        area = area + xs[length - 1] * ys[0] - ys[length - 1] * xs[0]
            
    return area >= 0

def transform_bounds(
        src_crs,
        dst_crs,
        left,
        bottom,
        right,
        top,
        densify_pts=21):
    if densify_pts < 0:
        raise ValueError('densify parameter must be >= 0')

    if src_crs == dst_crs:
        return [left, bottom, right, top]
    
    in_xs = []
    in_ys = []

    if densify_pts > 0:
        densify_factor = 1.0 / float(densify_pts + 1)
        
        # left_bottom to right_bottom
        in_xs.extend(
            left + np.arange(1, densify_pts + 1, dtype=np.float64) *
            ((right - left) * densify_factor)
        )
        in_ys.extend([bottom] * densify_pts)
        # right_bottom to right_top
        in_xs.extend([right] * (densify_pts + 2))
        in_ys.extend(
            bottom + np.arange(0, densify_pts + 2, dtype=np.float64) *
            ((top - bottom) * densify_factor)
        )
        # right_top to left_top
        in_xs.extend(
            right + np.arange(1, densify_pts + 1, dtype=np.float64) *
            ((left - right) * densify_factor)
        )
        in_ys.extend([top] * densify_pts)
        # left_top to left_bottom
        in_xs.extend([left] * (densify_pts + 2))
        in_ys.extend(
            top + np.arange(0, densify_pts + 2, dtype=np.float64) *
            ((bottom - top) * densify_factor)
        )        

    else:
        in_xs = np.array(left, right, right, left,left)
        in_ys = np.array(bottom, bottom, top, top, bottom)

    xs, ys = transform(src_crs, dst_crs, in_xs, in_ys)
    if dst_crs == 'EPSG:4326' and not is_clockwise(xs, ys):
        xs = [x + 360 if x < 0 else x for x in xs]
        return (min(xs), min(ys), max(xs) - 360, max(ys))
    
    return (min(xs), min(ys), max(xs), max(ys))

cogeo.bounds, cogeo.get_zooms, reader.part , etc., a lot of methods to create tile are affected.
Couldn't we improve?

[kylebarron]

This seems like an issue for rasterio, not rio-tiler?

[yonda-yonda]

This problem affects EPSG:4326 projection especially.
rasterio usage is not limited EPSG:4326 projection.
But Slippy map tiles always affected.
I thought this is rio-tiler's problem and should replace using rasterio.warp.transform_bounds.

If not, sorry for bothering you.

[vincentsarago]

@yonda-yonda thanks for the issue!

I will tend to think that this is a rasterio issue or a GDAL issue.

rio-tiler now uses morecantile to define it's bbox (which is now compatible with more Projection than only EPSG3857).

I've also encountered some weird behaviour of GDAL/Rasterio with coordinates outside the projection limits: https://github.com/developmentseed/morecantile/blob/master/tests/test_morecantile.py#L285-L287

We could patch rasterio.warp.transform_bounds in rio_tiler for EPSG4326 but I'm not sure this will be the best 🤷

I guess the issue is then tile_exists will raise TileOutsideBounds (or maybe not at all) 🤔

[vincentsarago]

FYI, instead of re-writting transform_bounds, we can use transform

import rasterio
from rasterio.warp import transform_bounds
from rasterio.warp import transform
print(rasterio.__gdal_version__)
>>> 2.4.4

print(transform_bounds('EPSG:32601', 'EPSG:4326', 365474, 7591951, 477094, 7696071))
>>> (-179.90689695793938, 68.40816857051695, 179.96969316611757, 69.37306968939305)


x = (365474, 477094)
y = (7591951, 7696071)
xs, ys = transform('EPSG:32601', 'EPSG:4326', x, y)
print((xs[0], ys[0], xs[1], ys[1])
>>> (179.72295317650574, 68.40816857051695, -177.58262682964732, 69.37306968939305)

Edits:

🤦 https://github.com/mapbox/rasterio/blob/master/rasterio/warp.py#L128 so yeah rasterio transform_bounds doesn't account for the antimeridian.

[yonda-yonda]

rasterio.warp.transform is good solution.
But, I think it will fail in the next case.

example (EPSG:3411)
Upper Left  (15243, -292109) 
Lower Left  (236647,-556794) 
Upper Right (279929,-75699)
Lower Right (496339,-337055)

Sample GeoJSON

{
	"type": "Polygon",
	"coordinates": [
		[
			[
				31.345,
				87.352
			],
			[
				11.171,
				84.435
			],
			[
				-21.937,
				84.433
			],
			[
				-42.182,
				87.346
			],
			[
				31.346,
				87.352
			]
		]
	]
}

rasterio.warp.transform_bound and my method

(-43.431841074419204, 83.1224363763864, 36.32837472999947, 89.28719867013199)

rasterio.warp.transform

(-43.431841074419204, 84.86156796173468, 36.32837472999947, 85.36769329194081)

Apart from the antimeridian problem, rasterio.warp.transform_bound is good work.

If improvetile_exists, how about this?

def tile_exists(self, tile_z: int, tile_x: int, tile_y: int) -> bool:
    tile = Tile(x=tile_x, y=tile_y, z=tile_z)
    tile_bounds = self.tms.bounds(*tile)
    bounds = self.bounds[:]
    if bounds[0] > bounds[2]:
        if tile_bounds[0] < 0:
            bounds[0] = -180
        elif tile_bounds[2] > 0:
            bounds[2] = 180
    return (
        (tile_bounds[0] < bounds[2])
        and (tile_bounds[2] > bounds[0])
        and (tile_bounds[3] > bounds[1])
        and (tile_bounds[1] < bounds[3])
    )

It has to be that self.bounds is always 'EPSG:4326', though.

[vincentsarago]

It has to be that self.bounds is always 'EPSG:4326', though.

this is not really defined anywhere in the docs, but yeah bounds should always be in WGS84.

Feel free to start a PR to improve tile_exits 🙏

[yonda-yonda]

I thought that's all I need are fixing tile_exits and replacing transform_bound.
But, get_zoom and reader.part are using calculate_default_transform and this method too doesn't account for the antimeridian.

Oh, dear. I'll think about it some more.

[vincentsarago]

Alright @yonda-yonda, first thanks for starting the PR I really appreciate it.

Working with bounds, crs and antimeridan is definitely not an easy task. I'd like for us to take a step back and define clearly:

• what's the problem (is there anything that make rio-tiler not working as expected)
• why we think the problem is in rio-tiler (and not in morecantile or rasterio)
• if there is a bug in rio-tiler, is it worth it to make rio-tiler more complicated for edge cases, instead of just raising a warning that rio-tiler is not optimized for this

[yonda-yonda]

Problem is incorrect bounds and maxzoom/minzoom in scenes crossing the antimeridian.

When found this problem, I thought this is a matter of how using rasterio and morecantile.
At first thought the solution would be simple, but when actually created PR, I felt the solution was getting too complicated (needing opposite side CRS, etc.).
I'm sure you're right, it makes rio-tiler overcomplicated.

We need to find more simpler solution in rio-tiler or improve rasterio/morecantile.
I guess I'll have to close PR once.