3

I'm currently trying to code a geodesic buffer from points (the one that is implemented in ArcGIS is very slow). The idea is to calculate 100 point coordinates (P_Buffer) around an initial point (P) and to shape a polygon from the P_Buffer Points (accordingly, the buffer polygon has 100 vertices).

The script is structured like this:

  1. Create Polygon Shapefile
  2. Open input file and get coordinates for P
  3. Loop: Calculate P_Buffer coordinates (Vincenty's formulae) and create an arcpy Point Geometry (100 times)
  4. Create an arcpy Polygon Geometry from 100 Point Geometries
  5. Write Polygon Geometry (Insert Cursor)

That's the script (simplified):

# Step 1

infile = "CRATER_Pickering.shp"
outfile = "C:\test.shp"
sr = arcpy.Describe(infile).spatialReference
arcpy.CreateFeatureclass_management("C:/", "test.shp", "POLYGON", "", "", "", sr)
cur = arcpy.da.InsertCursor(outfile, ["SHAPE@"])

# Step 2 & 3 work perfectly fine - I skip most part of the section. Each time the loop runs, 
# step 3 ends with adding the P_Buffer coordinates to an arcpy Point Array. After that, 
# the arcpy Point array is appended to an arcpy Array, which stores all the P_Buffer 
# coordinates / the buffer polygon's vertices.

...(get buffer distance and coordinates of P)...

    polygon = arcpy.Array()

        ...(coordinate transformation)... 

        point = arcpy.Point(P_Buffer_X, P_Buffer_Y)
        polygon.append(point)

# Step 4

    polygon_shape = arcpy.Polygon(polygon)

# Step 5

    cur.insertRow([polygon_shape])

Everything works fine until I get to step 5. The coordinates are calculated correctly (just as in the ArcGIS geodesic buffer) and all 100 P_Buffer points are polygon_shape vertices in Step 4. However, when writing the geometry in step 5, only few buffer polygons have 100 vertices. Big circular buffer polygons have 100 vertices but the smaller their diameter gets, the fewer vertices the polygons have (up to zero - meaning that smaller buffer polygons either show angular boundaries or that no geometry is visible).

Here's an example. The green dots in the image show the P_Buffer points from the polygon array in step 3. The buffer polygons written in step 5 are displayed in pinkish colour.

Huge diameters are fine:

Huge diameters are fine

Small diameters are not:

Small diameters are not

How is this possible? When I write a Point Shapefile from the polygon array (Step 3 - all 100 vertices are still there) or count the vertices in polygon_shape (step 4 - all 100 vertices are still there), everything looks fine. As soon as I call insertRow in Step 5, vertices are excluded. Do you have any idea what I did wrong?

Solved:

It was indeed a spatial reference issue. As suggested by klewis, I added the spatial reference to step 4. Now the polygons have all the vertices from the coordinate transformation included.

7
  • 2
    What SRS are you using? I've had similar issues when creating new geometry through arcpy when using WGS84. I'll have to see if I can find exactly what I did to correct the issue but I believe it revolves around the coordinate tolerance (ie. arcpy.env.XYTolerance). Jun 2, 2016 at 12:41
  • 1
    The thing with the coordinate tolerance is a good idea since the arcpy.polygon validates and simplifies the geometry before it creates the feature. Maybe thats why it only happens on small ones.
    – Matte
    Jun 2, 2016 at 12:48
  • 2
    Definitely a spatial reference issue. Not tolerance but precision, though.
    – Vince
    Jun 2, 2016 at 13:01
  • 2
    You should provide sr when creating the Polygon, like this - arcpy.Polygon(polygon, sr)
    – klewis
    Jun 2, 2016 at 13:02
  • Thank you guys. It was indeed a spatial reference issue. As suggested by klewis, I added the spatial reference to step 4. Now the polygons have all the vertices from the coordinate transformation included. To get back to your question, Evil Genious, I am working on Mars (GCS_Mars_2000).
    – C.Riedel
    Jun 2, 2016 at 13:12

1 Answer 1

3

Try adding sr to arcpy.Polygon(polygon, sr)

Here is the reason why:

I sometimes see significant coordinate movement when I do not supply a Spatial Reference to an Arcpy Geometry object such as a Polygon or Polyline. Through research I've learned the default XY Resolution (significant digits) is .0001. This is correct for Projected data but much too large for Geographic data. I try to provide a Spatial Reference whenever I work with Geometry objects, it is critical for Geographic data. I wish the ArcGIS geometry examples focused more on this common pitfall.

1
  • 1
    A tenth of a millimeter is actually usually significant overkill for projected data, causing a significant increase in geometry storage cost for complex lines and polygons for both file and enterprise geodatabases (except when using native storage). The best tool to explain spatial references is Understanding Coordinate Management in the Geodatabase
    – Vince
    Jun 2, 2016 at 23:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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