Ok. Less trivial as you wrote.
No vtk higher level function to do that ?
the orientation? i’m not sure… maybe a combination of vtkReverseSense and vtkPolyDataNormals with ConsistencyOn()…
This kind of specific settings (here made for only 1 cross) do not resist when you want to generalize the algorithm.
Hello,
Maybe you can try the Shapely library: python - Splitting self-intersecting polygon only returned one polygon in Shapely - Stack Overflow
best,
PC
Hello,
@marcomusy’s answer uses NumPy for performance critical operations (e.g. singular value decomposition), if that is your concern. And certainly the for loop can be avoided by vectorized calls to the cross product and the vector norm. NumPy is supposed to be a decent alternative to its Matlab equivalents. Vectorizing NumPy functions normally takes the standard of telling it how to iterate over the vectors according to how they are arranged in a matrix (e.g. np.linalg.norm(x, axis=1) if you have the point data as row-vectors). See this: python - How to apply numpy.linalg.norm to each row of a matrix? - Stack Overflow . The cross product function likely follows the same rationale.
regards,
PC
Yes indeed. Have tried but shapely does not work in 3D.
I wanted to have separate polylines from the vtk process.
So my concern is why there is 2 cells detected for the region 7.
I do not understand why.
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=6300E3, 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()
#r_choice = 7
for r in range(regCount):
#for r in [r_choice]:
connectivityFilter.SetInputData(edges)
connectivityFilter.SetExtractionModeToSpecifiedRegions()
connectivityFilter.InitializeSpecifiedRegionList()
connectivityFilter.AddSpecifiedRegion(r)
connectivityFilter.Update()
stripper.SetInputData(connectivityFilter.GetOutput())
stripper.SetJoinContiguousSegments(True)
stripper.Update()
cleaner = vtk.vtkCleanPolyData()
cleaner.SetInputData(stripper.GetOutput())
cleaner.Update()
reg = cleaner.GetOutput()
print(r, reg.GetNumberOfCells(), reg.GetNumberOfPoints())
if r == r_choice:
p = pv.PolyData(reg)
p.save('poly_%02d.vtk' %(r_choice))
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)
0 1 294
1 1 556
2 1 29
3 1 8
4 1 80
5 1 838
6 1 24
7 2 56 <----- 2 cells !
8 1 38
9 1 42
10 1 22
11 1 19
12 1 22
13 1 44
14 1 44
15 1 26
16 1 90
17 1 12
18 1 11
19 1 14
20 1 19
21 1 11
Your problem is not truly 3-D. Your data are locations on the Earth, so working with their XYZ positions is an unnecessary complication. I strongly recommend working with lat-long polar values or use some projection suitable for world maps (e.g. Mercator) if you want to work with Cartesian coordinates.
Not agree because finding convex polylines is easier with 3D cartesian coordinates.
Find for example the polyline that describes Antartica continent.
You can convert XYZ absolute postions to polar coordinates with this example: python - Faster numpy cartesian to spherical coordinate conversion? - Stack Overflow . The first solution is a vectorized (loop-free) example. Just throw away the radius (will be constant for points in the surface of a sphere) and you’ll reduce your problem to 2D.
From my experience I can tell you: that’s the wrong way to process geographic data. But if you want to stick to 3D data, be warned.
Always hard to have links to other approaches when I have tried to describe the problem with a concrete code ready to be tested.
See the problem described here bellow that comes from polygons that are folded.
My strategy was to find each individual convex polylines from vtk, then indeed use shapely to set the order. Then produces maps with cartopy.
They are not “other approaches”. They are steps to solve problems. I think you want a magic bullet to solve a complex problem in a single shot, but in reality this is rarely the case. Complex data processing must be broken down into small manageable problems. That’s what me and others have been trying to do. Then it’s your mission to put them together to engineer a workflow to suit your need. You should not expect others to do all the job for you (at least for free). No wonder that SO question got downvoted (that was not me, I swear).
I have really try to make my question accessible by separating it into smaller problems.
Giving also the big picture beside.
So my questions to VTK Discourse were clear enough to be analysed by more experimented developpers in my humble opinion.
Why I am not able to produce separate cells from the vtkStripper ? Why the region 7 produces 2 cells ?
Sorry No, I do not expect other do my work. That is really not correct to answer this.
I thank @marcomusy that indeed helps me to discover vedo library and really answered from the code I had prepared. Other answers were very far from what I was expecting indeed.
Why vtkStripper.SetJoinContiguousSegments() seems to have no differences when set to True or False ?
This is the kind of answers I would like to have…
Sorry, I mean no offense. But you reached SO with a “here’s my code, it has a bug, please, find it for me”. That’s not a question, that’s a challenge. The community frowns upon stuff like that, unless you use SO’s bounty-hunting mechanism. Place a bounty and maybe someone will accept the challenge… Just an advice.
Now back to the problem in hand:
Where’s the figure showing the region 7?
Really not easy to be understood. With a lot of caricatural behaviour you intend to put on me.
Again, I do not expect that someone solves my all problem but perhaps yes to have a look into some points I would have missed in the snippets of code each time prepared to illustrate the problem.
The snippet that produces 2 cells for region 7 is above.
The reply where I have put the output.
0 1 294
1 1 556
2 1 29
3 1 8
4 1 80
5 1 838
6 1 24
7 2 56 <----- 2 cells !
8 1 38
...
I hope these figures helps illustrating the importance of reducing your problem to 2D.
This is a simple example of a seemingly self intersecting polygon whose vertexes belongs to a sphere:
However, the fact that the vertexes are on the surface doesn’t mean the edges are:
As you can see, the edges may not cross in 3D. That likely explains the polygon-splitting algorithm fails here and there. Why don’t you work in 2D space (e.g. lat-long or projected XY) and try the algorithms again?
Surely Antarctica results in an open poly when projected in Mercator projection, for example. Use your creativity. It’s not difficult to write a routine that closes the Antarctica poly:
The red line shows how to add two more vertexes and close it.
Furthermore, you can use another projection. The azimuthal equidistant projection, for example, will result in all-closed polygons, even Antarctica.
Thank you again for the your investigation but really do no understand your suggestions.
Some replies however.
-
“Use your creativity.” is really not correct to write. Is it a discourse where to learn how to be ?
-
“It’s not difficult to write a routine that closes the Antarctica poly”, I don’t think you get the point between doing a process one time and write a general algorithm. I have configuration with paleo-continents where they are placed differently. I don’t want to do custom settings and solve specific problems for a particular polygon !
The ccrs.geodetic() projection handles correclty cross-date problems by the way. So do not need to add extra polygons to fill it. Please have a look to python - How to detect inner polygons from a multipolygon shapely object? - Stack Overflow where you will see polygons correctly filled without adding any extra points ! Just a matter of describing the polygon in the right order and with continuous lines.
And again, I would prefer to understand the vtk.vtkStripper() problem I have submitted rather than other things.
But thanks again.
Have a look to this code if you want. That was the start of my investigation.
And my choice to separate each polygons since they seems not ordered correctly when folded.
And I still think that this step can be achieve from the VTK processing with 3D cartesian coordinates.
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random
import pickle
#-----------------------------------------
! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle
with open('./polys.pickle', "rb") as poly_file:
polygons = pickle.load(poly_file)
#-----------------------------------------
def polygon2patch(poly, **kwargs):
path = Path.make_compound_path(
Path(poly.exterior.coords, closed=True),
*[Path(ring.coords) for ring in poly.interiors])
patch = PathPatch(path, **kwargs)
return patch
#-----------------------------------------
fig = plt.figure(figsize=(10,5))
#map_proj = ccrs.PlateCarree(-10)
#map_proj = ccrs.Orthographic(-10, 80) # fill not correct
map_proj = ccrs.Orthographic(-10, 60)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)
transform = ccrs.Geodetic()
holesNumber = []
for n,polygonA in enumerate(polygons.geoms):
# 2, 11
#if n in [2,11]: continue
if n != 5: continue
holes = []
for m,polygonB in enumerate(polygons.geoms):
if (m == n): continue
if polygonA.contains(polygonB):
holes.append(polygonB.exterior.coords)
holesNumber.append(m)
if n in holesNumber: continue # n is a hole
polygonNew = geometry.Polygon(polygonA.exterior.coords, holes=holes)
polygonNew = shapely.geometry.polygon.orient(polygonNew) # Orient to be oriented counter-clockwise
random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
patch = polygon2patch(polygonNew, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
ax.add_patch(patch)
print(n, len(polygonNew.exterior.coords))
ax.set_global()
ax.gridlines()
plt.show()
With map_proj = ccrs.Orthographic(-10, 60)
With map_proj = ccrs.Orthographic(-10, 80)
The previous example works with another input files where
all self intersected polygons have been separated.
wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys_03.pickle
See details from answer given in
https://discourse.vtk.org/t/create-separated-regions-from-polydata
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random
import pickle
#-----------------------------------------
! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys_03.pickle
with open('./polys_03.pickle', "rb") as poly_file:
polygons = pickle.load(poly_file)
#-----------------------------------------
def polygon2patch(poly, **kwargs):
path = Path.make_compound_path(
Path(poly.exterior.coords, closed=True),
*[Path(ring.coords) for ring in poly.interiors])
patch = PathPatch(path, **kwargs)
return patch
#-----------------------------------------
fig = plt.figure(figsize=(10,5))
map_proj = ccrs.Orthographic(-10, 80)
#map_proj = ccrs.Orthographic(-10, 60)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)
transform = ccrs.Geodetic()
holesNumber = []
for n,polygonA in enumerate(polygons):
if n != 16: continue
holes = []
for m,polygonB in enumerate(polygons):
if (m == n): continue
if polygonA.contains(polygonB):
holes.append(polygonB.exterior.coords)
holesNumber.append(m)
if n in holesNumber: continue # n is a hole
polygonNew = geometry.Polygon(polygonA.exterior.coords, holes=holes)
polygonNew = shapely.geometry.polygon.orient(polygonNew) # Orient to be oriented counter-clockwise
random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
patch = polygon2patch(polygonNew, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
ax.add_patch(patch)
ax.set_global()
ax.gridlines()
plt.show()
