Dan_Patterson

Errors in arcpy's Polygon class

Discussion created by Dan_Patterson Champion on Aug 5, 2010
Latest reply on Sep 11, 2018 by Dan_Patterson
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