# Errors in arcpy's Polygon class

**Dan_Patterson**Aug 5, 2010 2:55 AM

This code using arcpy's polygon class to calculate simple properties for simple shapes are not accurate.

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:

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

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

''' 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