Rasterio warp issue for antimeridian bounding box

[alexgleith]

-->

Expected behavior and actual behavior.

Transformation between two coordinate systems yields the correct result

Steps to reproduce the problem.

Notebook is here: https://gist.github.com/alexgleith/f94f3f3da2f43dd6cfab84eb9d8a13f7

Key code is this:

file = "https://soilspackage-useast.s3.amazonaws.com/EsriLandCover/01C_20200101-20210101.tif"
with rasterio.open(file) as f:
    bbox = f.bounds
    bbox_rasterio = warp.transform_bounds(f.crs, "epsg:4326", *f.bounds)
    
    transformer = Transformer.from_crs(f.crs, "EPSG:4326", always_xy=True)
    bbox_proj = transformer.transform(f.bounds[0], f.bounds[1])
    bbox_proj += transformer.transform(f.bounds[2], f.bounds[3])

Operating system

For example: Mac OS X 10.15.7

Rasterio version and provenance

Rasterio 1.2.6 installed on a mac with --no-binary against GDAL 3.3.0, released 2021/04/26

[snowman2]

#2144

[alexgleith]

It's not necessarily the same, @snowman2.

While the antimeridian results in coordinates that are confusing, they are usually correct.

From the notebook I linked above, coordinates from Proj and from Rasterio are as follows:

rio: -179.9831020235291 | proj: 175.7142130783133
rio: -80.51265674820361 | proj: -80.43693672756477
rio: 179.81681842082818 | proj: -170.9756575751534
rio: -77.57683559445971 | proj: -77.57683559445971

Note that both have issues, in that Proj coordinates have on X coord at 175 and the other at -170, which is painful to handle.

While the Rasterio coordinates are just plain wrong!

[snowman2]

Can you try Transformer.transform_bounds in pyproj 3.1?

[snowman2]

https://pyproj4.github.io/pyproj/stable/api/transformer.html#pyproj.transformer.Transformer.transform_bounds

[alexgleith]

Same result as before:

[snowman2]

Perfect, that is what I was hoping to see. When you run into that scenario, you will have to split the boundary by the antimeridian and end up with two bounding polygons.

[snowman2]

This is the logic to handle it:

if maxx < minx:
    return shapely.geometry.MultiPolygon(
        [
            shapely.geometry.box(minx, miny, 180, maxy),
            shapely.geometry.box(-180, miny, maxx, maxy),
        ]
    )
return shapely.geometry.box(minx, miny, maxx, maxy)

NOTE: Only for geographic coordinates.

[alexgleith]

Hmm. In the above example, the input coordinates are projected, so it's not until they're transformed that you know it crosses the antimeridian.

My solution is to do this, because I know that the features I'm working with are reasonable small:

    transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True)
    x_min, y_min = transformer.transform(bbox[0], bbox[1])
    x_max, y_max = transformer.transform(bbox[2], bbox[3])

    # Stupid dateline!
    if abs(x_max - x_min) > 90:
        if x_max < 0:
            x_max += 360
        if x_min < 0:
            x_min += 360

[snowman2]

it's not until they're transformed that you know it crosses the antimeridian.

Correct, I am suggesting to run the code above after projecting to geograpic. You can use the CRS.is_geographic property to check if the destination CRS is geographic if you aren't sure.

Side note: Transformer.transform_bounds code will handle the antimeridian for you in both directions ref.

[alexgleith]

So Proj is fine, although, handling transformed coords is still non-trivial.

I'm not sure if it's an outstanding issue in Rasterio, but the latest version installable via pip is broken, as shown above.

[snowman2]

I'm not sure if it's an outstanding issue in Rasterio, but the latest version installable via pip is broken, as shown above.

Added to rasterio #2208, targeted for 1.3 release.

[sgillies]

Thanks for the report @alexgleith, and thank you for the help @snowman2 . My understanding is that our geometry transformer rasterio.warp.transform_geom (based on GDAL's) does antimeridian cutting properly. Is that not the case?

[alexgleith]

I haven’t tested it, but I believe it’s fixed in 1.3