Account for the antimeridian.

[yonda-yonda]

closes #311

rasterio transform_bounds and calculate_default_transform dosen't account for the antimeridian.
COGReader's bouds, minzoom, maxzoom are incorrect, reader.part gets an extra range.

So, I replaced transform_bound and made some minimal changes related to create tiles.

Please tell me,

• If you know a more correct way to define WEB_MERCATOR_OPP_CRS.
• If you can think of anything else that needs for this problem.
• If you find poor English comments.

(I fixed tests, but couldn't improve coverage...)

Sample Landsat-8 scenes.

• normal: LC08_L1TP_063014_20201118_20201119_01_RT
• antimeridian: LC08_L1TP_091014_20200818_20200823_01_T1

import rio_tiler
from rio_tiler.io import COGReader
from rio_tiler.utils import get_vrt_transform
from morecantile import Tile
from PIL import Image
import numpy as np
def test_landsat8(path, z, x, y):
    with COGReader(path) as cog:
        print(cog.info())
        tile, mask = cog.tile(x,y,z, tilesize=256)
        pil_img_gray = Image.fromarray((np.array(tile[0])/16).astype(np.uint8))
        pil_img_gray.save('test_tile_{}_{}_{}.jpg'.format(z,x,y))
        dst_crs = cog.tms.crs
        if cog.bounds[0] > cog.bounds[2]:
            dst_crs = cog.tile_opp_crs
        vrt_transform, vrt_width, vrt_height = get_vrt_transform(
            cog.dataset, cog.tms.xy_bounds(*Tile(x=x, y=y, z=z)), dst_crs=dst_crs
        )
        print(vrt_transform, vrt_width, vrt_height)

normal

test_landsat8("{path to}/LC08_L1TP_063014_20201118_20201119_01_RT_B3.tiff", 7,14,33)

bounds=(-140.34456198888816, 64.36592734077415, -134.56144600353784, 66.75449423805574) center=(-137.453003996213, 65.56021078941495, 5) minzoom=5 maxzoom=11 band_metadata=[('1', {})] band_descriptions=[('1', '')] dtype='uint16' nodata_type='None' colorinterp=['gray'] scale=None offset=None colormap=None
| 72.41, 0.00,-15654303.39|
| 0.00,-72.41, 9705668.10|
| 0.00, 0.00, 1.00| 4324 4324

antimeridian

test_landsat8("{path to}/LC08_L1TP_091014_20200818_20200823_01_T1_B4.tiff", 7,0,33)

bounds=(176.73667686980406, 64.40396181006827, -177.67929428770051, 66.7168416575327) center=(179.52869129105176, 65.56040173380049, 5) minzoom=5 maxzoom=11 band_metadata=[('1', {})] band_descriptions=[('1', '')] dtype='uint16' nodata_type='None' colorinterp=['gray'] scale=None offset=None colormap=None
| 72.42, 0.00,-20037508.34|
| 0.00,-72.42, 9705668.10|
| 0.00, 0.00, 1.00| 4323 4323

Before changing, get_vrt_transform return large values of Affine in this case.

| 3366.52, 0.00,-20037508.34|
| 0.00,-3366.52, 9705668.10|
| 0.00, 0.00, 1.00| 93 93

[yonda-yonda]

In tests/test_io_cogeo.py

bbox = (
        -56.624124590533825,
        73.52687881825946,
        -56.530950796449005,
        73.50183615350426,
)

Dose rio-tiler allow bbox of [minlon, maxlat, maxlon, minlat]?

In GeoJSON, bbox is [leftlon, bottomlat, ritghtlon, toplat].
https://tools.ietf.org/html/rfc7946

In some libraries, bbox is treated as [minlon, minlat, maxlon, maxlat].

What is the stance of rio-tiler?

[vincentsarago]

I'll try to have a look at this PR later this week 🙏

nice catch @yonda-yonda
the bbox should always be in form of minlon, minlat, maxlon, maxlat

ref: https://github.com/mapbox/rasterio/blob/master/rasterio/coords.py#L5

[yonda-yonda]

Thanks @vincentsarago!
Ok, I change test cases as follows.

bbox = (
        -56.624124590533825,
        73.50183615350426,
        -56.530950796449005,
        73.52687881825946,
)

And, fix other test failures too before your review.

[vincentsarago]

I'm not sure to understand why this is needed?

[yonda-yonda]

It is Opposite side CRS.
For details, please see follows.
#322 (comment)

[vincentsarago]

Not sure to understand this. I think this will make working with other TMS quite complicated 🤷‍♂️

[yonda-yonda]

I too think to be quite complicated.
But, opposite side CRS was necessary to get get_vrt_transform to calculate Affine correctly.

rio-tiler/rio_tiler/reader.py

Line 246 in 982811a

vrt_transform, vrt_width, vrt_height = get_vrt_transform(

Do you have a better idea?

[vincentsarago]

I don't have a better idea 🙁 because I yet don't fully understand the underlying need to change this (I'll comment in the issue)

[vincentsarago]

Is this the only difference with rasterio? if yes I feel we should ask Sean why not add this in rasterio directly!

[yonda-yonda]

https://github.com/mapbox/rasterio/blob/master/rasterio/warp.py#L159
In rasterio, not account for a direction of in_xs and in_ys's rotation.
So, can't decide if it's clockwise or not.

Maybe rasterio isn't interested in direction of rotation, and the antimeridian.

[vincentsarago]

To be honest @yonda-yonda, I find the solution quite complex and I yet don't fully understand the problem of rio-tiler.

I understand that:

• the bounds are not defined properly because it's crossing the antimeridian (for epsg:4326) but still rio-tiler works as expected (meaning it returns information from rasterio).
• The maxzoom/minzoom are incorrectly defined

Still, rio-tiler can create tiles from a file crossing the antimeridian. I'm not saying we shouldn't improve rio-tiler, I'm just a bit worry that we over engineer this!

[yonda-yonda]

It makes rio-tiler overcomplicated.
Keep it simple.