How to extract edges from polydata as connected features?

I have a polydata structure and its extracted edges but computed with extract_feature_edges function as unconnected cells (separated lines).

Is it possible to connect those cells (lines) from their common points and then get the different features (lands, islands such as what you can see in the image - Antartica, Australia, ... - BTW they are paleo continents)?

In resume, I would like to extract from my grid and its edges the different land parts as separate polydata. I have tried with the python module shapely and the polygonize function, it works but not with 3D coordinates (https://shapely.readthedocs.io/en/latest/reference/shapely.polygonize.html).

import pyvista as pv

! wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/pyvista/mesh.vtk
mesh = pv.PolyData('mesh.vtk')
edges = mesh.extract_feature_edges(boundary_edges=True)

pl = pv.Plotter()

pl.add_mesh(pv.Sphere(radius=0.999, theta_resolution=360, phi_resolution=180))
pl.add_mesh(mesh, show_edges=True, edge_color="gray")
pl.add_mesh(edges, color="red", line_width=2)

viewer = pl.show(jupyter_backend='pythreejs', return_viewer=True)
display(viewer)

Any idea?

Solution

Here is a solution using vtk.vtkStripper() to join contiguous segments into polylines. See thread from https://discourse.vtk.org/t/get-a-continuous-line-from-a-polydata-structure/9864

import pyvista as pv
import vtk
import random

! wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/pyvista/mesh3D.vtk
mesh = pv.PolyData('mesh3D.vtk')
edges = mesh.extract_feature_edges(boundary_edges=True)

pl = pv.Plotter()

pl.add_mesh(pv.Sphere(radius=0.999, theta_resolution=360, phi_resolution=180))
pl.add_mesh(mesh, show_edges=True, edge_color="gray")

regions = edges.connectivity()
regCount = len(set(regions.get_array("RegionId")))

connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
stripper = vtk.vtkStripper()

for r in range(regCount):
    connectivityFilter.SetInputData(edges)
    connectivityFilter.SetExtractionModeToSpecifiedRegions()
    connectivityFilter.InitializeSpecifiedRegionList()
    connectivityFilter.AddSpecifiedRegion(r)
    connectivityFilter.Update()

    if r == 11:
        p = pv.PolyData(connectivityFilter.GetOutput())
        p.save('poly1.vtk')
    
    stripper.SetInputData(connectivityFilter.GetOutput())
    stripper.SetJoinContiguousSegments(True)
    stripper.Update()
    reg = stripper.GetOutput()
    
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    pl.add_mesh(reg, color=random_color, line_width=4)

viewer = pl.show(jupyter_backend='trame', return_viewer=True)
display(viewer)