3d support for PostGIS geometries?

[Tom Kralidis]

Hi all: does GeoAlchemy2 support 3d geometries? Requirement is managing 3d point geometries (create, insert, select).

Thanks

..Tom

[Eric Lemoine]

On Tue, Jul 25, 2017 at 8:50 PM, Tom Kralidis <tom.kr...@gmail.com> wrote:
> Hi all: does GeoAlchemy2 support 3d geometries? Requirement is managing 3d
> point geometries (create, insert, select).

Yes, GeoAlchemy2 should correctly handle 3D geometries. Would be good
to know more about your use-cases… Feel free to create GitHub issues
if you run into problems.

--
Eric

[Tom Kralidis]

Eric: thanks for the info and confirmation.  I didn't see mention in the docs or examples; I'll open a ticket.

Thanks

..Tom

[Tom Kralidis]

My specific use case is creating a Geometry type that can accept both POINT and POINTZ geometries.

The following example (https://gist.github.com/tomkralidis/3f11c88969e0c2d9d7a213827d75ab56) returns:

sqlalchemy.exc.DataError: (psycopg2.DataError) Geometry has Z dimension but column does not

So the insert works for the 2D geom but not the 3D geom.  This also happens when I add `dimension=2`
to the Geometry definition.

Any ideas?

Thanks

..Tom




On Wednesday, July 26, 2017 at 3:16:16 AM UTC-4, erilem wrote:

[Eric Lemoine]

On Wed, Jul 26, 2017 at 8:42 PM, Tom Kralidis <tom.kr...@gmail.com> wrote:
> My specific use case is creating a Geometry type that can accept both POINT
> and POINTZ geometries.
>
> The following example
> (https://gist.github.com/tomkralidis/3f11c88969e0c2d9d7a213827d75ab56)
> returns:
>
> sqlalchemy.exc.DataError: (psycopg2.DataError) Geometry has Z dimension but
> column does not
>
> So the insert works for the 2D geom but not the 3D geom. This also happens
> when I add `dimension=2`
> to the Geometry definition.
>
> Any ideas?

What you get here is a PostGIS error. It's not related to GeoAlchemy.

This is how you can reproduce the issue with SQL:

create table test (id int, geom geometry(geometry, 4326));
insert into test values (1, st_geomfromewkt('SRID=4326;pointz (1 2 3)'));

I guess you need to use a pointz column. So:

geom = Geometry(geometry_type='POINTZ', srid=4326)

And you'd insert pointz with z set to 0 for 2d points.

Would that work for you?

--
Eric

[Tom Kralidis]

That would work (i.e not throw error) but putting a 0 value could lead to the interpretation
that the elevation is actually 0 (as opposed to null).

I've added a psql equivalent of what I'm looking for in https://gist.github.com/tomkralidis/3f11c88969e0c2d9d7a213827d75ab56#file-test_zgeom-sql

Is there a way to achieve the same via GeoAlchemy?

Cheers

..Tom

[Eric Lemoine]

On Wed, Jul 26, 2017 at 11:25 PM, Tom Kralidis <tom.kr...@gmail.com> wrote:
> That would work (i.e not throw error) but putting a 0 value could lead to
> the interpretation
> that the elevation is actually 0 (as opposed to null).
>
> I've added a psql equivalent of what I'm looking for in
> https://gist.github.com/tomkralidis/3f11c88969e0c2d9d7a213827d75ab56#file-test_zgeom-sql
>
> Is there a way to achieve the same via GeoAlchemy?

Not sure.

One thing to note:

In your SQL test case you used the following to create the geometry column:

create table testptz (pkey integer primary key, geom geometry);

And you can then use this to insert a 3D point:

insert into testptz values (0, 'POINT(0 0 0)');

If you create the geometry column using "geom geometry(geometry, 0)",
which is what GeoAlchemy would do, the insertion does not work
anymore:

create table testptz (pkey integer primary key, geom geometry(geometry, 0));
insert into testptz values (1, 'POINT(0 0 0)');
ERROR: Geometry has Z dimension but column does not

I am not sure why "geom geometry" and "geom geometry(geometry, 0)"
behave differently.


--
Eric

[Tobias]

Hi Tom,

here is an example where GeoAlchemy 2 is used to store geometries of different types and dimensions:
https://github.com/c2corg/v6_api/blob/596440e/c2corg_api/models/document.py#L309-L320

The trick is to use `management=True` and `use_typmod=False` to avoid that a fixed dimension is enforced.

Other gotchas: If you want geometries with more than 3 dimensions (for example tracks with elevation and time), you can no longer use Shapely. In this project we are using geomet. In that case you will also have to force GeoAlchemy to use WKB instead of EWKB: https://github.com/c2corg/v6_api/blob/596440e/c2corg_api/models/document.py#L286-L287

I hope this helps a bit.

Tobias

[Eric Lemoine]

On Thu, Jul 27, 2017 at 9:55 AM, Tobias <to....@gmail.com> wrote:
> Hi Tom,
>
> here is an example where GeoAlchemy 2 is used to store geometries of
> different types and dimensions:
> https://github.com/c2corg/v6_api/blob/596440e/c2corg_api/models/document.py#L309-L320
>
> The trick is to use `management=True` and `use_typmod=False` to avoid that a
> fixed dimension is enforced.
>
> Other gotchas: If you want geometries with more than 3 dimensions (for
> example tracks with elevation and time), you can no longer use Shapely. In
> this project we are using geomet. In that case you will also have to force
> GeoAlchemy to use WKB instead of EWKB:
> https://github.com/c2corg/v6_api/blob/596440e/c2corg_api/models/document.py#L286-L287
>
> I hope this helps a bit.

Thanks Tobias. That makes some sense.



That being said, I don't undertand why this works:

create table testptz (pkey integer primary key, geom geometry);

insert into testptz values (0, 'POINT(0 0 0)');

while this does not:

create table testptz (pkey integer primary key, geom geometry(geometry, 0));

insert into testptz values (0, 'POINT(0 0 0)');

In both cases typmod is used.


--
Eric

[Tom Kralidis]

Thanks Tobias.  I've updated my gist example with your suggestion (adding management=True, use_typmod=False in https://gist.github.com/tomkralidis/3f11c88969e0c2d9d7a213827d75ab56#file-test_zgeom-py-L30).  The 2D point now inserts correctly but the 3D point does not and still results in error:

sqlalchemy.exc.IntegrityError: (psycopg2.IntegrityError) new row for relation "data_record_foo" violates check constraint "enforce_dims_location"
DETAIL:  Failing row contains (2, 01010000A0E61000000000000000C052C000000000008046400000000000C072...).
 [SQL: 'INSERT INTO data_record_foo (location) VALUES (ST_GeomFromEWKT(%(location)s)) RETURNING data_record_foo.identifier'] [parameters: {'location': 'SRID=4326;POINTZ(-75 45 300)'}]

Any idea on what I'm doing wrong here?

Thanks

..Tom

[Tobias]

Sorry, I missed this part:
https://github.com/c2corg/v6_api/blob/494d48d/alembic_migration/versions/38df9393c9a9_init.py#L1436-L1439

We are explicitly removing the constraints.

[Tom Kralidis]

Thanks Tobias. It would be great to have the same functionality as achieved by:

create table testptz (pkey integer primary key, geom geometry);

insert into testptz values (0, 'POINT(0 0 0)');

Should I open a GitHub issue?

[Eric Lemoine]

On Thu, Jul 27, 2017 at 2:05 PM, Tom Kralidis <tom.kr...@gmail.com> wrote:
> Thanks Tobias. It would be great to have the same functionality as achieved
> by:
>
> create table testptz (pkey integer primary key, geom geometry);
> insert into testptz values (0, 'POINT(0 0 0)');
>
> Should I open a GitHub issue?

Yep, but I'd like to undertand why "geom geometry(geometry, 0)" is not
the same as "geom geometry".


In the mean time you should be able to create your own Geometry type.
Something like this (untested):

import geoalchemy2.types

class Geometry(geoalchemy2.types.Geometry):
def get_col_spec(self):
if self.geometry_type == 'GEOMETRY' and self.srid == 0:
return self.name
return '%s(%s,%d)' % (self.name, self.geometry_type, self.srid)


--
Eric

[Tom Kralidis]

Thanks Eric.  This doesn't work either/yields the same results. I've issued a ticket in https://github.com/geoalchemy/geoalchemy2/issues/158
but if there is a workaround that would be great.

Thanks

..Tom

 

--
Eric

[Eric Lemoine]

On Thu, Jul 27, 2017 at 6:37 PM, Tom Kralidis <tom.kr...@gmail.com> wrote:
>
>
> On Thursday, July 27, 2017 at 11:49:58 AM UTC-4, erilem wrote:
>>
>> On Thu, Jul 27, 2017 at 2:05 PM, Tom Kralidis <tom.kr...@gmail.com> wrote:
>> > Thanks Tobias. It would be great to have the same functionality as
>> > achieved
>> > by:
>> >
>> > create table testptz (pkey integer primary key, geom geometry);
>> > insert into testptz values (0, 'POINT(0 0 0)');
>> >
>> > Should I open a GitHub issue?
>>
>> Yep, but I'd like to undertand why "geom geometry(geometry, 0)" is not
>> the same as "geom geometry".
>>
>>
>> In the mean time you should be able to create your own Geometry type.
>> Something like this (untested):
>>
>> import geoalchemy2.types
>>
>> class Geometry(geoalchemy2.types.Geometry):
>> def get_col_spec(self):
>> if self.geometry_type == 'GEOMETRY' and self.srid == 0:
>> return self.name
>> return '%s(%s,%d)' % (self.name, self.geometry_type, self.srid)
>>
>>
>
> Thanks Eric. This doesn't work either/yields the same results. I've issued
> a ticket in https://github.com/geoalchemy/geoalchemy2/issues/158
> but if there is a workaround that would be great.

I'll give it a try when I find time.


--
Eric

[Eric Lemoine]

This works for me:


-----
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy import Column, Integer, create_engine
from sqlalchemy.orm import sessionmaker
import geoalchemy2

class Geometry(geoalchemy2.types.Geometry):
def get_col_spec(self):
if self.geometry_type == 'GEOMETRY' and self.srid == 0:
return self.name
return '%s(%s,%d)' % (self.name, self.geometry_type, self.srid)

Base = declarative_base()


class Foo(Base):
__tablename__ = 'foo'
id = Column(Integer, primary_key=True)
geom = Column(Geometry(srid=0))


engine = create_engine('postgresql://elemoine@localhost/geoalchemy', echo=True)
Session = sessionmaker(bind=engine)
session = Session()


Foo.__table__.drop(engine, checkfirst=True)
Foo.__table__.create(engine)

foo = Foo()
foo.id = 1
foo.geom = 'POINT(0 0)'

session.add(foo)
session.commit()

foo = Foo()
foo.id = 2
foo.geom = 'POINTZ(0 0 0)'

session.add(foo)
session.commit()
------



As I said I'd really like to understand why "geom geometry" is more
permissive than "geom geometry(GEOMETRY, 0)". Have you asked the
question on the PostGIS mailing list or something? If not then I'd
suggest that we do it. Can you do that? Should I do it?

Thanks.




--
Eric