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:
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
myPolyhedral = wkt.dumps(myMultiPolygon).replace('MULTIPOLYGON', 'POLYHEDRALSURFACE')
and then insert it using:ST_GeomFromText(<geom>, <SRS>)
but I'm not happy with that solution.