5

Last time I used shapely, I really had this nice import-and-fly-feeling. However recently, I run into a rather non-intuitive behavior in this module, as I tried to find the intersection between a line segment and a triangle in 3D space. Let's define a segment and a triangle as follow:

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.]])

To get their point of intersection I used l.intersection(p), and expected a point, namely POINT Z (POINT Z (2 0.5 0.25)). It is illustrated with the blue dot below:

enter image description here

Instead what I got was LINESTRING Z (1.7 0.5 0.25, 2.8 0.5 1) - red line below - and frankly I am quite perplexed about what it is supposed to represent. enter image description here

Oddly enough, when the polygon/triangle is in the xz-plane and orthogonal to the line segment, the function behaves as one would expect. When the triangle is "leaning", however, it returns a line. This has temporarily lead me to believe that it returned the intersection between the line and the triangle's bounding box. The above red line proves otherwise.

So a workaround for this problem has been to read this very enlightening web-page, and adapt their C++ code to work with shapely objects. The intersection method works great to check whether the line goes through the polygon, and the function below finds the point of interest. 

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

Not very import-and-fly anymore...

So finally to my questions:

  • What does intersection return when applied to 3D-objects and why is it a line?

  • Is there a function in shapely that does what I want? or any optional argument, tweak or dark magic trick?

  • Is there any other libraries out there that would do this job while fulfilling my dreams of simplicity and laziness?

Georgy
  • 12,464
  • 7
  • 65
  • 73
snake_charmer
  • 2,845
  • 4
  • 26
  • 39

1 Answers1

7

unfortunately, as the documentation states:

Coordinate sequences are immutable. A third z coordinate value may be used when constructing instances, but has no effect on geometric analysis. All operations are performed in the x-y plane.

One can verify this with:

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)

or even:

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)
ewcz
  • 12,819
  • 1
  • 25
  • 47