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

[DanPatterson_Retired]

Yes it gets worse, even a two point line's properties can't be determined properly whether they are along the X or Y axis relative to the 0,0 origin...another example script and its output follow

'''
LineDemo2.py

Demonstrates errors in calculating various properties for polylines

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

polylines = []
vals = [0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0]
for a_val in vals:
  a_line = [[0.0,0.0],[a_val,0.0]]
  polylines.append(a_line)
array = arcpy.Array()
polys = []
for i in range(len(polylines)):
  a_poly = polylines
  print "\n", " 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

with the results in changing x from 0.001 to 1000 by a factor of 10

Coordinates:  [[0.0, 0.0], [0.001, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   0.0010000000000000  Y:   0.0000000000000000
Polyline properties:
length   0.0010986328125000
centroid 0.00054931640625 0 NaN NaN
true centroid 0.00054931640625 0 NaN NaN
first point 0 0 NaN NaN
last point 0.0010986328125 0 NaN NaN
extent 0 0 0.0010986328125 0 NaN NaN NaN NaN

Coordinates:  [[0.0, 0.0], [0.01, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   0.0100000000000000  Y:   0.0000000000000000
Polyline properties:
length   0.0100708007812500
centroid 0.005035400390625 0 NaN NaN
true centroid 0.005035400390625 0 NaN NaN
first point 0 0 NaN NaN
last point 0.01007080078125 0 NaN NaN
extent 0 0 0.01007080078125 0 NaN NaN NaN NaN

Coordinates:  [[0.0, 0.0], [0.10000000000000001, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   0.1000000000000000  Y:   0.0000000000000000
Polyline properties:
length   0.1000976562500000
centroid 0.050048828125 0 NaN NaN
true centroid 0.050048828125 0 NaN NaN
first point 0 0 NaN NaN
last point 0.10009765625 0 NaN NaN
extent 0 0 0.10009765625 0 NaN NaN NaN NaN

Coordinates:  [[0.0, 0.0], [1.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   1.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length   1.0001220703125000
centroid 0.50006103515625 0 NaN NaN
true centroid 0.50006103515625 0 NaN NaN
first point 0 0 NaN NaN
last point 1.0001220703125 0 NaN NaN
extent 0 0 1.0001220703125 0 NaN NaN NaN NaN

Coordinates:  [[0.0, 0.0], [10.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:  10.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length  10.0001220703125000
centroid 5.00006103515625 0 NaN NaN
true centroid 5.00006103515625 0 NaN NaN
first point 0 0 NaN NaN
last point 10.0001220703125 0 NaN NaN
extent 0 0 10.0001220703125 0 NaN NaN NaN NaN

Coordinates:  [[0.0, 0.0], [100.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X: 100.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length 100.0001220703125000
centroid 50.0000610351563 0 NaN NaN
true centroid 50.0000610351563 0 NaN NaN
first point 0 0 NaN NaN
last point 100.000122070313 0 NaN NaN
extent 0 0 100.000122070313 0 NaN NaN NaN NaN

Coordinates:  [[0.0, 0.0], [1000.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X: 1000.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length 1000.0001220703125000
centroid 500.000061035156 0 NaN NaN
true centroid 500.000061035156 0 NaN NaN
first point 0 0 NaN NaN
last point 1000.00012207031 0 NaN NaN
extent 0 0 1000.00012207031 0 NaN NaN NaN NaN

so percentage wise, the greater the percentage error as you suggest...hope this gets fixed before the "Power of ArcPy" lecture gets delivered 🙂

[DanPatterson_Retired]

I have reported this indirectly through Jim Barry since faculty for our site license, don't have the "magic numbers" to report bugs and I can't find any generic "bug report" site that allows users (and not administrators) to post bugs....perhaps Jim could add MVP's and frequent contributors etc to post bugs, or provide sites without a "magic number" to report bugs.  Perhaps putting a "bug" icon on the forum threads could serve as a direct link to notify the correct group, in this case, Geoprocessing.

And on an aside, it would be inappropriate for any of these metrics to apply to geographic coordinates, which is even more disturbing.

I will investigate other issues further and try the long slog through the official channels if no one here picks up on this...but lecture and manual notes can't wait until SP1 arrives, I am probably going to have to demonstrate the principles outside of the arcpy environment.

[ChrisFox3]

I searched and eventually found solid confirmation of your failure: there is no way for mere mortals to report a bug, no matter how big or how pressing it might be.  I am stunned.

Anyone can report a software defect to Esri regardless of maintenance status. To do so you can contact Esri Support through any of our media channels; phone, e-mail or chat. The most efficient way would be to go to http://support.esri.com and click "New Support Request". You will then be asked to enter some information to validate you as a customer, but don't let this deter you as even if you are not an authorized maintenance caller you will still be able to submit a defect. On the next page there is a submission type option where you can select "Bug".

Once you submit the defect report it will be routed to a Technical Support Analyst who will work with you to reproduce the issue and if it is confirmed to be a defect will go ahead and submit it to the Development team for review.

-Chris

[DanPatterson_Retired]

That route sent me to ESRI Canada, who would have to process it and forward it on to ESRI Redlands.  I have never been able to file a bug report without the "magic" customer numbers which are housed within technical support somewhere within our university. 

There has got to be a simpler solution to get to Redlands without the corporate/country travel stuff.  I have to unfortunately had to rely on "insiders" grace for a number of years.  As suggested, perhaps, MVP/Frequent Contributors/Pests be given a "Go Away" number to report bugs, particularly when they are more important that a color palette discrepencies and the like.  Pardon my frustration, but as you can see, the important issue being brought forward within this thread has been relegated to a lower status than that of how to report bugs.  Should you wish to contact me directly, I would appreciate it, in the interim, can the major issue here please be addressed
Regards

[ChrisSnyder]

Hasn't this always been like this though? Even in v9.3 if you manually create a square having coordinate pairs of:
0,0
0,2
2,2
2,0

It reports the centroid as being:
1.000061, 1.000061

and the area as:
4.000488

and the perimeter as:
8.000488

Could be wrong here, but I was led to believe what Bill hinted at: ESRI's data models store/process coordinates as integers rather that double precision floats (seemingly to save space and speed processing).

Seems like the larger the feature (relative distances) that, at least in v9.3, things become more "precise". For example, if you make a larger square of:
0,0
0,100000
100000,100000
100000,0
the centroid is now "exactly":
50000, 50000
and the area:
10000000000

[DanPatterson_Retired]

That means that a solution, if that was the case, hasn't been provided after all these years.  I am trying to demonstrate to students that the new powers of GIS and its incorporation of Python within its framework will make a simpler work environment for them...so I gave a simple demonstration of calculating the area of a unit square (etc) using the new "arcpy" environment and all its useful utilities, then I go on to show them how to do it in other environments, particular, pure Python...I had hoped to show that if you have ArcGIS installed use "arcpy" to get basic metric properties without the need to code your own procedures or use other modules (like Shapely etc)...unfortunately, those notes are burning on the bonfire outback.

The more important issue is that there are many instances where feature geometry is manipulated (ie translated, rotated and scaled) to facilitate determinations of other geometric properties such as minimum area bounding circles, convex hulls, minimum area bounding ellipses, densifying features etc etc) since working in native projected coordinates is not appropriate (for a variety of reasons).  One would expect that any software would allow one to compute simple metrics for standard shapes centered at the origin and/or the minimum extent being at the origin.  If this is not the case (I never found an issue with this in AV3.x and/or calculations in 9.x that interfaced with Python), then a clear statement should be made that if you aren't using large coordinates, then you are going to get large errors. 

ESRI has made the move to Python as your field calculator language de jour (VB Script as well) and ESRI strongly supports NumPy for a lot of the SA stuff (I presume, but I am not finished testing), so why can't I get a simple area for a unit square with the lower left extent at the origin from a module called "ARCPY", I thought, perhaps erroneously, that a synergy had developed...please someone, fix the issue, rather than defending it or bringing up the arcinfo past or call the module "ARCPYish"...
rant over...thanks for the comments Chris

PS
Don't take offense, but this has serious ramifications for anyone that wants to use "a" GIS for any computational purposes.  The limitations have to be spelled out or addressed, there is no "silent" option


By the way intersection of two polylines doesn't work for simply line features either

'''
IntersectDemo.py

Demonstrates errors in calculating various properties for polylines

'''
import arcpy
pnt = arcpy.Point()
#
#NOTE vary the epsilon parameter 1.0/8192 = 1.0/(2**13
#  Try 1.0/2.0**14 etc etc
epsilon = 1.0/(2.0**13)
polylines = [[[0.0,0.0], [1.0,0.0]], [[1.0 + epsilon, 0.0],[2.0,0.0]]]
array = arcpy.Array()
polys = []
for i in range(len(polylines)):
  a_poly = polylines
  print "\n", " 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)
print "\n Does polyline 0 touch polyline 1?"
print "answer ", polys[0].touches(polys[1])
inter_pnt_dist = polys[0].lastPoint.X - polys[1].firstPoint.X
print "distance between last of 0 and 1st of 1 ", inter_pnt_dist

results, vary the epsilon value to see different outcomes.  The results are wrong,
the line segments should be separated by epsilon instead of 0.0

>>>
Coordinates:  [[0.0, 0.0], [1.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   1.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length   1.0001220703125000
centroid 0.50006103515625 0 NaN NaN
true centroid 0.50006103515625 0 NaN NaN
first point 0 0 NaN NaN
last point 1.0001220703125 0 NaN NaN
extent 0 0 1.0001220703125 0 NaN NaN NaN NaN

Coordinates:  [[1.0001220703125, 0.0], [2.0, 0.0]]
X:   1.0001220703125000  Y:   0.0000000000000000
X:   2.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length   1.0000000000000000
centroid 1.5001220703125 0 NaN NaN
true centroid 1.5001220703125 0 NaN NaN
first point 1.0001220703125 0 NaN NaN
last point 2.0001220703125 0 NaN NaN
extent 1.0001220703125 0 2.0001220703125 0 NaN NaN NaN NaN

Does polyline 0 touch polyline 1?
answer  True
distance between last of 0 and 1st of 1  0.0

[ChrisFox3]

Should you wish to contact me directly, I would appreciate it, in the interim, can the major issue here please be addressed

Hi Dan,

I would really like to talk to you about this and see what we can do to address the technical issue in the short term and the process for submitting defects in the long term. Unfortunately, I don't have your contact information, could you contact me at chris_fox@esri.com.

Thanks,
-Chris

[DanPatterson_Retired]

Contacted by ESRI (esri?), still awaiting the final means of officially reporting the bug, I can report that the simple line example is also directional, for example whether a line radiates outwards from the 0,0 origin in the positive direction, inward toward the origin and their complements in the negative direction along the x-axis (outward and inward).  Haven't finished varying along the y-axis or at angles, but the following code will allow others to investigate on their own.

'''
LineDemo2.py

Demonstrates errors in calculating various properties for polylines

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

polylines = []
#vals = [0.001, 0.01, 0.1, 1.0, math.sqrt(2.0), 2.0]
vals = [0.001, 0.01, 1.0 ]
for a_val in vals:
  polylines.append( [[0.0,0.0],[a_val,0.0]] )  #positive out
  polylines.append( [[0.0,0.0],[-a_val,0.0]] ) #negative out
  polylines.append( [[a_val,0.0],[0.0,0.0]] )  #positive in
  polylines.append( [[-a_val,0.0],[0.0,0.0]] ) #negative in
array = arcpy.Array()
polys = []
for i in range(len(polylines)):
  a_poly = polylines
  print "\n", " 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)
  cal_leng = poly.length
  real_leng = max(abs(a_poly[1][0]), abs(a_poly[0][0]))
  diff_leng = cal_leng - real_leng
  print "Polyline properties:"
  print " length:        %20.16f" % (cal_leng)
  print " true length:   %20.16f" % (real_leng)
  print " difference:    %20.16f" % (diff_leng)
  print " percent error: %20.16f" % (diff_leng/real_leng * 100)
  array.removeAll()
  polys.append(poly)  #for further use

results for 3 simple lengths, outward from the origin, inward towards the origin in both the positive and negative x directions

>>>
Coordinates:  [[0.0, 0.0], [1.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:   1.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length:          1.0001220703125000
true length:     1.0000000000000000
difference:      0.0001220703125000
percent error:   0.0122070312500000

Coordinates:  [[0.0, 0.0], [-1.0, 0.0]]
X:   0.0000000000000000  Y:   0.0000000000000000
X:  -1.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length:          1.0000000000000000
true length:     1.0000000000000000
difference:      0.0000000000000000
percent error:   0.0000000000000000

Coordinates:  [[1.0, 0.0], [0.0, 0.0]]
X:   1.0000000000000000  Y:   0.0000000000000000
X:   0.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length:          1.0001220703125000
true length:     1.0000000000000000
difference:      0.0001220703125000
percent error:   0.0122070312500000

Coordinates:  [[-1.0, 0.0], [0.0, 0.0]]
X:  -1.0000000000000000  Y:   0.0000000000000000
X:   0.0000000000000000  Y:   0.0000000000000000
Polyline properties:
length:          1.0000000000000000
true length:     1.0000000000000000
difference:      0.0000000000000000
percent error:   0.0000000000000000

the gp group is on this and I will report when I can report

[DavidWynne]

Hi Dan,
I appreciate your concern, and thanks for bringing light to the problem.  Looks like the geometry is essentially low-precision when created in this way.  As a workaround, I was able to get much better precision if I applied a SpatialReference object to the class' 2nd parameter.  Not ideal of course, but better than the alternative.

# sr being a SpatialReference
poly = arcpy.Polyline(array, sr)

I will update this thread shortly when I have a proper Nimbus ID to track on.  I'm marking this issue as candidate for 10.0 sp1.

-Dave