# Thread: Errors in arcpy's Polygon class

1. ## Errors in arcpy's Polygon class

This code using arcpy's polygon class to calculate simple properties for simple shapes are not accurate.
Code:
```'''
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]
print " X %20.16f  Y %20.16f" % (pnt.X, pnt.Y)
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:

Code:
```'''
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
Code:
```'''
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]
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

2. ## Re: Errors in arcpy's Polygon class

Here are some hints, Dan:

1. Apropos the numbers in "hull 1.0001220703125 2.22044604925031E-16 0.50006103515625 -0.50006103515625 -0.50006103515625 0.50006103515625 0 1.0001220703125":
1.0001220703125 = 1 + 2^(-13)
2.22044604925031E-16 is double-precision error (most likely)
0.50006103515625 = (1 + 2^(-13)) / 2
2. Apropos the numbers in "centroid 0.3333740234375 0.3333740234375...":
0.3333740234375 = (1 + 2^(-13)) / 3
3. Apropos the number in "area 0.5001220777630806":
0.5001220777630806 = 1/2 + 2^(-13) + 2^(-27) = (1 + 2^(-13))^2 / 2
4. Apropos the number in "length 3.41463033649":
3.41463033649 = (1 + 2^(-13)) + (1 + 2^(-13)) + Sqrt(2) * (1 + 2^(-13)) + powers of 2 smaller than 10^(-13).
(I use "=" in the sense of exactly equal to, not just approximately equal to!)

From these we can glean that ArcGIS is likely to be representing your original floating coordinates as scaled integers within a discrete space wherein 2^13 = 8192 units correspond to one of your units, and in so doing it has rounded some of the coordinates outwards. E.g., 1.00 looks like it got internally represented 8193/8192. (I suspect the zeros were accurately represented.)

One suspects an (internal, invisible) error in floating point to integer conversion has occurred.

As a check, what happens when you scale all your coordinates up to, say, 1024.00 = 2^10 instead of 1.00? Is the error still 2^(-13) or does it scale up to 2^(-3) = 0.125 also? Regardless, this discrepancy seems to obliterate 19 binary digits from the expected 32-bit precision of an integer coordinate and almost 40 binary digits from a true double precision coordinate. In base 10, these errors translate to six and 13 orders of magnitude, respectively: they're huge!

Even if this is not technically an "error" and turns out to be our misinterpretation of how Python and ArcGIS interact, it seems like it's going to be the predominant interpretation by even expert programmers, as your code so clearly shows. In a system (to borrow from the theme of this year's UC) that works hard to reproject points to submillimeter accuracy and routinely is used to position life- and mission-critical things to within a few meters or better (e.g., military targets, surveyed boundaries, buried gas and electric power lines), a positional discrepancy this large has to be considered a show-stopper. It needs immediate investigation and a clear, convincing resolution.

3. ## Re: Errors in arcpy's Polygon class

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
Code:
```'''
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[i]
print "\n", " Coordinates: ", a_poly
for pair in a_poly:
pnt.X = pair[0]
pnt.Y = pair[1]
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

4. ## Re: Errors in arcpy's Polygon class

Dan, 0.0010986328125000 equals 0.001 rounded up to the nearest multiple of 2^(-13).

If these numbers were decimal degrees, the rounding would be creating up to a 14 meter error. For a pair of coordinates, that would introduce a net positional shift between 0 and 20 meters (varying from point to point). In effect, the arcpy procedures appear to be snapping all points onto a 0.44 arcsecond grid. (0.44 = 3600/2^13.) That might be ok at scales of 1:50,000 or smaller but obviously is unacceptable when you zoom in more than that or when you do precision analyses.

Have you reported this yet?

5. ## Re: Errors in arcpy's Polygon class

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.

6. ## Re: Errors in arcpy's Polygon class

This potentially goes way beyond geoprocessing, Dan. How do we know this problem doesn't lie in some core component of the entire system?

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.

7. ## Re: Errors in arcpy's Polygon class

Originally Posted by whuber
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

8. ## Re: Errors in arcpy's Polygon class

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

9. ## Re: Errors in arcpy's Polygon class

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

10. ## Re: Errors in arcpy's Polygon class

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

Code:
```'''
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[i]
print "\n", " Coordinates: ", a_poly
for pair in a_poly:
pnt.X = pair[0]
pnt.Y = pair[1]
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?"
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?
distance between last of 0 and 1st of 1 0.0

11. ## Re: Errors in arcpy's Polygon class

Originally Posted by Dan_Patterson
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

12. ## Re: Errors in arcpy's Polygon class

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.
Code:
```'''
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[i]
print "\n", " Coordinates: ", a_poly
for pair in a_poly:
pnt.X = pair[0]
pnt.Y = pair[1]
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

13. ## Re: Errors in arcpy's Polygon class

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.

Code:
```# 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

14. ## Re: Errors in arcpy's Polygon class

Originally Posted by dwynne
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.
The Nimbus ID is NIM059845

15. ## Re: Errors in arcpy's Polygon class

David
One would expect high precision as the default... perhaps next time, in the interim

poly = arcpy.Polyline(array,arcpy.SpatialReference())

seemed to work for polylines and I presume polygons etc...
perhaps, documentation on, arcpy.SpatialReference would be useful or go the high route by default. If one isn't creating any output geometry (next lesson) the SR is not necessary. A fix would be nice...thanks for getting on this.

16. ## Re: Errors in arcpy's Polygon class

For polygons, the spatial reference (SR) needs to be explicit. The original GeometryErrorsDemo.py is shown below with an example SR and its new output for a triangle, square and rectangle
Code:
```'''
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 = float(pair[0])
pnt.Y = float(pair[1])
print " X %20.16f  Y %20.16f" % (pnt.X, pnt.Y)
#poly = arcpy.Polygon(array)   #low precision
#
SR = arcpy.SpatialReference("c:/temp/NAD 1983 UTM Zone 18N.prj")
#
poly = arcpy.Polygon(array, SR)
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 below

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.5000000000000000
length 3.41421356237
centroid 0.333333333333333 0.333333333333333 NaN NaN
true centroid 0.333333333333333 0.333333333333333 NaN NaN
first point 0 0 NaN NaN
last point 0 0 NaN NaN
extent 0 0 1 1 NaN NaN NaN NaN
hull 1 2.22044604925031E-16 0.5 -0.5 -0.5 0.5 0 1

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.0000000000000000
length 4.0
centroid 0.5 0.5 NaN NaN
true centroid 0.5 0.5 NaN NaN
first point 0 0 NaN NaN
last point 0 0 NaN NaN
extent 0 0 1 1 NaN NaN NaN NaN
hull 6.12323399573677E-17 1 1 1 1 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.0000000000000000
length 6.0
centroid 1 0.5 NaN NaN
true centroid 1 0.5 NaN NaN
first point 0 0 NaN NaN
last point 0 0 NaN NaN
extent 0 0 2 1 NaN NaN NaN NaN
hull 6.12323399573677E-17 1 2 1 2 -2.22044604925031E-16 0 0

17. ## Re: Errors in arcpy's Polygon class

Originally Posted by dwynne
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.
I've had some offline discussions recently about the best way to create a SpatialReference object for these circumstances; specifically regarding IsLowPrecision, IsHighPrecision.

If you're describing data, most data is now high precision (unless you're working with an earlier version personal geodatabase or SDE geodatabase). So if describing a 9.1 personal geodatabase feature class, and using it's spatialreference property, and you would get a low precision geometry. You could try and manipulate the SpatialReference object, but adjusting a single property on a SpatialReference is fraught with perils (and likewise manipulating the string equivalent of the SpatialReference). Its not as easy as just modifying a property, as many properties of the SpatialReference are interconnected. A better way of ensuring HighPrecision would be creating a spatial reference not via data source, but a .prj file.

Code:
`sr = arcpy.SpatialReference(r"C:\Program Files\ArcGIS\Desktop10.0\Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj")`

-Dave

18. ## Re: Errors in arcpy's Polygon class

I've always used the text contained in the .prj file. That way you don't have to rely on having a valid path to the system-stored .prj files, that, for 64 bit OS, are stored in a different directory that a 32 bit OS!

19. ## Re: Errors in arcpy's Polygon class

I've always used the text contained in the .prj file. That way you don't have to rely on having a valid path to the
What I normally do to make my own code adaptable for location of the installed .prj files (or for that matter any of the installed files), is to use GetInstallInfo like below:

Code:
```prjFile = os.path.join(arcpy.GetInstallInfo()['InstallDir'],
"Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj")
sr = arcpy.SpatialReference(prjFile)```
-Dave

20. ## Re: Errors in arcpy's Polygon class

NIM059845 didn't make it into 10 SP1 and there is no indication as to when it is going to be fixed, nor does the NIM site indicate the temporary workaround. Any further information?

21. ## Re: Errors in arcpy's Polygon class

Originally Posted by Dan_Patterson
NIM059845 didn't make it into 10 SP1 and there is no indication as to when it is going to be fixed, nor does the NIM site indicate the temporary workaround. Any further information?
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

Thanks David

23. ## Re: Errors in arcpy's Polygon class

Thank you David for the update.

24. ## Re: Errors in arcpy's Polygon class

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 Code:
``` geom = arcpy.SearchCursor(polygonFC, "FID = 0").next().getValue("Shape") geomGeoJSON = geom.__geo_interface__ geom2 = arcpy.AsShape(geomGeomJSON) geom.equals(geom2)     # This fails  ```
Sample representations of the geometry (note at least their representations seem identical):
PHP Code:
```     >>> 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)]]}  ```
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 Code:
``` 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]))  ```
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.

25. ## Re: Errors in arcpy's Polygon class

Originally Posted by dwynne
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?