Opened 3 years ago

Closed 3 years ago

#18919 closed Bug (fixed)

GEOSGeometry transform method is dropping Z attribute on 3D geometries

Reported by: luizvital Owned by: nobody
Component: GIS Version: master
Severity: Normal Keywords:
Cc: Triage Stage: Accepted
Has patch: yes Needs documentation: no
Needs tests: no Patch needs improvement: no
Easy pickings: no UI/UX: no

Description

Here goes an example:
https://gist.github.com/3655316

The Z info is being lost in the conversion of GEOSGeometry to OGRGeometry, due to the use of GEOSGeometry.wkb in the constructor of OGRGeometry.

https://github.com/django/django/blob/master/django/contrib/gis/geos/geometry.py#L516

On a GEOSGeometry, only ewkb and hexewkb properties returns Z data.

But they are not understood by OGRGeometry constructor, but it accepts WKT3D.

Any hints to make GEOS to OGR 3D geometries conversion work?

Attachments (3)

18919-1.diff (3.1 KB) - added by claudep 3 years ago.
Keep 3D data when converting from GEOS to OGR
18919-2.diff (3.5 KB) - added by claudep 3 years ago.
Use more modern assertions
18919-3.diff (11.7 KB) - added by claudep 3 years ago.
Include Z value in wkb of geometries

Download all attachments as: .zip

Change History (16)

comment:1 Changed 3 years ago by DavidEklund

  • Needs documentation unset
  • Needs tests unset
  • Patch needs improvement unset

See also ticket #16594. The reason for the missing z-coordinate in the admin discussed there is due to the above problems, or related ones.

comment:2 Changed 3 years ago by claudep

  • Triage Stage changed from Unreviewed to Accepted

I've been able to reproduce. I made some encouraging tests in providing ewkb to OGRGeometry. Hopefully I can provide a patch soon.

Changed 3 years ago by claudep

Keep 3D data when converting from GEOS to OGR

comment:3 Changed 3 years ago by claudep

  • Has patch set

Testing welcome

comment:4 Changed 3 years ago by charettes

  • Patch needs improvement set

Included tests could use the assertIsInstance and assertIsNone methods for clarity.

Changed 3 years ago by claudep

Use more modern assertions

comment:5 follow-up: Changed 3 years ago by claudep

  • Patch needs improvement unset

comment:6 in reply to: ↑ 5 Changed 3 years ago by georger.silva@…

Hi guys!

I'll test this and let you know! This bug has been bothering me and the original reporter for some time :D

comment:7 Changed 3 years ago by georger.silva@…

I've just tested this and this is giving me an absurd result:

from django.contrib.gis.geos import *
p = Point(x=604567.28, y=8615808.59, z=221.14, srid=31980)
print p.x,p.y,p.z
# 604567.28 8615808.59 221.14
p.transform(4674)
print p.x
# Out[12]: -62.999999999999986
p.y
# Out[13]: -89.99999999999997
int(p.z))
# Out[15]: -41614946576139420832491006031366654402922094247483221784079256392454853820736417123190160119685530720854980944546305785406432088872556911665447459015376718790656L

This is how PostGIS processes it:

SELECT
        ST_AsEWKT(ST_Transform(ST_SetSRID(ST_MakePoint(604567.28,8615808.59,221.14),31980),4674));
-- "SRID=4674;POINT(-62.0375842938472 -12.5195091838259 221.14)"

This is a windows box, GEOS 3.3.5, python 2.7.

My linux box is broken, so I can't test that on linux. I could screwed something up while merging, but I don't think so...

comment:8 Changed 3 years ago by claudep

  • Patch needs improvement set

You are totally right, the test of my current patch is broken (Z should not be changed to 0).

comment:9 Changed 3 years ago by claudep

  • Patch needs improvement unset

While working on #16594, I realized that the 3D dimension was stripped out from the wkb of geometries. I think that difference between wkb and ewkb should only be the srid, not the Z dimension. Unfortunately, I'm unable to download the latest OGC spec to confirm this, but at least the PostGIS documentation for ST_AsBinary mentions that it supports 3D.

The patch is then more invasive that initially thought. Please test and comment!

Changed 3 years ago by claudep

Include Z value in wkb of geometries

comment:10 Changed 3 years ago by georger.silva@…

Hi Claude,

I don't think the patch needs to be too invasive. I did not test this (yet - I'll let you know in a few), but I think the solution is as simple as:

1) identify if the geometry has a z coordinate;
2) copy the points z coordinates;
3) set the z coordinates after the conversion is done;

PostGIS does not touch Z values:

{{{!python
SELECT ST_AsEWKT(ST_Transform(ST_GeomFromText('POLYGON((0 0 10,0 1 10,1 1 10,1 0 10,0 0 10))',4326),31981))

"SRID=31981;POLYGON((8234117.01325455 10000000 10,8231607.63895047 10201616.0688736 10,8432755.23239255 10206800.1025772 10,8435402.30459162 10000000 10,8234117.01325455 10000000 10))"
}}}

Nevertheless, I'll test your most welcome patch.

I'll let you know soon.

comment:11 Changed 3 years ago by luizvital

The patch worked like a charm for me!

I think the way it's implemented has a better performance than reconstruct the Z coordinates of geometries at point level.

But I really don't know if there are possible side-effects of updating thread_context.wkb_w.outdim depending on a method argument.

comment:12 Changed 3 years ago by claudep

I think it's important to find out if the Z dimension should be included in the WKB or not.

I've just been able (at last!) to download a PDF copy of the latest OpenGis Standard (Simple feature access) version 1.2.1, 2011-05-28. The WKB does include the 3D/4D (Z/M) dimensions. See also http://en.wikipedia.org/wiki/Well-known_text#Well-known_binary The original spec from 1999 did not include it, which probably explains the current implementation.

So if tests are positive (and noone opposes with good arguments), I will proceed with the current patch.

comment:13 Changed 3 years ago by Claude Paroz <claude@…>

  • Resolution set to fixed
  • Status changed from new to closed

In ffdd6595ea2220f8e8a6fb3aacd3213b751d982f:

Fixed #18919 -- Stopped dropping Z attribute when transforming geometries

Previously, the wkb of geometries was dropping the Z attribute.
Thanks luizvital for the report and tests and georger.silva@…
for the tests.

Note: See TracTickets for help on using tickets.
Back to Top