Errors in arcpy's Polygon class

[DanPatterson_Retired]

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   print "\n", labels," 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     p1 = pnts     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     p1 = pnts     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   print "\n", labels," 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   print "\n", labels," 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

[MichalisAvraam]

Has there actually been a fix for this in 10.0? I am on 10.0 SP5 here, and I have a failure that probably relates, but I think it is the way the geometries.py file is built (using the arcobjectconversion and _base - perhaps it actually goes even deeper into the core C++ code). Note this failure is a simplified version as shown by Victor velarde at gis.stackexchange. Basically, the geometry exported by arcpy cannot be imported as the same on by the same session of arcpy.

[PHP]
geom = arcpy.SearchCursor(polygonFC, "FID = 0").next().getValue("Shape")
geomGeoJSON = geom.__geo_interface__
geom2 = arcpy.AsShape(geomGeomJSON)
geom.equals(geom2)     # This fails
[/PHP]

Sample representations of the geometry (note at least their representations seem identical):
[PHP]
    >>> geom
    {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
    >>> geom2
    {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
[/PHP]

Incidentally, the __geo_interface__ fails miserably if the polygon has a hole, and you need to change the geometries.py to actually handle the holes yourself (the Polygon class). Just in case anyone notices at ESRI and wants to release this as a patch for 10.0:

Change the Polygon class' return values for the following methods:
[PHP]
def __geo_interface__(self):
    return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]}
def _fromGeoJson(cls, data):
    ...
    return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates]))
[/PHP]

Also, a side-note: Will we ever get a [PYTHON] tag for highlighting? The vBulletin forums have plenty of examples for syntax highlighting in multiple languages.

[MichalisAvraam]

Sorry Dan, I was disappointed myself.  We have had it fixed in our 10.1 code base, and we will be pushing it through to service pack 2 in the next short while.  The projected date for SP2 I believe is still early 2011.

The NIM report has also been updated, and should get pushed to the web on the next flush.

-Dave

Dave,

It seems like 10 SP 5 still doesn't have a fix. Can you verify that?

[FC_Basson]

Same issue still present in 10.5.1, but applying a Spatial Reference as suggested by David Wynne‌ solved it for me.

[DanPatterson_Retired]

I know...  it still applies