Dan_Patterson

Errors in arcpy's Polygon class

Discussion created by Dan_Patterson Champion on Aug 5, 2010
Latest reply on Nov 8, 2012 by michalis
This code using arcpy's polygon class to calculate simple properties for simple shapes are not accurate.
'''
GeometryErrorsDemo.py

Demonstrates calculating various properties for polygons

'''
import arcpy
pnt = arcpy.Point()

triangle = [[0,0],[0,1],[1,0]]
square = [[0,0],[0,1],[1,1],[1,0]]
rectangle = [[0,0],[0,1],[2,1],[2,0]]
polygons = [triangle, square, rectangle]
labels = ["Triangle", "Square", "Rectangle"]
array = arcpy.Array()
polys = []
for i in range(len(polygons)):
  a_poly = polygons[i]
  print "\n", labels[i]," Coordinates: ", a_poly
  for pair in a_poly:
    pnt.X = pair[0]
    pnt.Y = pair[1]
    array.add(pnt)
    print " X %20.16f  Y %20.16f" % (pnt.X, pnt.Y)
  array.add(array.getObject(0))
  poly = arcpy.Polygon(array)
  print "Polygon properties:"
  print " area %20.16f\n length %s" % (poly.area, poly.length)
  print " centroid %s\n true centroid %s " % (poly.centroid, poly.trueCentroid) 
  print " first point %s\n last point %s " % (poly.firstPoint, poly.lastPoint)
  print " extent %s " % (poly.extent)
  print " hull %s\n " % (poly.hullRectangle)
  array.removeAll()
  polys.append(poly)


Results for a simple triangle, square and rectangle are as follows:

Triangle  Coordinates:  [[0, 0], [0, 1], [1, 0]]
X   0.0000000000000000  Y   0.0000000000000000
X   0.0000000000000000  Y   1.0000000000000000
X   1.0000000000000000  Y   0.0000000000000000
Polygon properties:
area   0.5001220777630806
length 3.41463033649
centroid 0.3333740234375 0.3333740234375 NaN NaN
true centroid 0.3333740234375 0.3333740234375 NaN NaN
first point 0 0 NaN NaN
last point 0 0 NaN NaN
extent 0 0 1.0001220703125 1.0001220703125 NaN NaN NaN NaN
hull 1.0001220703125 2.22044604925031E-16 0.50006103515625 -0.50006103515625 -0.50006103515625 0.50006103515625 0 1.0001220703125


Square  Coordinates:  [[0, 0], [0, 1], [1, 1], [1, 0]]
X   0.0000000000000000  Y   0.0000000000000000
X   0.0000000000000000  Y   1.0000000000000000
X   1.0000000000000000  Y   1.0000000000000000
X   1.0000000000000000  Y   0.0000000000000000
Polygon properties:
area   1.0002441555261612
length 4.00048828125
centroid 0.50006103515625 0.50006103515625 NaN NaN
true centroid 0.50006103515625 0.50006103515625 NaN NaN
first point 0 0 NaN NaN
last point 0 0 NaN NaN
extent 0 0 1.0001220703125 1.0001220703125 NaN NaN NaN NaN
hull 6.12398146082414E-17 1.0001220703125 1.0001220703125 1.0001220703125 1.0001220703125 0 0 0


Rectangle  Coordinates:  [[0, 0], [0, 1], [2, 1], [2, 0]]
X   0.0000000000000000  Y   0.0000000000000000
X   0.0000000000000000  Y   1.0000000000000000
X   2.0000000000000000  Y   1.0000000000000000
X   2.0000000000000000  Y   0.0000000000000000
Polygon properties:
area   2.0003662258386612
length 6.00048828125
centroid 1.00006103515625 0.50006103515625 NaN NaN
true centroid 1.00006103515625 0.50006103515625 NaN NaN
first point 0 0 NaN NaN
last point 0 0 NaN NaN
extent 0 0 2.0001220703125 1.0001220703125 NaN NaN NaN NaN
hull 6.12398146082414E-17 1.0001220703125 2.0001220703125 1.0001220703125 2.0001220703125 -2.22044604925031E-16 0 0

close but no cigar, and not readily attributable to floating point representation (ie a unit square should yield an area of 1.0 or exceptionally close to it.

Results using more classic pure Python methods from:

'''
AreaCentroid.py

original source:  http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/python.txt

'''

def polygon_area(pnts):
  '''determine the area given a list of points'''
  area = 0.0
  n = len(pnts)
  j = n - 1
  i = 0
  for point in pnts:
    p0 = pnts[i]
    p1 = pnts[j]
    area += p0[0] * p1[1]
    area -= p0[1] * p1[0]
    j = i
    i += 1
  area /= 2.0
  return area

def polygon_centroid(pnts):
  '''return the centroid of a polygon'''
  n = len(pnts)
  x = 0.0
  y = 0.0
  j = n -1
  i = 0
  for point in pnts:
    p0 = pnts[i]
    p1 = pnts[j]
    f = p0[0] * p1[1] - p1[0] * p0[1]
    x += (p0[0] + p1[0]) * f
    y += (p0[1] + p1[1]) * f
    j = i
    i += 1
  f = polygon_area(pnts) * 6
  return [x/f, y/f]

import os, sys

triangle = [[0,0],[0,1],[1,0]]
square = [[0,0],[0,1],[1,1],[1,0]]
rectangle = [[0,0],[0,1],[2,1],[2,0]]
polygons = [triangle, square, rectangle]
labels = ["Triangle", "Square", "Rectangle"]

for i in range(len(polygons)):
  a_poly = polygons[i]
  print "\n", labels[i]," Coordinates: ", a_poly
  print "polygon area: ", polygon_area(a_poly)
  print "polygon centroid: ", polygon_centroid(a_poly)



yield correct results.
Triangle  Coordinates:  [[0, 0], [0, 1], [1, 0]]
polygon area:  0.5
polygon centroid:  [0.33333333333333331, 0.33333333333333331]

Square  Coordinates:  [[0, 0], [0, 1], [1, 1], [1, 0]]
polygon area:  1.0
polygon centroid:  [0.5, 0.5]

Rectangle  Coordinates:  [[0, 0], [0, 1], [2, 1], [2, 0]]
polygon area:  2.0
polygon centroid:  [1.0, 0.5]

Is this a bug? and can anyone else confirm?

This also extends to the polyline class as shown in the following
'''
LineDemo.py

Demonstrates calculating various properties for polylines

'''
import arcpy
pnt = arcpy.Point()

simple = [[0,0],[1,1]]
two_piece = [[0,0],[0,1],[1,1]]
three_piece = [[0,0],[0,1],[1,1]]
polylines = [simple, two_piece, three_piece]
labels = ["simple", "two_piece", "three_piece"]
array = arcpy.Array()
polys = []
for i in range(len(polylines)):
  a_poly = polylines[i]
  print "\n", labels[i]," Coordinates: ", a_poly
  for pair in a_poly:
    pnt.X = pair[0]
    pnt.Y = pair[1]
    array.add(pnt)
    print " X: %20.16f  Y: %20.16f" % (pnt.X, pnt.Y)
  poly = arcpy.Polyline(array)
  print "Polyline properties:"
  print " length %20.16f" % (poly.length)
  print " centroid %s\n true centroid %s " % (poly.centroid, poly.trueCentroid) 
  print " first point %s\n last point %s " % (poly.firstPoint, poly.lastPoint)
  print " extent %s " % (poly.extent)
  print " hull %s\n " % (poly.hullRectangle)
  array.removeAll()
  polys.append(poly)  #for further use


which for a simple two point line yields

simple  Coordinates:  [[0, 0], [1, 1]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   1.0000000000000000  Y:   1.0000000000000000
Polyline properties:
length   1.4143861958645956
centroid 0.50006103515625 0.50006103515625 NaN NaN
true centroid 0.50006103515625 0.50006103515625 NaN NaN
first point 0 0 NaN NaN
last point 1.0001220703125 1.0001220703125 NaN NaN
extent 0 0 1.0001220703125 1.0001220703125 NaN NaN NaN NaN
hull 0 0 0 0 1.0001220703125 1.0001220703125 1.0001220703125 1.0001220703125

as output.  The only common thread seems to be the extent/last point is in error according to

>>> print 1 + (1.0/2.0**13)
1.00012207031

so is something being added to the coordinate(s)?  I will have to report this thread as a bug, I presume

Outcomes