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