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:
- Create Polygon Shapefile
- Open input file and get coordinates for P
- Loop: Calculate P_Buffer coordinates (Vincenty's formulae) and create an arcpy Point Geometry (100 times)
- Create an arcpy Polygon Geometry from 100 Point Geometries
- 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:
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.
arcpy.env.XYTolerance
).