2

I have an irregular polygon, the perimeter of which is defined by an array of Cartesian coordinates.

I'm looking for a method to find the minimum and maximum y coordinate given an x coordinate that would work for a continuous range from x-min to x-max.

I would imagine one way to do this would be to define equations of the straight lines between each point and then apply the two of these whose range is satisfied by a given x coordinate.

Is there a quicker method to calculate/implement? Or is there a module that would allow me to do this easily?

I have been using matplotlib.path to plot the polygon - perhaps this class can help?

Thanks

BJH
  • 443
  • 1
  • 5
  • 11
  • By array of Cartesian coordinates - do you mean the polygon is defined by lines `array[i] -> array[i+1]`? – michaelrccurtis Aug 28 '15 at 16:52
  • Yes exactly, there are not even any curves, so the lines between points can be defined as equations trivially. – BJH Aug 28 '15 at 16:55
  • I think your suggestion is probably easiest - bear in mind that you can just do a simple interpolation, so it's pretty straightforward. – michaelrccurtis Aug 28 '15 at 17:00
  • 1
    See `shapely` https://pypi.python.org/pypi/Shapely which is a python package for doing geometric computation like this. – tacaswell Aug 28 '15 at 17:30

2 Answers2

3

Per @tcaswell's suggestion, shapely makes this fairly easy. In particular, the approach is to compute the .intersection (i.e. overlap) between the boundary of your polygon, and a vertical line at your chosen x-location.

A key thing to note here is that shapely polygons are filled objects, such that computing the overlap of these polygons with a line will return the segment(s) of the line that are inside the polygon. Therefore, to get the points where the vertical line crosses the edge of the polygon, it's probably simplest to compute the intersection between the .boundary attribute of a polygon and the vertical line.

Implementation

Here is a function returns the intersection points of an input polygon with a vertical line at a specified x value:

import shapely.geometry as sg

def polygon_intersect_x(poly, x_val):
    """
    Find the intersection points of a vertical line at
    x=`x_val` with the Polygon `poly`.
    """
    if x_val < poly.bounds[0] or x_val > poly.bounds[2]:
        raise ValueError('`x_val` is outside the limits of the Polygon.')
    if isinstance(poly, sg.Polygon):
        poly = poly.boundary
    vert_line = sg.LineString([[x_val, poly.bounds[1]],
                               [x_val, poly.bounds[3]]])
    pts = [pt.xy[1][0] for pt in poly.intersection(vert_line)]
    pts.sort()
    return pts

Example Usage

First create an example shapely.geometry.Polygon:

p = sg.Polygon([(0, 0),
                (10, 0),
                (30, 10),
                (70, -50),
                (80, 30),
                (40, 40),
                (60, 80),
                (50, 100),
                (15, 20),
                (0, 0)])

Note that you can also create shapely objects by: a) wrapping NumPy arrays to shapely types, b) Converting MPL paths to shapely polygons as p = sg.Polygon(my_path.vertices).

The function above can now be used to calculate min/max intersection points as:

x_val = 45
points = polygon_intersect_x(p, x_val)
minmax = [points[0], points[-1]]  # The values are already sorted
print minmax
# [-12.5, 88.5714286]

Here's a simple plot that demonstrates the approach:

import matplotlib.pyplot as plt

ax = plt.gca()

ax.fill(*p.boundary.xy, color='y')
ax.axvline(x_val, color='b')

ax.plot([x_val, ] * len(points), points, 'b.')
ax.plot([x_val, ] * 2, minmax, 'ro', ms=10)

Example Figure: Red dots are the 'minmax' values, blue dots are the other intersections.

farenorth
  • 10,165
  • 2
  • 39
  • 45
1

Here's a quick and dirty example of the interpolation:

def y_on_line(x, pos1, pos2):
    if(x < min(pos1[0], pos2[0])) or (x > max(pos1[0], pos2[0])):
        return None
    else:
        if(pos1[0] == pos2[0]):
            return pos1[1]
        else:
            return pos1[1] + (pos2[1] - pos1[1]) * (x - pos1[0]) / (pos2[0] - pos1[0])

def min_max_y(x, array):
    miny = None
    maxy = None
    for i in range(len(array)):
        y = y_on_line(x, array[i], array[(i+1)%len(array)])
        if(y is None):
            continue
        if(y < miny) or (miny is None):
            miny = y
        if(y > maxy) or (maxy is None):
            maxy = y
    return [miny, maxy]

print min_max_y(0.5, [[0,0],[0,2],[1,1]])

output: [0.5,1.5]

Note that if your polygon is concave then there's no reason that (x, miny) - (x,maxy) is all contained within it, in case that matters.

michaelrccurtis
  • 1,172
  • 3
  • 14
  • 16