2014-03-08

Classify Polygon Shape with Python

(FME 2014 build 14235)

From this thread > Community Answers: Parcels geometric shape

A Python script example for the task.
I think the getVertices function in this script can be used for general purpose.
-----
# PythonCaller Script Example: Classify Polygon Shape
# 2013-03-09 Added triangle classification.
# Classify a polygon into triangle, quadrilateral, and other.
# Triangle will be further classified into right, isosceles, and equilateral.
# Quadrilateral will be further classified into trapezoid, parallelogram, rectangle, and square.
import fmeobjects

def classifyPolygon(feature):
    shape = 'Invalid'
    g = feature.getGeometry()
    if isinstance(g, fmeobjects.FMEPolygon) and g.isBoundaryLinear():
        vertices = getVertices(feature)
        if len(vertices) == 3:
            shape = ','.join(classifyTriangle(vertices))
        elif len(vertices) == 4:
            shape = ','.join(classifyQuadrilateral(vertices))
        elif 4 < len(vertices):
            shape = 'Other'
    feature.setAttribute('_shape', shape)

X, Y, PRECISION = 0, 1, 1.0e-6
sqrLength = lambda a: a[X]**2 + a[Y]**2
isNearZero = lambda v: abs(v) < PRECISION
isZeroLength = lambda a: isNearZero(sqrLength(a))
isSameLength = lambda a, b: isNearZero(sqrLength(a) - sqrLength(b))
isParallel = lambda a, b: isNearZero(a[X] * b[Y] - a[Y] * b[X]) # outer product
isPerpendicular = lambda a, b: isNearZero(a[X] * b[X] + a[Y] * b[Y]) # inner product

# Remove excess coordinates and return tuple (x, y) list of polygon vertices.
def getVertices(feature):
    coords = feature.getAllCoordinates()
    if feature.getDimension() == fmeobjects.FME_THREE_D:
        coords = [(x, y) for x, y, z in coords]
    if len(coords) < 3:
        return coords
    x0, y0 = coords[0][X], coords[0][Y]
    x1, y1 = coords[1][X], coords[1][Y]
    vertices = [(x0, y0)]
    for x2, y2 in coords[2:]:
        a = (x1 - x0, y1 - y0)
        b = (x2 - x1, y2 - y1)
        if not isZeroLength(a) and not isZeroLength(b) and not isParallel(a, b):
            vertices.append((x1, y1))
            x0, y0 = x1, y1
        x1, y1 = x2, y2
    return vertices

# Return string list of triangle class names.
def classifyTriangle(vertices):
    # Create verctors of three edges.
    p1, p2, p3 = vertices[0], vertices[1], vertices[2]
    v1 = (p2[X] - p1[X], p2[Y] - p1[Y])
    v2 = (p3[X] - p2[X], p3[Y] - p2[Y])
    v3 = (p1[X] - p3[X], p1[Y] - p3[Y])
 
    # Classify the triangle.
    shape = ['Triangle']
    if isPerpendicular(v1, v2) or isPerpendicular(v2, v3):
        shape.append('Right')
    if isSameLength(v1, v2):
        shape.append('Isosceles')
        if isSameLength(v2, v3):
            shape.append('Equilateral')
    elif isSameLength(v2, v3):
        shape.append('Isosceles')
    return shape

# Return string list of quadrilateral class names.
def classifyQuadrilateral(vertices):
    # Create vectors of four edges.
    p1, p2, p3, p4 = vertices[0], vertices[1], vertices[2], vertices[3]
    v1 = (p2[X] - p1[X], p2[Y] - p1[Y])
    v2 = (p3[X] - p2[X], p3[Y] - p2[Y])
    v3 = (p4[X] - p3[X], p4[Y] - p3[Y])
    v4 = (p1[X] - p4[X], p1[Y] - p4[Y])
 
    # Classify the quadrilateral.
    shape = ['Quadrilateral']
    if isParallel(v1, v3):
        shape.append('Trapezoid')
        if isParallel(v2, v4):
            shape.append('Parallelogram')
            if isPerpendicular(v1, v2):
                shape.append('Rectangle')
                if isSameLength(v1, v2):
                    shape.append('Square')
    elif isParallel(v2, v4):
        shape.append('Trapezoid')
    return shape
-----

No comments:

Post a Comment