Python: How to define a new cross section from an arbitrary closed shape cut by a horizontal line

Question

I am trying to write some code that takes a user specified closed polygon, and then a user specified horizontal line to cut the polygon with the horizontal line to create a new shape. I also need to determine the area of the shape, and its perimeter. This is what I have so far:

I do not know how to create the new shape as I need to insert new points to form a new polygon? PS.. Have tried to insert this as a code block, but for some reason it re-formats it all the time ??? I have found bits of code to do some bits on Stack Overflow... Any assistance much appreciated...

enter code here
"""
ROUTINE TO FIND MULTIPLE INTERSECTION POINTS FOR A HORIZONTAL LINE CUTTING THROUGH AN  ARBITRARY CROSS SECTIONAL SHAPE

This routine should start with a poly defined by the user, 
take the Y_specified, to adjust the polygon shape, (Effectively cut the top off)

and then calculate the area of the polygon

Since Y_Specified may be between points of the defined polygon, the intersection points need to be found

"""
import pylab as pl
import numpy as np

def line(p1, p2):
    # This is how a line is defined
    A = (p1[1] - p2[1])
    B = (p2[0] - p1[0])
    C = (p1[0]*p2[1] - p2[0]*p1[1])
    return A, B, -C
def intersection(L1, L2):
    # This finds the intersection of 2 lines
    # Should I limit the Range here ???
    D  = L1[0] * L2[1] - L1[1] * L2[0]
    Dx = L1[2] * L2[1] - L1[1] * L2[2]
    Dy = L1[0] * L2[2] - L1[2] * L2[0]
    if D != 0:
        x = Dx / D
        y = Dy / D
        return x,y
    else:
        return False

# This defined a closed Hexagonal sort of shape... 30 wide and 15 high
poly = [[0,10],[10,0],[20,0],[30,10],[18,15],[12,15],[0,10]]
# This defines a W shape with a lid !!
poly = [[0,10],[10,0],[12,0],[15,13],[18,0],[20,0],[30,10],[18,15],[12,15],[0,10]]

x = tuple(x[0] for x in poly)
y = tuple(x[1] for x in poly)
print x
print y

pl.plot( x,y, 'go-')
pl.suptitle('CROSS SECTIONAL SHAPE', fontsize=14, fontweight='bold')    
pl.title('To be cut and replot it')
pl.xlabel('X-Co-ord')
pl.ylabel('Y-Co-ord')
pl.show()

print len(poly)
# This is the cut line
Y_Specified = 12.0
Orig_ymax = max(poly, key=lambda x: x[1])
Orig_ymax = Orig_ymax[1]
Orig_xmin = min(poly, key=lambda x: x[0])
Orig_xmin = Orig_xmin[0]
Orig_xmax = max(poly, key=lambda x: x[0])
Orig_xmax = Orig_xmax[0]
Orig_max_width = Orig_xmax - Orig_xmin

# DEFINE THE HORIZONTAL LINE
L2_PT1 = [Orig_xmin,Y_Specified]
L2_PT2 = [Orig_xmax,Y_Specified]
L2 = line(L2_PT1, L2_PT2 )
x2=[]
y2=[]
for n,i in enumerate(poly):
     if n+1< len(poly):
         L1_PT1 = poly[n]
         L1_PT2 = poly[n+1]
         L1 = line(L1_PT1, L1_PT2)
         R = intersection(L1, L2) # This correctly identifies the intersection points   within the width... but also outside ???
         if R:
             print "Intersection detected:", R
             print R[0],R[1]
             x2.append(R[0])
             y2.append(R[1])
         else:
             print "No single intersection point detected"
     print x2
     print y2
     # Now need to remove points from the original poly above the Y_Specified and insert the new intersection points

# That is if Poly[1] > Y_Specified (Delete it) and insert new point to create the newly formed smaller hexagonal shape ???

pl.plot( x,y, 'go-')
pl.plot(x2,y2, 'ro-')
# The new polygon needs to be defined by the green line cut off by the red line...???
# HOW CAN THIS BE DONE ???

pl.show()

Answer 1

You had most of it worked out already. There are two significant things I noticed.

Lines vs. Segments

I believe you got the line + intersection system from here? If so, then my understanding is that it is modeling lines instead of segments as you need which would result in the phantom intersections you are getting.

You will need to fix this with a different or updated intersection method before you can trust your clipped polygons.

Construction of New Polygon

I guess there is a more automatic / built-in way to do this but I couldn't find one.

Here is what I did: Rather than finding all the intersection points at one time this builds a new polygon by going around the vertexes and:

• include the vertex if it qualifies (for you, being over the line)
• also include any intersection between the vertex and the next vertex

If you can get an alternate line/intersection implementation, then you should be fine. Here is the result:

Copy-Paste Code

import pylab as pl
# important note - this system is X-major (x-width first, y-height second)

# line and intersection system (I didn't touch this part)
# looks like it came from here: https://stackoverflow.com/a/20679579/377366
def line(p1, p2):
    # This is how a line is defined
    A = (p1[1] - p2[1])
    B = (p2[0] - p1[0])
    C = (p1[0]*p2[1] - p2[0]*p1[1])
    return A, B, -C


def intersection(L1, L2):
    # This finds the intersection of 2 lines
    # Should I limit the Range here ???
    D  = L1[0] * L2[1] - L1[1] * L2[0]
    Dx = L1[2] * L2[1] - L1[1] * L2[2]
    Dy = L1[0] * L2[2] - L1[2] * L2[0]
    if D != 0:
        x = Dx / D
        y = Dy / D
        return x,y
    else:
        return False

# a polygon (cleaned up a little)
poly_original = [[0,10], [10,0], [12,0], [15,13], [18,0], [20,0], [30,10],
                 [18,15],[12,15],[0,10]]
x, y = zip(*poly_original)

# This is the cut line (I didn't touch this part)
Y_Specified = 12.0
Orig_ymax = max(poly_original, key=lambda x: x[1])
Orig_ymax = Orig_ymax[1]
Orig_xmin = min(poly_original, key=lambda x: x[0])
Orig_xmin = Orig_xmin[0]
Orig_xmax = max(poly_original, key=lambda x: x[0])
Orig_xmax = Orig_xmax[0]
Orig_max_width = Orig_xmax - Orig_xmin

L2_PT1 = [Orig_xmin, Y_Specified]
L2_PT2 = [Orig_xmax, Y_Specified]
L2 = line(L2_PT1, L2_PT2 )


# from here is a different algorithm than you used
# rather than finding all the intersection points at one time
# this builds a new polygon by going around the vertexes
# include each vertex if it qualifies (for you, being over the line)
# also include any intersection between the vertex and the next vertex
def is_qualified_vertex(v):
    # this could be changed depending on requirements
    # in your case, it is simply whether or not the y value is above (or equal?)
    # the horizontal line
    return v[1] >= Y_Specified


# I don't recommend including the same point twice to close a polygon since
# you are really defining vertices instead of edges, but this works with your
# original setup of duplicating the first/last point.
poly_new = list()
for i in range(0, len(poly_original)-1):
    vertex_a, vertex_b = poly_original[i:i+2]
    # decide to include the original vertex in the new poly or not
    # for your requirements, this just needs to be on or above a horizontal line
    if is_qualified_vertex(vertex_a):
        poly_new.append(vertex_a)
    # before moving on to the next vertex, check if there was an intersection
    # with the divider to be included in the new poly
    edge_line = line(vertex_a, vertex_b)
    # problem - I think you are using a line intersection algorithm
    #           rather than a segment intersection algorithm.
    #           That would be why you are getting phantom intersections.
    #           Also how do you want to handle overlapping segments?
    intersection_point = intersection(edge_line, L2)
    if intersection_point:
        # intersection - build new poly with it
        print('------------------')
        print(vertex_a)
        print(intersection_point)
        print(vertex_b)
        poly_new.append(intersection_point)
x2, y2 = zip(*poly_new)

pl.suptitle('CROSS SECTIONAL SHAPE', fontsize=14, fontweight='bold')
pl.title('Original and cut.')
pl.xlabel('X-Co-ord')
pl.ylabel('Y-Co-ord')
pl.plot(x, y, 'go-')
pl.plot(x2, y2, 'ro-')
pl.show()