7

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.

  • 1
    What's the point of this exercise? – Vince Nov 11 '16 at 3:08
  • 2
    No it's not obvious. We don't know what spatial reference you're using in your mxd. Just trying to help. You don't set the SR of buf in the code. – jbalk Nov 11 '16 at 4:47
  • 2
    maybe add some more detail to the question. Like point out the problem you're having, the fact that the SR of the mxd is same as SR of POLYGONS.shp – jbalk Nov 11 '16 at 4:49
  • 1
  • 1
    I think you should report this to ESRI as a bug. – Hornbydd Nov 11 '16 at 16:15
7
+200

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!

  • 2
    Nice catch. Have tried GeoJSON as well? As far as I remember it could specify beziers/arcs separate from multipolygon. – fatih_dur Nov 17 '16 at 23:48
  • 1
    @fatih_dur, great point and that is the winner! – Paul Nov 18 '16 at 0:32
8

This is not an answer but some extra info for this question. The sinister behaviour is not occurring at arcpy.Polygon as suggested by @FelixIP, it is occuring at the getPart() bit of code. I tweaked your code and ran it on a dataset with just a triangle. I had set the coordinate system to be British National Grid when I had created it.

import arcpy
arcpy.env.addOutputsToMap = True
infc = "tri"
outBuffer=r'c:\scratch\OUT_BUFFER.shp'
outPart=r'c:\scratch\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)
            # Print out sequence of XY points
            print(prt)
            for p in prt:
                print str(p.X ) + "," + str(p.Y)
            pgon=arcpy.Polygon(prt,SR)
            arcpy.CopyFeatures_management(pgon, outPart)
            break

Zooming into one corner of the buffer of the triangle in edit mode we can see many vertices:

Vertices of buffer

Infact there are 42 vertices for the whole geometry. My adjusted code reading the output of the getPart() returns only 11.

Clipped part

So the bug is occuring when the getPart() is called to return the Array of Points.

  • 1
    Thank you for confirming this bug +1. You missed fun, at some stage this question had 4 down votes, ha-ha, love it! I am thinking of moving it to featured, perhaps esri people respond, using their site is no pleasure. I'll use avenue to solve this, expect it will work 3 times faster. Almost removed 3.3 when 10.2 became available – FelixIP Nov 11 '16 at 18:31
  • If you look at the WKT, you'll see that all vertices are maintained. buf.WKT.count(",") == 42. arcpy.FromWKT(buf.WKT, SR).area – Paul Nov 11 '16 at 22:03
  • @FelixIP Initially you did not ask a question and you were being arrogant in both question and comments. Probably why it was downvoted – BERA Nov 14 '16 at 17:44

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.