I narrowed a 3D problem to a GeoDataFrame G containing four Shapely polygons, each is a 3D flat triangle. In 3 dimensions, these polygons all share a common summit and also a common edge, 2-by-2. They never overlap.

Here's how it looks like:

In [2]: G
0  POLYGON Z ((2848000. 1091400. 1500., 2865000. 1111500. 2000., ...
1  POLYGON Z ((2848000. 1091400. 1500., 2841500. 1121700.    0., ...
2  POLYGON Z ((2848000. 1091400. 1500., 2845000. 1109000.  100., ...
3  POLYGON Z ((2848000. 1091400. 1500., 2837000. 1101000.  700., ...

For visualization purposes, I want to represent these on a 2D map. As my coordinates are already in a projected CRS (2056), it's OK to simply drop the Z value. But then, when keeping only the (x,y) coordinates, these polygons sometime overlap:

overlapping 3D triangles

This is fine as shown by this QGIS screenshot.
Here, I loaded a GeoPackage file to which I exported my GeoDataFrame, each triangle being a single record in that file, hence 4 features.

But now I would like to stack them all together into a single MultiPolygon (either 2D or 3D), because, let's say, they are members of a single family, representing a 3D object.

If I do so using a shapely.geometry.MultiPolygon instance:

from shapely.geometry import MultiPolygon
M = MultiPolygon(list(G['geometry']))

it results in an invalid MultiPolygon (M.is_valid returns False). I don't know why because each of the original individual shape is a valid polygon. But I actually suspect that it's simply because Shapely doesn't really take the Z value into consideration. If it were a true MultiPolygonZ, I hope it would have been valid, but Shapely is not meant to deal with 3D data. And that's fine.
Please note that the resulting MultiPolygon is also not valid in the case where no triangle overlaps with the others.
Also, M.has_z always evaluates to True.

That's why I'm asking; do the overlaps cause the invalidity?
In such case, is it possible to build:

  • either a valid 2D MultiPolygon having overlapping features (because of the projection from 3D non-overlapping features to 2D) / coalesce each row of the GeoDataFrame into a single one with such 2D MultiPolygon?
  • or a valid 3D MultiPolygonZ / coalesce each row of the GeoDataFrame into a single one with such 3D MultiPolygonZ?

In most 2D cases these overlaps don't make sense, but in this particular one where these polygons are actual projections of a family of 3D features, this definitely makes sense.

In any case, if I simply cast the MultiPolygon instance back into a GeoDataFrame (gpd.GeoDataFrame(M)), it's automatically split in its parts. And in the doc I found no words on multipolygons: https://geopandas.org/docs/reference/api/geopandas.GeoDataFrame.html

The final idea would be to write the 2D MultiPolygon either to a GeoPackage having a single record or to a PostGIS table which has a MultiPolygon geometry field also as a single record. Same things for the 3D use case, in which case the PostGIS table has a MultiPolygonZ geometry field (and I don't know how does GeoPackage handle 3D objects.)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.