Shapely: punto de intersección entre línea y polígono en 3D

La última vez que utilicé bien, tuve este bonito sentimiento de importación y vuelo . Sin embargo, recientemente, me encuentro con un comportamiento bastante no intuitivo en este módulo, ya que traté de encontrar la intersección entre un segmento de línea y un triángulo en el espacio 3D. Vamos a definir un segmento y un triángulo como sigue:

l = LineString([[1,0.5,0.5],[3,0.5,0.5]]) p = Polygon([[1.2,0.0,0.],[2.2,1.0,0.],[2.8,0.5,1.]]) 

Para obtener su punto de intersección utilicé l.intersection(p) , y esperé un punto, a saber POINT Z (POINT Z (2 0.5 0.25)) . Está ilustrado con el punto azul abajo:

introduzca la descripción de la imagen aquí

En su lugar, lo que obtuve fue LINESTRING Z (1.7 0.5 0.25, 2.8 0.5 1) – línea roja debajo – y francamente estoy bastante perplejo acerca de lo que se supone que representa. introduzca la descripción de la imagen aquí

Por extraño que parezca, cuando el polígono / triángulo está en el plano xz y ortogonal al segmento de línea, la función se comporta como uno esperaría. Cuando el triángulo está “inclinado”, sin embargo, devuelve una línea. Esto me ha llevado temporalmente a creer que devolvió la intersección entre la línea y el cuadro delimitador del triángulo. La línea roja anterior demuestra lo contrario.

Así que una solución para este problema ha sido leer esta página web muy esclarecedora , y adaptar su código C++ para trabajar con objetos bien formados. El método de intersection funciona a la perfección para verificar si la línea atraviesa el polígono y la siguiente función encuentra el punto de interés.

 def intersect3D_SegmentPlane(Segment, Plane): # Points in Segment: Pn Points in Plane: Qn P0, P1 = np.array(Segment.coords) Q0, Q1, Q2 = np.array(Plane.exterior)[:-1] # vectors in Plane q1 = Q1 - Q0 q2 = Q2 - Q0 # vector normal to Plane n = np.cross(q1, q2)/np.linalg.norm(np.cross(q1, q2)) u = P1 - P0 # Segment's direction vector w = P0 - Q0 # vector from plane ref point to segment ref point ## Tests parallelism if np.dot(n, u) == 0: print "Segment and plane are parallel" print "Either Segment is entirely in Plane or they never intersect." return None ## if intersection is a point else: ## Si is the scalar where P(Si) = P0 + Si*u lies in Plane Si = np.dot(-n, w) / np.dot(n, u) PSi = P0 + Si * u return PSi 

Ya no importa mucho y vuela …

Así que finalmente a mis preguntas:

  • ¿Qué devuelve la intersection cuando se aplica a objetos 3D y por qué es una línea?

  • ¿Hay una función en forma que hace lo que quiero? ¿O cualquier argumento opcional, truco o magia oscura?

  • ¿Hay alguna otra biblioteca por ahí que haga este trabajo mientras cumple mis sueños de simplicidad y pereza?

Desafortunadamente, como indica la documentación:

Las secuencias de coordenadas son inmutables. Se puede usar un tercer valor de coordenada z al construir instancias, pero no tiene efecto en el análisis geométrico. Todas las operaciones se realizan en el plano xy.

Uno puede verificar esto con:

 from shapely.geometry import LineString, Polygon l = LineString([[1,0.5,0.5],[3,0.5,0.5]]) p = Polygon([[1.2,0.0,0.],[2.2,1.0,0.],[2.8,0.5,1.]]) print(l.intersection(p)) #LINESTRING Z (1.7 0.5 0.25, 2.8 0.5 1) l = LineString([[1,0.5],[3,0.5]]) p = Polygon([[1.2,0.0],[2.2,1.0],[2.8,0.5]]) print(l.intersection(p)) #LINESTRING (1.7 0.5, 2.8 0.5) 

o incluso:

 from shapely.geometry import LineString, Polygon l = LineString([[1,0.5,0],[3,0.5,0]]) p = Polygon([[1.2,0.0,1],[2.2,1.0,1],[2.8,0.5,1]]) print(l.intersects(p)) #True (even though the objects are in different z-planes)