antimeridian_cutting problem

[dzanaga]

Following the discussion at geopandas/geopandas#448 I was trying the solution proposed by @snowman2 (based on fiona) to transform the Sentinel 2 epsg zone 32601 to lat lon:

epsg_number = 32601
epsg_gdf = load_epsg_gdf(epsg_number)

fig, ax = plt.subplots()

epsg_gdf.plot(ax=ax, facecolor='none', edgecolor='m')

Projecting it to lat lon using the solution proposed yields this:

which looks good overall, except for one tile (01RBK):
(zoom on the left side)

(zoom on the right)

There seems to be an issue with the projection of the bottom left point.

To reproduce:

from functools import partial

import geopandas
import fiona
from fiona.transform import transform_geom
from shapely.geometry import mapping, shape, box

%matplotlib inline

gdf = gpd.GeoDataFrame([{'tile': '01RBK', 'geometry': box(199980., 2890200.,  309780., 3000000.)}], crs={'init': 'epsg:32601'})

def base_transformer(geom, src_crs, dst_crs):
    return shape(
        transform_geom(
            src_crs=src_crs,
            dst_crs=dst_crs,
            geom=mapping(geom),
            antimeridian_cutting=True,
        )
    )

forward_transformer = partial(base_transformer, src_crs=gdf.crs, dst_crs="epsg:4326")
reverse_transformer = partial(base_transformer, src_crs="epsg:4326", dst_crs=gdf.crs)
with fiona.Env(OGR_ENABLE_PARTIAL_REPROJECTION="YES"):
    gdf2 = gdf.set_geometry(gdf.geometry.apply(forward_transformer))

Env info:

fiona                     1.8.6            py37hf242f0b_3    conda-forge
gdal                      2.4.1           py37h5f563d9_10    conda-forge
libgdal                   2.4.1               hc4f5fd6_10    conda-forge

[sgillies]

@dzanaga does this problem persist with newer versions of fiona and GDAL? Have you checked to see if this is a consequence of enabling partial reprojection?

[dzanaga]

@sgillies I still have to test this with the latest GDAL and fiona. No I haven't checked if this is caused by enabling partial reprojection. Will check it out.

[dzanaga]

The problem persists. To reproduce:

from fiona.crs import from_epsg

from functools import partial

import geopandas
import fiona
from fiona.transform import transform_geom
from shapely.geometry import mapping, shape, box
from pyproj.crs import CRS

%matplotlib inline

gdf = gpd.GeoDataFrame([{'tile': '01RBK', 'geometry': box(199980., 2890200.,  309780., 3000000.)}],
                       crs=CRS.from_epsg(32601))

def base_transformer(geom, src_crs, dst_crs):
    return shape(
        transform_geom(
            src_crs=src_crs,
            dst_crs=dst_crs,
            geom=mapping(geom),
            antimeridian_cutting=True,
        )
    )

# forward_transformer = partial(base_transformer, src_crs=gdf.crs, dst_crs=from_epsg(4326))
# reverse_transformer = partial(base_transformer, src_crs=from_epsg(4326), dst_crs=gdf.crs)

# with fiona.Env(OGR_ENABLE_PARTIAL_REPROJECTION="YES"):
#     gdf2 = gdf.set_geometry(gdf.geometry.apply(forward_transformer))
gdf['geometry'] = gdf.apply(lambda x: base_transformer(x.geometry, from_epsg(32601), from_epsg(4326)), axis=1)

fig, ax = plt.subplots()

gdf.plot(ax=ax, facecolor='lime', edgecolor='m')
ax.set_xlim(-181, -178)

Versions:

gdal                      3.0.4           py37h4b180d9_10    conda-forge
libgdal                   3.0.4               he6a97d6_10    conda-forge
fiona                     1.8.13           py37h0492a4a_1    conda-forge
pyproj                    2.6.1.post1      py37h34dd122_0    conda-forge
shapely                   1.7.0            py37hc88ce51_3    conda-forge

I also had to switch to using fiona.crs.from_epsg because using pyproj CRS.from_epsg or strings like dst_crs="epsg:32628"
was causing the operation to throw the following error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
AttributeError: 'CRS' object has no attribute 'encode'
Exception ignored in: 'fiona._transform._crs_from_crs'
AttributeError: 'CRS' object has no attribute 'encode'
---------------------------------------------------------------------------

See the commented out part of the example.

[dzanaga]

Found the workaround to the CRS issue in #580, and that works, but shouldn't that be implemented already?

[dzanaga]

For some reason, three triangles are created instead of one triangle (lon +180) and one square (lon -180). I think maybe it has to do with the fact that the upper left point is projected to exactly the antimeridian line at lon = -180 and lon = +180 which means that the polygon on the right will have only 3 vertices, and it's correctly generated, but then somehow it creates two triangles also on the left instead of one square.

'MULTIPOLYGON (
((179.9744868751185 27.08988180310208, 180 27.09042912371288, 180 26.11638182898121, 179.9744868751185 27.08988180310208)),

((-180 27.09042912371288, -179.9995670363036 26.09986130673256, -180 26.11638182898121, -180 27.09042912371288)),

((-180 27.09042912375563, -178.9188954790474 27.1093614166162, -178.9024214134282 26.1185212331391, -180 27.09042912375563))
)'

[sgillies]

@dzanaga thanks for the update. Let me see if I understand -- box(199980., 2890200., 309780., 3000000.) when reprojected to long/lat exactly touches the antimeridian?

[dzanaga]

That's what I thought but it's not like that, I projected the sinlge points:

from fiona.transform import transform_geom
from shapely.geometry import mapping, shape, box, Point
from rasterio.crs import CRS


def base_transformer(geom, src_crs, dst_crs):
    return shape(
        transform_geom(
            src_crs=src_crs.to_string(),
            dst_crs=dst_crs.to_string(),
            geom=mapping(geom),
            antimeridian_cutting=True,
        )
    )

b = box(199980., 2890200., 309780., 3000000.)

points = [Point(xy[0], xy[1]) for xy in zip(*b.exterior.coords.xy)][:-1]

points2 = [base_transformer(p, CRS.from_epsg(32601), CRS.from_epsg(4326)) for p in points]

for p1, p2 in zip(points, points2):
    print(f"p1: {p1.x}, {p1.y} - p2: {p2.x}, {p2.y}")

p1: 309780.0, 2890200.0 -> p2: -178.9024214134282, 26.1185212331391
p1: 309780.0, 3000000.0 -> p2: -178.91889547904742, 27.109361416616196
p1: 199980.0, 3000000.0 -> p2: 179.97448687511846, 27.08988180310208
p1: 199980.0, 2890200.0 -> p2: -179.9995670363036, 26.099861306732564

So it seems that only one vertex is reprojected to the other side and no points end up on the antimeridian line.

[sgillies]

@dzanaga although I don't understand why we see a problem with GDAL 2.4, I see a solution for you with GDAL 3+. Use OGC:CRS84 (strictly long/lat coordinate order) instead of EPSG:4326 (which is strictly speaking lat/long coordinate order). GDAL versions before 3.0 used to treat these as more or less the same reference system. For GDAL 3.0 they are different. We have a bug in Fiona 1.8.x (#919) such that we don't do the same thing at the antimeridian as ogr2ogr does for EPSG:4326. But if you use OGC:CRS84, which is strictly more correct, you'll get the results you expect.

from fiona.transform import transform_geom
from shapely.geometry import mapping, shape, box


def base_transformer(geom, src_crs, dst_crs):
    return shape(
        transform_geom(
            src_crs,
            dst_crs,
            geom=mapping(geom),
            antimeridian_cutting=True,
        )
    )


geom = box(199980., 2890200., 309780., 3000000.)
print(base_transformer(geom, "EPSG:32601", "OGC:CRS84"))

The result is a two part multipolygon, one part a triangular sliver, one a quadrilateral.

[sgillies]

I'm sorry it took me so long to realize that this is another manifestation of #919. Thank you so much for the data and code!

[rbuffat]

@dzanaga is this bug still present in Fiona 1.8.16?

[snowman2]

With fiona 1.8.17 the problem still persists:

Using "OGC:CRS84":

Using "EPSG:4326":

[sgillies]

Potentially related: OSGeo/gdal#2942.