2

Context

In a Python script, I have a 3D square pyramid composed of 4 triangles (a list of 4 Shapely 3D polygons) that I need to store within a PostgreSQL database (+PostGIS).

Here's what this pyramid looks like when projected on a cartographic plane as a Shapely MultiPolygon:

enter image description here

I am able to store this MultiPolygon as such in PostGIS even so from a cartographic perspective, it is obviously not a valid MultiPolygon because of the massive overlaps between triangles. But still, it's a valid 3D pyramid.

Therefore, I was thinking of using the POLYHEDRALSURFACE geometry type, which is more suited to this use case. But then, this data type doesn't seem to be available in Shapely, hence the following error when trying to write the MultiPolygon object to a POLYHEDRALSURFACE field in the database:

InvalidParameterValue:
    Geometry type (MultiPolygon) does not match column type (PolyhedralSurface)

Here's the type of SQL query I'm running using psycopg2 for this purpose (here with dummy data):

import psycopg2
from shapely import wkt

myMultiPolygon = wkt.loads("""MULTIPOLYGON Z (
  ((0 0 0,1 1 1,1 2 1,0 0 0)),
  ((0 0 0,1 2 1,3 2 4,0 0 0)),
  ((0 0 0,3 2 4,2 2 1,0 0 0)),
  ((0 0 0,2 2 1,1 1 1,0 0 0)))"""
)

sql = f"""UPDATE mytable set poly3d = st_transform(
              st_geomfromtext('{mymultipoly}', 2056),
              4326)
          WHERE id=1;"""

conn = psycopg2.connect(**db_params)
curs = conn.cursor()
curs.execute(sql, vars={"mymultipoly": wkt.dumps(myMultiPolygon)})

Question

So, is there an elegant (by elegant I mean, without manually tweaking the WKT string...) solution to store those 4 triangles as a whole, namely as a single valid 3D object, in a PostgreSQL database as a "POLYHEDRALSURFACE Z"?

Just in case, I also have the four arrays containing each 3 vertices for each triangle.

EDIT:

What I'm actually trying to do:

Usually, one goes this way to convert a JSON (dict) to a WKT string if it's a GeoJSON like object:

from shapely import wkt
from shapely.geometry import mapping, shape

s = """MULTIPOLYGON Z (
((0 0 0, 1 1 1, 1 2 1, 0 0 0)),
((0 0 0, 1 2 1, 3 2 4, 0 0 0)),
((0 0 0, 3 2 4, 2 2 1, 0 0 0)),
((0 0 0, 2 2 1, 1 1 1, 0 0 0)))"""

d = mapping(wkt.loads(s))

# is equivalent to:
d = {'type': 'MULTIPOLYGON',
 'coordinates': [(((0.0, 0.0, 0.0),
    (1.0, 1.0, 1.0),
    (1.0, 2.0, 1.0),
    (0.0, 0.0, 0.0)),),
  (((0.0, 0.0, 0.0), (1.0, 2.0, 1.0), (3.0, 2.0, 4.0), (0.0, 0.0, 0.0)),),
  (((0.0, 0.0, 0.0), (3.0, 2.0, 4.0), (2.0, 2.0, 1.0), (0.0, 0.0, 0.0)),),
  (((0.0, 0.0, 0.0), (2.0, 2.0, 1.0), (1.0, 1.0, 1.0), (0.0, 0.0, 0.0)),)]}

wkt.dumps(shape(d))

will return the original WKT string (actually with a bunch more of 0 after the digit):

'MULTIPOLYGON Z (
((0.000 0.000 0.000, 1.000 1.000 1.000, 1.000 2.000 1.000, 0.000 0.000 0.000)),
((0.000 0.000 0.000, 1.000 2.000 1.000, 3.000 2.000 4.000, 0.000 0.000 0.000)),
((0.000 0.000 0.000, 3.000 2.000 4.000, 2.000 2.000 1.000, 0.000 0.000 0.000)),
((0.000 0.000 0.000, 2.000 2.000 1.000, 1.000 1.000 1.000, 0.000 0.000 0.000)))'

Easy.

But now, with a custom geometry type given in the type key of the dict, it won't work anymore:

d = {'type': 'POLYHEDRALSURFACE Z',
 'coordinates': [(((0.0, 0.0, 0.0),
    (1.0, 1.0, 1.0),
    (1.0, 2.0, 1.0),
    (0.0, 0.0, 0.0)),),
  (((0.0, 0.0, 0.0), (1.0, 2.0, 1.0), (3.0, 2.0, 4.0), (0.0, 0.0, 0.0)),),
  (((0.0, 0.0, 0.0), (3.0, 2.0, 4.0), (2.0, 2.0, 1.0), (0.0, 0.0, 0.0)),),
  (((0.0, 0.0, 0.0), (2.0, 2.0, 1.0), (1.0, 1.0, 1.0), (0.0, 0.0, 0.0)),)]}

wkt.dumps(shape(d))

error: GeometryTypeError: Unknown geometry type: polyhedralsurface z

So the problem could be reduced to: dumping this dict to a WKT string, but by "skipping" the geometry type validation in order to get a valid WKT string out of it, whatever the geometry type is:

'POLYHEDRALSURFACE Z (
((0.000 0.000 0.000, 1.000 1.000 1.000, 1.000 2.000 1.000, 0.000 0.000 0.000)),
((0.000 0.000 0.000, 1.000 2.000 1.000, 3.000 2.000 4.000, 0.000 0.000 0.000)),
((0.000 0.000 0.000, 3.000 2.000 4.000, 2.000 2.000 1.000, 0.000 0.000 0.000)),
((0.000 0.000 0.000, 2.000 2.000 1.000, 1.000 1.000 1.000, 0.000 0.000 0.000)))'

I wish I could make shapely "blind" about the geometry type.

Version info:
– postgresql: 14.1 (Debian 14.1-1.pgdg110+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 10.2.1-6) 10.2.1 20210110, 64-bit
– postgis_version(): 3.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
– postgis_sfcgal_version(): 1.3.8

1
  • For the moment I'm using this workaround: myPolyhedral = wkt.dumps(myMultiPolygon).replace('MULTIPOLYGON', 'POLYHEDRALSURFACE') and then insert it using: ST_GeomFromText(<geom>, <SRS>) but I'm not happy with that solution. Commented Sep 22, 2022 at 20:27

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.