+ Reply to Thread
Results 1 to 25 of 25

Thread: Errors in arcpy's Polygon class

  1. #1
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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]
        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:

    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]
        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
    Last edited by Dan_Patterson; 08-05-2010 at 07:01 AM.
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  2. #2
    William Huber

    Join Date
    Apr 2010
    Posts
    694
    Points
    73
    Answers Provided
    2


    1

    Thumbs up 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.
    --Bill Huber
    Quantitative Decisions
    For more help, visit the worldwide community at http://gis.stackexchange.com

  3. #3
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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]
        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
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  4. #4
    William Huber

    Join Date
    Apr 2010
    Posts
    694
    Points
    73
    Answers Provided
    2


    1

    Default 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?
    --Bill Huber
    Quantitative Decisions
    For more help, visit the worldwide community at http://gis.stackexchange.com

  5. #5
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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.
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  6. #6
    William Huber

    Join Date
    Apr 2010
    Posts
    694
    Points
    73
    Answers Provided
    2


    1

    Default 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.
    --Bill Huber
    Quantitative Decisions
    For more help, visit the worldwide community at http://gis.stackexchange.com

  7. #7
    Chris Fox

    Join Date
    Oct 2009
    Posts
    523
    Points
    302
    Answers Provided
    36


    1

    Default Re: Errors in arcpy's Polygon class

    Quote Originally Posted by whuber View Post
    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. #8
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    0

    Default 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
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  9. #9
    Chris Snyder

    Join Date
    May 2010
    Posts
    1,263
    Points
    418
    Answers Provided
    38


    1

    Default 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. #10
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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]
        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
    Last edited by Dan_Patterson; 08-05-2010 at 04:13 PM.
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  11. #11
    Chris Fox

    Join Date
    Oct 2009
    Posts
    523
    Points
    302
    Answers Provided
    36


    0

    Default Re: Errors in arcpy's Polygon class

    Quote Originally Posted by Dan_Patterson View Post
    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. #12
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    0

    Default 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]
        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
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  13. #13
    David Wynne
    Join Date
    Oct 2009
    Posts
    97
    Points
    42
    Answers Provided
    3


    0

    Default 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. #14
    David Wynne
    Join Date
    Oct 2009
    Posts
    97
    Points
    42
    Answers Provided
    3


    0

    Default Re: Errors in arcpy's Polygon class

    Quote Originally Posted by dwynne View Post
    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. #15
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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.
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  16. #16
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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])
        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
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  17. #17
    David Wynne
    Join Date
    Oct 2009
    Posts
    97
    Points
    42
    Answers Provided
    3


    1

    Default Re: Errors in arcpy's Polygon class

    Quote Originally Posted by dwynne View Post
    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
    Last edited by dwynne; 08-11-2010 at 02:41 PM.

  18. #18
    Chris Snyder

    Join Date
    May 2010
    Posts
    1,263
    Points
    418
    Answers Provided
    38


    0

    Default 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!

    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]]'
    Last edited by csny490; 08-11-2010 at 02:18 PM.

  19. #19
    David Wynne
    Join Date
    Oct 2009
    Posts
    97
    Points
    42
    Answers Provided
    3


    0

    Default Re: Errors in arcpy's Polygon class

    Quote Originally Posted by chris.snyder@wadnr.gov View Post
    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. #20
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default 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?
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  21. #21
    David Wynne
    Join Date
    Oct 2009
    Posts
    97
    Points
    42
    Answers Provided
    3


    1

    Default Re: Errors in arcpy's Polygon class

    Quote Originally Posted by Dan_Patterson View Post
    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
    Last edited by dwynne; 11-18-2010 at 12:56 PM. Reason: wrong year

  22. #22
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,756
    Points
    443
    Answers Provided
    41


    1

    Default Re: Errors in arcpy's Polygon class

    Thanks David
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  23. #23
    Ted Cronin

    Join Date
    Oct 2009
    Posts
    3,151
    Points
    29
    Answers Provided
    0


    1

    Default Re: Errors in arcpy's Polygon class

    Thank you David for the update.
    Regards,

    Ted C
    County of Riverside - Assessor/Clerk/Recorder
    ValueGIS

  24. #24
    Michalis Avraam
    Join Date
    May 2010
    Posts
    21
    Points
    0
    Answers Provided
    0


    0

    Default 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.842348155999947.060497293000083), (-122.8423975559999247.059262423000064), (-122.8441691359998947.059309693000046), (-122.8441691359998947.060497293000083), (-122.842348155999947.060497293000083)]]}
        >>> 
    geom2
        
    {'type''Polygon''coordinates': [[(-122.842348155999947.060497293000083), (-122.8423975559999247.059262423000064), (-122.8441691359998947.059309693000046), (-122.8441691359998947.060497293000083), (-122.842348155999947.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.Xpt.Y) if pt else None) for pt in part] for part in self]}
    def _fromGeoJson(clsdata):
        ...
        return 
    cls(Array([map(lambda pPoint(*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. #25
    Michalis Avraam
    Join Date
    May 2010
    Posts
    21
    Points
    0
    Answers Provided
    0


    0

    Default Re: Errors in arcpy's Polygon class

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

+ Reply to Thread

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts