Friday, 13 November 2015

arcpy - getPart() method returns incorrect geometry from buffer in ArcGIS


I created 4 most basic polygons:


Sample data


And ran this script:


import arcpy
infc = r'D:\Scratch\POLYGONS.shp'
d=arcpy.Describe(infc)
SR=d.spatialReference

with arcpy.da.SearchCursor(infc, ("SHAPE@","LABEL")) as rows:
for shp,label in rows:
buf=shp.buffer(5)
n=buf.partCount
for i in xrange (n):
prt=buf.getPart(i)
pgon=arcpy.Polygon(prt,SR)
arcpy.AddMessage((int(buf.area),int(pgon.area),label))

OUTPUT:



[68735, 68707, u'RECTANGLE']
[29892, 29846, u'TRIANGLE']
[47832, 47832, u'CIRCLE']
[55126, 55060, u"'COMPLEX'"]

Anybody could possibly help me to understand result.


UPDATE:


I modified the code:


import arcpy
infc = r'....POLYGONS.shp'

outBuffer=r'...OUT_BUFFER.shp'
outPart=r'...SCRARCH\OUT_PART.shp'
d=arcpy.Describe(infc)
SR=d.spatialReference
with arcpy.da.SearchCursor(infc, ("SHAPE@","LABEL")) as rows:
for shp,label in rows:
buf=shp.buffer(5)
arcpy.CopyFeatures_management(buf, outBuffer)
n=buf.partCount
for i in xrange (n):

prt=buf.getPart(i)
pgon=arcpy.Polygon(prt,SR)
arcpy.CopyFeatures_management(pgon, outPart)
break

to check actual shapes. When working with geometries it would be expected that the first part of a single part polygon would be identical to the polygon itself.


This is an overlay of input and outputs and you can see where mismatch happened - at the corners:


enter image description here


Something very sinister is happenning inside arcpy.Polygon, it simplifies the shape using some sort of 'logic'. Handling true curves(?). Why result is so ugly?


NOTE:




  • same result obtained with exporting to geodatabase

  • both outputs came up with input's spatial reference


This buffer issue is a part of a bigger job, where I have to deal with potentially multi-part polygons. The workaround I am using is explosion of buffer to single parts, using arcgis tool (!). It works, but ridiculously slow due to number of iterations. Why on earth the abc of geometry handling is not working?


If you know solution different from my workaround, I'll award 200 points to your answer.


Well, as it happens, after posting I realised I shouldn't be bothered with parts of the buffer, because label point of polygon is placed inside larger part of multipart polygon.



Answer



Instead of using getPart(), you can use WKT on the buffered polygon itself. It's not the exact shape of the buffer (due to approximating the 3 true curves on the triangle with 15 vertices each), but it's quite close. From the help:




Any true curves in the geometry will be densified into approximate curves in the WKT string.



import json
from pandas import DataFrame
import arcpy


header = ("type", "area", "vertices")
sr = arcpy.SpatialReference(3857) #WMAS
coords = [(-10136090, 3460507),

(-10135313, 3461773),
(-10134909, 3460488),
(-10136090, 3460507)]
poly = arcpy.Polygon(arcpy.Array(arcpy.Point(x,y) for x,y in coords), sr)

buffer = poly.buffer(5)
part = arcpy.Polygon(buffer.getPart(0), sr)
wkt = arcpy.FromWKT(buffer.WKT, sr)
wkb = arcpy.FromWKB(buffer.WKB)
densify = arcpy.Polygon(buffer.densify("DISTANCE", 1, .00001).getPart(0), sr)

esriJSON = arcpy.AsShape(json.loads(buffer.JSON), True)

print(DataFrame([("buffer", buffer.area, buffer.pointCount)], columns=header))
df = DataFrame([("getPart", part.area, part.pointCount),
("WKT", wkt.area, wkt.pointCount),
("WKB", wkb.area, wkb.pointCount),
("densify", densify.area, densify.pointCount),
("JSON", esriJSON.area, esriJSON.pointCount)],
columns=header)
df["delta"] = [v - buffer.area for v in df["area"]]


print(df.sort(columns='delta',ascending=False))

And here's the output:


     type            area  vertices
0 buffer 775085.643363 9
type area vertices delta
4 JSON 775085.643363 9 1.641456e-08
3 densify 775085.633439 4235 -9.923947e-03
1 WKT 775085.297882 43 -3.454809e-01

0 getPart 775041.921399 9 -4.372196e+01
2 WKB 775041.921399 9 -4.372196e+01

Based on comment below from @faith_dur, the way to go is esri JSON!


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...