21 Replies Latest reply on Nov 8, 2012 7:40 AM by michalis

    Errors in arcpy's Polygon class

    Dan_Patterson
      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[i]
        print "\n", labels[i]," 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[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
      '''
      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]
          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
          Dan_Patterson
          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[i]
            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
            Dan_Patterson
            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
              cfox-esristaff
              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
                Dan_Patterson
                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
                  csny490
                  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
                    Dan_Patterson
                    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[i]
                      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
                      cfox-esristaff
                      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
                        Dan_Patterson
                        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[i]
                          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
                          dwynne-esristaff
                          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
                            dwynne-esristaff
                            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
                              Dan_Patterson
                              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
                                Dan_Patterson
                                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[i]
                                  print "\n", labels[i]," 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
                                  dwynne-esristaff
                                  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
                                    csny490
                                    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
                                      dwynne-esristaff
                                      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
                                        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?
                                        • Re: Errors in arcpy's Polygon class
                                          dwynne-esristaff
                                          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
                                            Riverside
                                            Thank you David for the update.
                                            • Re: Errors in arcpy's Polygon class
                                              michalis
                                              Has there actually been a fix for this in 10.0? I am on 10.0 SP5 here, and I have a failure that probably relates, but I think it is the way the geometries.py file is built (using the arcobjectconversion and _base - perhaps it actually goes even deeper into the core C++ code). Note this failure is a simplified version as shown by Victor velarde at gis.stackexchange. Basically, the geometry exported by arcpy cannot be imported as the same on by the same session of arcpy.

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

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

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

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

                                              Also, a side-note: Will we ever get a [PYTHON] tag for highlighting? The vBulletin forums have plenty of examples for syntax highlighting in multiple languages.
                                              • Re: Errors in arcpy's Polygon class
                                                michalis
                                                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?