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

[DavidWynne]

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

[DanPatterson_Retired]

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.

[DanPatterson_Retired]

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

[DavidWynne]

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

[ChrisSnyder]

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

[DavidWynne]

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

[DanPatterson_Retired]

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?

[DavidWynne]

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

[DanPatterson_Retired]

Thanks David

[TedCronin]

Thank you David for the update.