Rasterio: warp.transform_geom not working for Point geometries

[youngpm]

Expected behavior and actual behavior.

transform_geom of a point is returning the wrong x coordinate. Works for multipoints and via transform when you feed it the coordinates directly.

Steps to reproduce the problem.

In [1]: from rasterio import crs, warp
   ...: x, y = -122.51403808499907, 38.06106733107932
   ...: print(warp.transform(crs.CRS.from_epsg(4326), crs.CRS.from_epsg(32610), [x], [y]))
   ...: print(warp.transform_geom(crs.CRS.from_epsg(4326), crs.CRS.from_epsg(32610), {'type': 'Point', 'coordinates': (x, y)}))
   ...: print(warp.transform_geom(crs.CRS.from_epsg(4326), crs.CRS.from_epsg(32610), {'type': 'MultiPoint', 'coordinates': ((x, y),)}))
   ...: 
   ...: 
([542630.876613667], [4212702.061096219])
{'type': 'Point', 'coordinates': (110.87661366700195, 4212702.061096219)}
{'type': 'MultiPoint', 'coordinates': [(542630.876613667, 4212702.061096219)]}

The x coordinate on the middle call is wrong.

Operating system

Ubuntu 18.04, python 3.6

Rasterio version and provenance

Rasterio 1.0.2

I think this has been broken for a while, just came across in some pretty old code that hasn't been ran for quite some time (a year at least!). Pretty sure it worked at one point (ha!).

[sgillies]

Confirmed.

[sgillies]

@youngpm can you confirm that the reverse transformation works?

>>> warp.transform_geom(crs.CRS.from_epsg(32610), crs.CRS.from_epsg(4326), {'type': 'Point', 'coordinates': (542630.876613667, 4212702.061096219)})
{'type': 'Point', 'coordinates': (-122.51403808499907, 38.061067331079315)}

Weird.

[youngpm]

Yeah works for me as well, bizarre! I did the reverse transform test and confirmed the bug above on osx, so it spans platforms.

[sgillies]

Here's a solution: internally, we convert all points to 1-part multipoints 😂

Just kidding.

@youngpm do you know if this affects the GDAL Python bindings?

[youngpm]

@sgillies It appears to be okay. On OSX with GDAL 2.3.1:

from osgeo import osr, ogr
geo_src = osr.SpatialReference()
geo_src.ImportFromEPSG(4326)
utm_src = osr.SpatialReference()
utm_src.ImportFromEPSG(32610)
transform = osr.CoordinateTransformation(geo_src, utm_src)
inv_transform = osr.CoordinateTransformation(utm_src, geo_src)

point = ogr.CreateGeometryFromJson('{"type":"Point", "coordinates": [-122.51403808499907, 38.06106733107932]}')
print(f'input point: {point.ExportToJson()}')
point.Transform(transform)
print(f'to utm: {point.ExportToJson()}')

Gives

input point: { "type": "Point", "coordinates": [ -122.514038084999072, 38.061067331079322 ] }
to utm: { "type": "Point", "coordinates": [ 542630.876613667001948, 4212702.061096219345927 ] }

[drnextgis]

$ git bisect bad
5d037bf102bcba78557a43efc98d442ea0e3c080 is the first bad commit
commit 5d037bf102bcba78557a43efc98d442ea0e3c080
Author: Sean Gillies <[email protected]>
Date:   Tue May 2 14:08:36 2017 +0200

    Antimeridian splitting True

If set antimeridian_cutting=False then we get an error:

In [5]: print(warp.transform_geom(crs.CRS.from_epsg(4326), crs.CRS.from_epsg(32610), {'type': 'Point', 'coordinates': (x, y)}, antimeridian_cutting=False))
---------------------------------------------------------------------------
GDALVersionError                          Traceback (most recent call last)
<ipython-input-5-d3a324217c1b> in <module>()
----> 1 print(warp.transform_geom(crs.CRS.from_epsg(4326), crs.CRS.from_epsg(32610), {'type': 'Point', 'coordinates': (x, y)}, antimeridian_cutting=False))

~/git/rasterio/rasterio/env.py in wrapper(*args, **kwds)
    361         else:
    362             with Env.from_defaults():
--> 363                 return f(*args, **kwds)
    364     return wrapper
    365 

~/git/rasterio/rasterio/env.py in wrapper(*args, **kwds)
    574                             'parameter "{0}={1}" requires '
    575                             'GDAL {2} {3}{4}'.format(
--> 576                                 param, full_kwds[param], inequality, version, reason))
    577 
    578             return f(*args, **kwds)

GDALVersionError: parameter "antimeridian_cutting=False" requires GDAL <= 2.1
Antimeridian cutting is always enabled on GDAL >= 2.2

But if remove this part of code: https://github.com/mapbox/rasterio/blob/fd412a6ea6f21f2d92e06091daad496425f690d8/rasterio/warp.py#L61-L64 then all works as expected:

In [1]: from rasterio import crs, warp

In [2]: x, y = -122.51403808499907, 38.06106733107932

In [3]: print(warp.transform_geom(crs.CRS.from_epsg(4326), crs.CRS.from_epsg(32610), {'type': 'Point', 'coordinates': (x, y)}, antimeridian_cutting=False))
{'type': 'Point', 'coordinates': (542630.8766136674, 4212702.06109646)}

[sgillies]

🙏 @drnextgis . I'm inclined to override the keyword argument and make antimeridian_cutting=False always for points – a point can't really be cut and trying to do so is clearly a problem. Afterwards, I'll report this in the GDAL tracker.

[AndrewAnnex]

I think I am seeing this issue with linestring features after upgrading to 1.0.2, oddly it happens rarely as it appears to be related to crossing the center longitude in a given projection.
My coords are going form EPSG:4326 to a sinusoidal projection centered at 350.85607029556 longitude
which is -9.143929 longitude on the first figure below.

Before Transform

After Transform

notice how the third line from the left (green dots and teal line) now is now on the far right

although looking at the plot it does not make all that much sense that only one line in that collection would be affected.

[sgillies]

@AndrewAnnex this is very helpful, thanks! Clearly, Rasterio is misusing the splitting/wrapping option of GDAL (which changed a bit in 2.2).

[AndrewAnnex]

@sgillies I partially updated my note as my idea of my observed issue being related to the center longitude does not entirely make sense to me anymore. According to that initial idea, the only line in that collection in danger of being split is the right most hockey-stick shape line (it is the closest to -9.14392..) which looks fine in the second image.

[perrygeo]

@AndrewAnnex That is strange result! Can you share the data+code to reproduce it?

[AndrewAnnex]

@perrygeo will do, probably in a few hours/tomorrow, thanks!
I should also say reverting back to 1.0a8 and python 3.6.5 produces correct results.

[AndrewAnnex]

@perrygeo here is the gist: https://gist.github.com/AndrewAnnex/417dfe3f742ecbbf0976d3b4d59fee6a
The script file plots the transformed geometry(s) on the right with the source on the left, if I switch between python 3.6.5 with 1.0a8 rasterio the correct transformation occurs, with the issue apparent when running with python 3.7 and 1.0.2 rasterio.

If you are curious about the crs I am using, I am mapping features on Mars. By convention a lot of planetary scientists / the IAU? have adopted 0-360 degree longitudes.

[perrygeo]

Thanks for the example, I can confirm the behavior you're seeing.

Proj expects data in the domain of -180 to 180 degrees. There are settings to override the wrapping behavior (https://proj4.org/usage/projections.html#longitude-wrapping) - might be worth experimenting with these.

In terms of the proj definition, I found some strange results coming from the origin longitude. Instead of using "lon_0": 350.85607029556, I tried -9.1439 to side-step the longitude wrapping issue. Same result, the middle (blue) linestring shifts right relative to the others, coordinates are x>0

But try again with "lon_0": -9.14. The result is correct. The precision of the longitude matters:

It appears that either the rasterio.warp.transform_geom function, the underlying GDAL or proj functions are numerically unstable.

@AndrewAnnex Let's move this to a separate issue since the changes in #1447 did not improve the situation and might not be related. I'll write up a ticket with more details tomorrow.

[sgillies]

If you were fetching 1.0a8 wheels from PyPI, those included GDAL 2.1.3. The 1.0a10 wheels were the first to include GDAL 2.2, which is where GDAL's geometry transformer wrapping/splitting behavior changed.

[sgillies]

So, setting the WRAPDATELINE transform option with GDAL 2.2 seems to be the culprit. I've used @AndrewAnnex's awesome gist with Rasterio 1.0.3 to make the figures below:

The green line is thrown eastward.

If we don't use WRAPDATELINE with GDAL 2.2 we get the following

Which looks correct, yes?

[sgillies]

PR #1468 undoes #1447 on the way to resolving this issue here.

[AndrewAnnex]

just got caught up with the progress on this. Thanks!