Opened 12 years ago

Closed 12 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: dev
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 Claude Paroz 12 years ago.
Keep 3D data when converting from GEOS to OGR
18919-2.diff (3.5 KB ) - added by Claude Paroz 12 years ago.
Use more modern assertions
18919-3.diff (11.7 KB ) - added by Claude Paroz 12 years ago.
Include Z value in wkb of geometries

Download all attachments as: .zip

Change History (16)

comment:1 by David Eklund, 12 years ago

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 by Claude Paroz, 12 years ago

Triage Stage: UnreviewedAccepted

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

by Claude Paroz, 12 years ago

Attachment: 18919-1.diff added

Keep 3D data when converting from GEOS to OGR

comment:3 by Claude Paroz, 12 years ago

Has patch: set

Testing welcome

comment:4 by Simon Charette, 12 years ago

Patch needs improvement: set

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

by Claude Paroz, 12 years ago

Attachment: 18919-2.diff added

Use more modern assertions

comment:5 by Claude Paroz, 12 years ago

Patch needs improvement: unset

in reply to:  5 comment:6 by georger.silva@…, 12 years ago

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 by georger.silva@…, 12 years ago

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 by Claude Paroz, 12 years ago

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 by Claude Paroz, 12 years ago

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!

by Claude Paroz, 12 years ago

Attachment: 18919-3.diff added

Include Z value in wkb of geometries

comment:10 by georger.silva@…, 12 years ago

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 by luizvital, 12 years ago

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 by Claude Paroz, 12 years ago

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 by Claude Paroz <claude@…>, 12 years ago

Resolution: fixed
Status: newclosed

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