weird issue using geometry objects in python

[AndrewHaley]

Hi,

I am writing a script in python to generate polygons from lists of coordinates that I have saved in csv files.

I have followed the proper method of assigning the coordinates as point objects, adding them to an array, and then "arcpy.Polygon(array)" to get it as a polygon object.

I am having a weird issue when I try to make the polygon from the array.

I have gone into the shell and step by stepped it, when i print my array, i get this:

<Array [<Point (147.7359, -35.5386, #, #)>, <Point (147.7366, -35.539, #, #)>, <Point (147.737, -35.539, #, #)>, <Point (147.7375, -35.539, #, #)>, <Point (147.7374, -35.5397, #, #)>, <Point (147.7368, -35.5401, #, #)>, <Point (147.7354, -35.5403, #, #)>, <Point (147.7343, -35.5405, #, #)>, <Point (147.732, -35.5408, #, #)>, <Point (147.7304, -35.5409, #, #)>, <Point (147.7299, -35.5414, #, #)>, <Point (147.7287, -35.5422, #, #)>, <Point (147.7272, -35.5423, #, #)>, <Point (147.726, -35.5433, #, #)>, <Point (147.7255, -35.5439, #, #)>, <Point (147.7247, -35.5445, #, #)>, <Point (147.724, -35.5452, #, #)>, <Point (147.7228, -35.5432, #, #)>, <Point (147.7235, -35.5421, #, #)>, <Point (147.7249, -35.5416, #, #)>, <Point (147.7256, -35.5411, #, #)>, <Point (147.7269, -35.541, #, #)>, <Point (147.7275, -35.541, #, #)>, <Point (147.7283, -35.5409, #, #)>, <Point (147.7287, -35.5408, #, #)>, <Point (147.7296, -35.5399, #, #)>, <Point (147.7308, -35.5396, #, #)>, <Point (147.7323, -35.5395, #, #)>, <Point (147.7333, -35.5392, #, #)>, <Point (147.734, -35.5387, #, #)>, <Point (147.7347, -35.5386, #, #)>, <Point (147.7359, -35.5386, #, #)>]>

which looks good to me, they are the exact coordinates from my csv, same order, and i don't require Z or M values so the "#" marks are expected.

when I use Polygon.getPart(), i get this:

<Array [<Array [<Point (147.727905273, -35.5413818359, #, #)>, <Point (147.725097656, -35.5441894531, #, #)>, <Point (147.723327637, -35.5426025391, #, #)>, <Point (147.725280762, -35.5413208008, #, #)>, <Point (147.727905273, -35.5413818359, #, #)>]>]>

I have done this numerous times, and i am definitely using the correct array in the arcpy.Polygon() function...

I should note that the polygon i am trying to draw is quite a simple shape also... no "donut holes", and its not multi-part...

I have added an attachment comparing the two in arcmap.

has anyone ever come across this before?

[NeilAyres · Accepted solution]

Andrew,
when you create the polygon object, be sure to specify the spatial reference of the object as well.
like :

polygeom = arcpy.Polygon(array, sr)

Usually I get the sr object from Describe of the feature class I am adding to, or you can use the WKID of the coord sys.
Cheers and good luck,
Neil

[JoshuaBixby]

Although this is an older thread/message, I stumbled onto this problem recently and noticed some odd behavior with the ArcPy Geometry Classes and spatial references.

As Neil Ayres‌ points out, including a spatial reference is a good practice.  The documentation of ArcPy Geometry Classes is a bit weak or confusing on this specific issue because spatial_reference is listed as an optional argument with a default value of None.  According to the documentation, there is no reason your original code shouldn't work.  What is happening might be related to a similar issue that was discussed in detail years back by Dan Patterson‌ and David Wynne in the Errors in arcpy's Polygon class thread.

I don't know what exactly is going on, but I would say there is a bug in the ArcPy Geometry Classes.  Looking at the same set of coordinates with four different spatial reference scenarios:

coords = [(147.7359, -35.5386), (147.7374, -35.5397), (147.7304, -35.5409),
          (147.7304, -35.5409), (147.7255, -35.5439), (147.7249, -35.5416)]

SR_list = [['noneSR', None],
           ['emptySR', arcpy.SpatialReference()],
           ['webmercSR', arcpy.SpatialReference(3857)],
           ['wgs84', arcpy.SpatialReference(4326)]]

results = [['SR_Name', 'Type', 'PCSName', 'GCSName', 'XYResolution', 'pointCount']]

for SR_name, SR in SR_list:
    pn = arcpy.Polygon(
        arcpy.Array([arcpy.Point(*coord) for coord in coords]),
        SR
    )
    pnSR = pn.spatialReference
    if pnSR is None:
        results.append([SR_name, SR.type, SR.PCSName,
                        SR.GCSName, SR.XYResolution, pn.pointCount])
    else:
        results.append([SR_name, pnSR.type, pnSR.PCSName,
                        pnSR.GCSName, pnSR.XYResolution, pn.pointCount])

for result in results:
    print '{:<12}{:<12}{:<36.34}{:<12.10}{:<12.10}{:<4}'.format(*result)

Yields:

SR_Name     Type        PCSName                             GCSName     XYResoluti  pointCount
noneSR      Unknown                                                     0.0001      0   
emptySR                                                                 0.0         6   
webmercSR   Projected   WGS_1984_Web_Mercator_Auxiliary_Sp              0.0001      0   
wgs84       Geographic                                      GCS_WGS_19  1e-09       6   

To me, this behavior seems odd.  Using no spatial reference (None), which is the default, arcpy.Polygon does not generate an error, but it effectively doesn't generate a polygon either (from a practical perspective, I don't consider a polygon with zero points a polygon, especially when it should have six points).  Using an empty spatial reference, arcpy.Polygon generates a 6-pt polygon using a spatial reference that has no type, no names, and no XY resolution.  Interestingly enough, although a polygon is generated with six points and all the Polygon methods work, none of the geoprocessing tools (CopyFeatures_management, etc...) that I tried would work with it, likely due to its nonexisting spatial reference.  Using the Web Mercator projection, arcpy.Polygon fails to create the polygon in question even though the coordinates are valid within the Web Mercator projection.  And finally using a geographic coordinate system, arcpy.Polygon generates the polygon we have been after the whole time.

From my perspective, the none/default and Web Mercator spatial references should generate a 6-pt polygon but don't while the empty spatial reference should generate an error instead of generating a 6-pt polygon.

[DanPatterson_Retired]

I can confirm that the behaviour I noticed several years ago still exists and your recommendation that a SR always be used still holds.  This not only affects the ability to generate geometric objects but also their 'precision'. 

User beware...this makes working with a fixed accuracy and precision difficult to control.  The errors may seem insignificant, with the exception of no features generated, but it affects working with more math specific modules such as numpy etc.  I have been comparing numpy generated point coordinates with those generated with arcpy and the spatial reference omission has always generated the differences.  I now use numpy recarrays to do most of my geometric operations using arcpy solely to transfer geometry and attribute information back and forth between the applications.

[JoshuaBixby]

Opened a case with Esri Support back on 12/19.  Esri Support was able to narrow the issue down even further, but the case remains open until Esri Development decides they have some time to look into it more.  I will post the bug number when it gets formally issued.

[JoshuaBixby]

Esri Support issued a formal bug:

Bug BUG-000084665:

Attempting to create an equilateral triangle with lengths of 0.0027 gives, "The feature could not be created. The geometry is invalid." However, the distance between each vertex is greater than the XY Tolerance (0.001) and the XY Resolution (0.0001).

Alternate Solution:

Ensure the distance between each vertex is at least 0.0029 units.

I can't say I am a fan of the Summary text.  I think it downplays the scope and impact of the bug, i.e., someone might think it is about triangles and not the much larger issue of all polygons and lines.

The workaround is basically to not trust the XY Tolerance or XY Resolution of spatial references.  Caveat utilitor.

[DanPatterson_Retired]

wow...works fine with a shapefile...I guess no one else works on projects with millimeter references and sub-meter extents

[JoshuaBixby]

Tough part about forums and not knowing people personally, comments can be read multiple ways....

Anyhow, two other bugs were issued related to the aforementioned Support case:  BUG-000084666 and BUG-000084667.  Even though the impact might seem limited on the surface, cases like this can point to larger issues with managing precision deeper in the software stack.  What isn't stated in the public-facing Support system is that the "bug is not coordinate specific and it is not spatial reference specific."  Although the illustrative example involved working in a coordinate system that presented an unusual real-world case, it has been demonstrated the issue is abstracted from a single coordinate system.

In the end, the true impact may be small but maybe not.  It isn't always easy to tell until developers dig into the bug deeper.

[DanPatterson_Retired]

Which I was surprised find out, that the original errors that I posted on when using geometry objects without a SR hadn't been addressed in several releases and still is an issue, from what I am finding.

If you are working in small areas and need a 0,0 origin...it is a pain that you have to specify a coordinate system in the first place...

When people ask where I do my most work I say ... about 600km south of Accra and 750 km west of Sao Tome