How to perform union on polygons in pygplates

Hello everyone,

May I ask for suggestion on how to perform union on polygon features in pygplates? I understand that pygplates provides wonderful functions for PolygonOnSphere such as pygplates.PolygonOnSphere.partition to find a relationship between a PolygonOnSphere and any other GeometryOnSphere. I want to perform a GIS union on polygons that are intersecting with each other and get the coordinates of the union polygon. I have tried to apply different modules such as pyproj and shapely in Python scientific package to perform union. However, I ended up with different artifacts because some geometries crossing the dateline. From what I understand, most GIS modules do not recognize that -180 and +180 longitude is the same line which causes unusual behavior.

Many thanks for your time and your help.

Hi Lavie,

The pygplates.PolygonOnSphere.partition function is really just meant for intersections (as opposed to unions). One way to use it to find the union of two polygons is to create a uniform, dense multi-point (MultipointOnSphere) in the area covered by the two polygons and find which points are inside either (or both) polygons (using pyGPlates to partition the multi-point). And then use some kind of outline algorithm to find the boundary of those points and create a polygon (possibly with interior holes) from that. However this approach is a bit cumbersome and has the problem that the final union polygon would look a bit gridded (depending on the density of the initial points).

Alternatively, what you want may be possible with the Boost geometry library - I’m not sure - they mention a union operation. There’s mention of supporting geographical coordinate systems (ie, points on a sphere) but they mention that that is an extension where they say:

Note that the library extensions are not distributed in the official Boost releases, but only available in the Boost.Geometry (develop branch) and that they are subject to change.

I would assume that the result you should get from these packages would be equivalent to the final union polygon (calculated without any dateline clipping) followed by dateline clipping that final intersection polygon. Because presumably they would start with the two initial polygons (but that are each implicitly dateline clipped because they are in latitude/longitude coordinates) and then perform the union of those which I think should give the same result. But maybe that’s not happening.

On a side note: Currently pygplates.PolygonOnSphere.partition does not properly partition a polygon (in other words when you intersect a polygon with another polygon you only get the clipped polyline pieces rather than the actual intersected polyon - there’s a warning in the docs here that goes into more detail). In a future release soon, proper polygon-polygon intersection (with interior holes) will be implemented, actually it’s currently implemented in the dateline wrapper but that’s a simpler case. Also intersection is just one of the outputs of doing this properly, other boolean operations like union and difference can be made available with minor modifications.

Regards,
John

Thank you so much for your advice. Your advice has guided me to solve the problem :slight_smile: