topic Re: Errors in arcpy's Polygon class in Geoprocessing Questions

Errors in arcpy's Polygon class

DanPatterson_Retired — Thu, 05 Aug 2010 09:55:39 GMT

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

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Sat, 11 Dec 2021 19:04:10 GMT

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 🙂

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Thu, 05 Aug 2010 19:47:00 GMT

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.

Re: Errors in arcpy's Polygon class

ChrisFox3 — Thu, 05 Aug 2010 21:44:58 GMT

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

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Thu, 05 Aug 2010 22:17:05 GMT

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

Re: Errors in arcpy's Polygon class

ChrisSnyder — Thu, 05 Aug 2010 22:53:49 GMT

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

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Sat, 11 Dec 2021 19:04:12 GMT

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

Re: Errors in arcpy's Polygon class

ChrisFox3 — Thu, 05 Aug 2010 23:18:40 GMT

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

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Sat, 11 Dec 2021 19:04:15 GMT

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

Re: Errors in arcpy's Polygon class

DavidWynne — Sat, 11 Dec 2021 19:04:18 GMT

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

Re: Errors in arcpy's Polygon class

DavidWynne — Mon, 09 Aug 2010 23:57:39 GMT

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

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Tue, 10 Aug 2010 03:07:18 GMT

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.

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Sat, 11 Dec 2021 19:04:20 GMT

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

'''
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 = float(pair[0])
    pnt.Y = float(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)   #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

Re: Errors in arcpy's Polygon class

DavidWynne — Wed, 11 Aug 2010 20:57:28 GMT

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.

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

-Dave

Re: Errors in arcpy's Polygon class

ChrisSnyder — Wed, 11 Aug 2010 21:14:04 GMT

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!

sr = 'PROJCS["NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",45.83333333333334],PARAMETER["Standard_Parallel_2",47.33333333333334],PARAMETER["Latitude_Of_Origin",45.33333333333334],UNIT["Foot_US",0.3048006096012192]]'

Re: Errors in arcpy's Polygon class

DavidWynne — Sat, 11 Dec 2021 19:04:23 GMT

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:

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

-Dave

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Sat, 13 Nov 2010 12:50:15 GMT

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?

Re: Errors in arcpy's Polygon class

DavidWynne — Thu, 18 Nov 2010 16:28:18 GMT

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

Re: Errors in arcpy's Polygon class

DanPatterson_Retired — Thu, 18 Nov 2010 19:31:57 GMT

Thanks David

Re: Errors in arcpy's Polygon class

TedCronin — Thu, 18 Nov 2010 21:01:49 GMT

Thank you David for the update.