GEOSGeometry incorrectly transforms geometry.

[Yury Ryabov]

Lets demonstrate the bug with a point geometry in WGS84 Point(37 51).

Here is the plain-sql example of how PostGIS transforms a single point:

SELECT ST_Transform(ST_GeomFromEWKT('SRID=4326;Point(37 51)'), 3857) AS result
--------------------------------------------
result                                     |
-------------------------------------------|
POINT (4118821.159351122 6621293.722740169)|

Lets try to do the same with GEOSGeometry:

>>> from django.contrib.gis.geos import GEOSGeometry
>>> GEOSGeometry('srid=4326;Point(37 51)').transform(3857, clone=True).ewkt
'SRID=3857;POINT (5677294.030456952 4439106.787250583)'

The result is quite different.

Not sure if this bug belongs to Django or GEOS or it is something wrong with my setup. GEOS version is 3.7.2.

[Mariusz Felisiak]

It works for me in the same way in Django 2.2 and GEOS 3.6.2.

>>> GEOSGeometry('srid=4326;Point(37 51)').transform(3857, clone=True).ewkt
`SRID=3857;POINT (4118821.159351123 6621293.722740169)`

Can you check with the previous versions of libgeos? e.g. 3.7.1.

[Yury Ryabov]

Tried with GEOS 3.7.1 - transformation didn't even work - code crashed with exit signal 11. Tried with GEOS 3.6.2 - the result is the same for me as with GEOS 3.7.2. Tried on 2 computers with openSUSE Leap 15 on both machines.

[Claude Paroz]

Tested with Django master, GEOS 3.7.1, Debian 10, unable to reproduce the bug. I'm getting the correct result (like felixxm).

[Yury Ryabov]

Digging further into this issue I found out that GEOSGeometry.transform() method actually uses GDAL/OGR under the hood and not GEOS. On one machine I have both GDAL 2.4 and GDAL 3.0.1 installed and on the other I experimented and with available GDAL versions having only one version installed at a time. From what I observe it seems that incorrect transformation result I discovered the first time is returned when GDAL 3.0.1 is in use. If GDAL 2.4 is in use I get OGR error.

Checked versions of geospatial libraries supported by Django and it turns out that neither GDAL 3.x or 2.4 are supported. I wish Django raised exceptions or at least a warning about unsupported versions of geospatial libraries.

[Sergey Fedoseev]

This is caused by breaking changes in GDAL 3: ​https://github.com/OSGeo/gdal/blob/f6a8ae04245e00927e5d1d70ea97c32b83923614/gdal/MIGRATION_GUIDE.TXT#L7-L12

In [1]: from osgeo import gdal, ogr
   ...: from osgeo.osr import SpatialReference, CoordinateTransformation
   ...:
   ...: g = ogr.Geometry(wkt='POINT (37 51)')
   ...:
   ...: source = SpatialReference()
   ...: source.ImportFromEPSG(4326)
   ...:
   ...: dest = SpatialReference()
   ...: dest.ImportFromEPSG(3857)
   ...:
   ...: g.Transform(CoordinateTransformation(source, dest))
   ...: print(g)
   ...:
   ...: source.SetAxisMappingStrategy(0)
   ...: dest.SetAxisMappingStrategy(0)
   ...: g = ogr.Geometry(wkt='POINT (37 51)')
   ...: g.Transform(CoordinateTransformation(source, dest))
   ...: print(g)
POINT (5677294.03045695 4439106.78725058)
POINT (4118821.15935112 6621293.72274017)