4

Currently building a tool that automates the creation of some polygons that we need done quarterly. Below is just a set script I am using that is giving me a similar issue when I try to run the script. The current output only gives me a polygon with 4 vertices and it is not making a point for each coordinate I have listed. This full script is much larger and complex this is just the one part I am stuck on for some reason. This sample is one of my variations I have tried.

Code

import arcpy

workspace = r"....\Workspace\WritingGeometries"
coordinateSystem = 4326
arcpy.env.workspace = workspace
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = coordinateSystem

fcPointTest = "TestPolygon"
fcPointPath = dsPath + "\\" + fcPointTest
arcpy.CreateFeatureclass_management(dsPath, fcPointTest, "POLYGON")

testCord3 = [arcpy.Point(-95.542806, 30.720361),
             arcpy.Point(-95.479499, 30.805696),
             arcpy.Point(-95.499301, 30.814402),
             arcpy.Point(-95.509831, 30.817513),
             arcpy.Point(-95.516297, 30.818963),
             arcpy.Point(-95.520657, 30.819754),
             arcpy.Point(-95.523791, 30.820232),
             arcpy.Point(-95.526151, 30.820543),
             arcpy.Point(-95.527991, 30.820756),
             arcpy.Point(-95.529465, 30.820909),
             arcpy.Point(-95.530674, 30.821021),
             arcpy.Point(-95.531681, 30.821107),
             arcpy.Point(-95.532535, 30.821174),
             arcpy.Point(-95.533267, 30.821227),
             arcpy.Point(-95.533901, 30.821270),
             arcpy.Point(-95.534457, 30.821305),
             arcpy.Point(-95.542806, 30.821560),
             arcpy.Point(-95.534457, 30.821305),
             arcpy.Point(-95.533901, 30.821270),
             arcpy.Point(-95.533267, 30.821227),
             arcpy.Point(-95.532535, 30.821174),
             arcpy.Point(-95.531681, 30.821107),
             arcpy.Point(-95.530674, 30.821021),
             arcpy.Point(-95.529465, 30.820909),
             arcpy.Point(-95.527991, 30.820756),
             arcpy.Point(-95.526151, 30.820543),
             arcpy.Point(-95.523791, 30.820232),
             arcpy.Point(-95.520657, 30.819754),
             arcpy.Point(-95.516297, 30.818963),
             arcpy.Point(-95.509831, 30.817513),
             arcpy.Point(-95.499301, 30.814402),
             arcpy.Point(-95.606112, 30.805696)]

#Begin Processing
array = arcpy.Array(testCord3)
polygon = arcpy.Polygon(array)
cursor = arcpy.da.InsertCursor(fcPointPath, ["SHAPE@"])
cursor.insertRow([polygon])

del cursor

Below is a screenshot of what my output is.

enter image description here

Solution

Once Paul helped me out and noticed an error in my coordinates I was able to get the script to run through with 0 issues.

enter image description here

  • Not sure why the other vertices are not showing up. Maybe use copy features tool and see what happens – ziggy Oct 10 '16 at 15:59
  • 1
    This is a coordinate reference issue. You should always extract the coordinate reference from the source feature class and apply it at geometry construction, because the default constructor assumes a projected coordinate system and millimeter (0.001) precision – Vince Oct 10 '16 at 16:21
1

Your point locations are what's causing the issue. If you add them in as points, you'll see that 14 of the points are coincident, meaning you are trying to create a polygon part that has 0 area (a line essentially). If you remove all those points, you'll see the result is actually the same.

Don't forget to close your polygons and set the spatial_reference when constructing the Polygon object:

testCord3.append(testCord3[0])
polygon = arcpy.Polygon(array, arcpy.env.outputCoordinateSystem)

Take a look at polygon.getPart(0) and array in your IDE to see what's happening. Nothing wrong with the cursor at all.

  • 1
    Thank you Paul and everyone else who took time to reply to my post. I am happy to say the problem is now resolved. Paul you made the comment that my points are coincident and it never even occurred to me that I should QC that. Come to find out earlier in the script my math had an addition sign where it should have had a subtraction sign to calculate the points on the other side of the sector. Once I fixed that and got my coordinates correct the polygon mapped out as it should. Also great point on making sure I close my polygons, everything is working like a champ. Thanks again everyone! – HiMyNameIzMatt Oct 10 '16 at 17:10
  • @HiMyNameIzMatt, huzzah! Welcome to GIS.SE :) – Paul Oct 10 '16 at 17:11
  • @HiMyNameIzMatt If Paul's answer has solved your problem, please mark it as the answer. See What should I do when someone answers my question? – Midavalo Oct 10 '16 at 18:10
1

I converted your coordinates into a POLYGON in Well-Known Text and got an error (polygon does not close properly), so I added the first point to the end of the point list, and generated this shape: enter image description here

Then I created a LINESTRING in Well-Known Text and from there into a shapefile and got this:

enter image description here

Then I shifted to POINT and got this (draped over the line, with vertex order labeled):

enter image description here

So it appears that you have a zero-area extension between the third and final vertices, which ArcGIS has "helpfully" removed.

Guidelines for constructing geometry in ArcPy:

  • Always extract the SpatialReference from a feature class (this is not the same as a coordinate system) before constructing geometries:

    sr = arcpy.Describe(tabPath).spatialReference

    You can use the SpatialReference used to create the feature class, but then it's wise to specify the origin and units which correspond to your data (in this case, with 6 digits precision, one millionth of a degree)

    sr = arcpy.SpatialReference(4326)

    sr.setFalseOriginAndUnits(-400.0,-400.0,1000000.0)

    but be forewarned that in certain cases, the created feature class may have different parameters than what was requested, so extracting from a Describe is safer.

  • Always apply the coordinate reference in the geometry constructor

    polygon = arcpy.Polygon(array,sr)

  • Always include the closing point in polygons (yes, ArcGIS, will do it for you, but this may mask issues)

         arcpy.Point(-95.499301, 30.814402),
         arcpy.Point(-95.606112, 30.805696),
         arcpy.Point(-95.542806, 30.720361)]
    
  • If something goes wrong, try plotting your coordinates, first as lines, then as points, if necessary, to review correctness.

  • Many thanks for the reply, a bad list of coordinates ended up being the main contributor for the problem. I was neglectful and didn't QC the list I was testing from the beginning and once I revised the list my issue was resolved. Thanks again! – HiMyNameIzMatt Oct 10 '16 at 18:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.