fix: handle 180th meridian case for tms.tiles()

[ljstrnadiii]

👋

I ran into a strange case when my bbox for of a custom tms crossed the 180th meridian. Happy to add a test, but I now at least see expected results. Here is my example:

import pyproj
import geopandas as gpd
from pyproj import CRS
from morecantile import TileMatrixSet
from shapely.geometry import box


def get_utm_tiles(bbox) -> gpd.GeoDataFrame:
    geo = gpd.GeoDataFrame(geometry=[bbox], crs=4326)
    utm = geo.estimate_utm_crs()
    print(utm)
    rs_extent = utm.area_of_use.bounds
    tms = TileMatrixSet.custom(crs=utm, extent_crs=CRS.from_epsg(4326), extent=list(rs_extent))
    tiles = list(tms.tiles(*geo.total_bounds, zooms=4))

    # transform back in order to investigate with gpd + .explore()
    transform = pyproj.Transformer.from_crs(utm, 4326, always_xy=True)
    geos = [box(*transform.transform_bounds(*list(tms.xy_bounds(t)), densify_pts=21)) for t in tiles]
    return gpd.GeoDataFrame(geometry=geos, crs=4326)


# all relatively small; should be covered by at most 4 tiles w/ zoom 8
# works as expected
bbox1 = box(-110.80, 44.58, -110.72, 44.66)
# case: utm=32750 which generates tms.bbox that overlaps 180th meridian
bbox2 = box(119.1, -32.86, 119.2, -32.82)

print(len(get_utm_tiles(bbox1)))  # >>> 1
print(len(get_utm_tiles(bbox2)))  # >>> 150!

With the change here, we now see only 1 intersected tile and no longer see the warning that a point is not in the bounds of the tms. Happy to add this as an edge case for testing if this seems like a reasonable change! Open to other approaches as well. Would there be another place in the code base that should also handle this case?

[ljstrnadiii]

@vincentsarago thanks for unblocking tests. I'll take a closer look and think it through a bit since this change seemed to have some side effects.

[ljstrnadiii]

@vincentsarago how does this look?

[vincentsarago]

could we just do a tests to see if the w > e

[ljstrnadiii]

The specific case I am hoping to handle is when the tms.bbox lons contain the antimeridian, but also the bounds provided for tiles is on either side of the antimeridian or overlaps. This requires we swap the underlying operator (min to max or max to min on a w/e basis) in order to truly clamp and avoid generating many many tiles that don't intersect the provided bounds.

Here is an example:
Given a tms.bbox (110, 0, -158, 90) and bounds used to find tiles of (-179, 45, -178, 46), the final clamped box we want is the bounds itself since it is fully contained inside the tms.bbox i.e. (-179, 45, -178, 46).

• Assuming you mean checking if w > e for bounds only: w > e evaluates to False, and we keep w = max(self.bbox.left, w) giving us a new w of 110, but we should have -179 so that won't work (hence the need to swap max to min).
• Assuming you mean checking w > e for tms.bbox only: w > e evaluates to True, so we would swap max to min and get -179, which is what we want! However, consider the same tms.bbox, but new bounds that are completely west of the antimeridian (150, 45, 155, 46) and still contained inside tms.bbox, then applying the swapped max to min would give us 110, which is not what we want. We would actually need to keep max (not min) here to get the expected 150. That is the reason for needing to check if the anti-meridian is inside interval (tms.bbox.west, bounds.west) to swap the operator in order to truly clamp to intersection of the provided bounds and tms.bbox.

Sorry for all the words, happy to illustrate this better if the added test does not highlight this well enough.