topic Bisect Polygon w/ Arcpy in Python Questions

Bisect Polygon w/ Arcpy

NoahHuntington — Wed, 02 Jul 2014 18:59:31 GMT

[ATTACH=CONFIG]35043[/ATTACH]

How might I bisect the attached polygon in a NE/SW direction using arcpy? Looking for a line feature class as output also. Any suggestions?

Re: Bisect Polygon w/ Arcpy

T__WayneWhitley — Tue, 08 Jul 2014 22:46:29 GMT

...first post on GeoNet, thought I'd use this thread (hope you don't mind Noah) to try out the experience on the new forum format. Since this seems to be a one-off, I'll even try posting some really rough code below that was constructed interactively (IDLE), but 1st I'll explain the so-called solution. This may not be good, but it's what came to mind - here it is in basically 3 steps:

1- In the Features toolset (Data Management toolbox), use Minimum Bounding Geometry with the RECTANGLE_BY_AREA geometry type (basic licence level)...

(This part is what is actually shown in the code section below, working on the Min Bound Geom result.)

2- With the known 4 corners of the 'bounding polygon' output, extract the midpoints along the major axis...

3- use some geometry methods and some kind of overlay processing to 'split' the orig poly.

I apologize firsthand for what I resorted to doing in the final step 3, kind of lost the handle on what I was doing and cheated a little using an 'advanced' (and probably inefficient for this use case) couple of heavyduty Feature To Polygon and Feature To Line tool executions (since you wanted the 'oval' outline in 2 line 'mirrored' parts, if I understood correctly?)...yeah, I know, failed to do it all with the lower license and I know well there's a workaround in order to complete the task, but frankly it's beer-thirty and I'm thirsty. If it cripples you at the last step I can help out tomorrow, but it does the trick with just ver 10.0 (easily adapted to more advanced techniques via 10.2.x). Okay, remember this was only meant for a one-off, no polish whatsoever, and only meant as a demo to help furnish you an idea how to go about this problem. Of course this could be dramatically shortened, and improved esp where I failed to include the final hack converting the geom without requiring a higher-than-Basic Desktop license...maybe someone else will do that. Enjoy, Wayne

...so let's see if I can figure out how to paste this rough-hacked solution:

>>> rows = arcpy.SearchCursor("the MinimumBoundingGeometry result on your oval shape")
>>> row = rows.next()
>>> feat = row.Shape


>>> for point in feat.getPart(0):
...     print point.X, point.Y
...    
535920.293977   123231.784748
537048.195264   124114.261001
537569.582345   123447.87093
536441.681058   122565.395005
535920.293977   123231.784748


>>> anArray = feat.getPart(0)
>>> for i in range(4):
...     print anArray.getObject(i)
...    
535920.293976922 123231.784747839 NaN NaN
537048.195264421 124114.261001259 NaN NaN
537569.582345005 123447.870929845 NaN NaN
536441.681057505 122565.395004511 NaN NaN


>>> newArray = arcpy.Array()
>>>
>>> # coords of bisect line


>>> pt1X = anArray.getObject(0).X + (anArray.getObject(3).X - anArray.getObject(0).X)/2
>>> print pt1X
536180.987517


>>> pt1Y = anArray.getObject(3).Y + (anArray.getObject(0).Y - anArray.getObject(3).Y)/2
>>> print pt1Y
122898.589876
>>>
>>> # now pt2 to form the other line endpoint...
>>> pt2X = anArray.getObject(1).X + (anArray.getObject(2).X - anArray.getObject(1).X)/2
>>> print pt2X
537308.888805


>>> pt2Y = anArray.getObject(2).Y + (anArray.getObject(1).Y - anArray.getObject(2).Y)/2
>>> print pt2Y
123781.065966
>>>
>>> # create the points...
>>> pt1 = arcpy.Point(pt1X,pt1Y)
>>> pt2 = arcpy.Point(pt2X,pt2Y)
>>>
>>> # add the points to the new array...
>>> newArray.add(pt1);newArray.add(pt2)
>>>
>>> # create the line geometry...
>>> bisect = arcpy.Polyline(newArray)


>>> arcpy.env.workspace = r'C:\Users\whitley-wayne\Documents\ArcGIS\Default.gdb'
>>> arcpy.env.overwriteOutput = True


>>> arcpy.CopyFeatures_management(bisect, 'bisectingLine')
<Result 'C:\\Users\\whitley-wayne\\Documents\\ArcGIS\\Default.gdb\\bisectingLine'>
>>>
>>> # define in_features, out_feature_class
>>> in_features = ['ovalTendingNEtoSW', 'bisectingLine']
>>> out_feature_class = 'bisectedPoly'
>>>
>>> # try the union again with geom copied to line fc:
>>> arcpy.Union_analysis (in_features, out_feature_class)
Runtime error <class 'arcgisscripting.ExecuteError'>: ERROR 000366: Invalid geometry type


>>> # wrong gp tool with ERROR 000366, try feat to poly...
>>> # FeatureToPolygon_management (in_features, out_feature_class, {cluster_tolerance}, {attributes}, {label_features})


>>> arcpy.FeatureToPolygon_management (in_features, out_feature_class)
<Result 'C:\\Users\\whitley-wayne\\Documents\\ArcGIS\\Default.gdb\\bisectedPoly'>
>>>
>>> # didn't work; try increasing the cluster tolerance...
>>> arcpy.FeatureToPolygon_management (in_features, out_feature_class, '1')
<Result 'C:\\Users\\whitley-wayne\\Documents\\ArcGIS\\Default.gdb\\bisectedPoly'>
>>>
>>> # that worked...
>>> arcpy.FeatureToLine_management("bisectedPoly","bisectedPolyToLine")
<Result 'C:\\Users\\whitley-wayne\\Documents\\ArcGIS\\Default.gdb\\bisectedPolyToLine'>
>>>

Re: Bisect Polygon w/ Arcpy

T__WayneWhitley — Wed, 09 Jul 2014 22:59:34 GMT

I got more time to look at this problem again - 3 things I really did not like about the last solution I posted are:

- This was first tried as a 10.0 solution, so I didn't have access to the 10.2.x geometry methods, boundary() and cut(cutter).

- I couldn't think of a 'lower license' solution to 'reconstruct' the oval polygon into 2 bisected parts and resorted to using Feature To Polygon.

- Unfortunately the bisecting line geometry constructed didn't work until I 'fiddled' with the cluster tolerance - this introduced some error.

I have another machine with 10.2.1 installed that yields a better solution (lower license) without introducing error. Just being lazy before, but a simple solution to overcoming the cluster tolerance challenge is to use a little algebra on the bisecting line to extend it enough to effectively use it as a 'cutter'. This bit of code (within IDLE) then 'corrects' the 'bisect' line geometry:

>>> # recall (from the previous post) the point coordinates derived from the MBG output
>>> print pt1X;pt1Y;pt2X;pt2Y
536180.987517
122898.58987617493
537308.88880471326
123781.06596555188
>>> 
>>> slope = (pt2Y - pt1Y)/(pt2X - pt1X)
>>> print slope
0.782405427813
>>> 
>>> # Y2 - Y1 = slope(X2 - X1)
>>> # Then, say, delta X is 1 unit...then Y2 - Y1 = slope.
>>> 
>>> pt1Xext = pt1X -1.0
>>> pt2Xext = pt2X +1.0
>>> print pt1Xext;pt2Xext
536179.987517
537309.88880471326
>>> 
>>> # Y2 = slope + Y1
>>> pt1Yext = -slope + pt1Y
>>> print pt1Yext
122897.807471
>>> 
>>> pt2Yext = pt2Y + slope
>>> print pt2Yext
123781.848371
>>> 
>>> # now load the 'corrected' coords into point objs
>>> pt1 = arcpy.Point(pt1Xext,pt1Yext)
>>> pt2 = arcpy.Point(pt2Xext,pt2Yext)
>>> newArray = arcpy.Array()
>>> newArray.add(pt1);newArray.add(pt2)
>>> 
>>> # from the array, create the line geom
>>> bisect = arcpy.Polyline(newArray)

...here's where it gets interesting - I access the original oval shape to be cut, but unfortunately I run into a spatial reference required on 'bisect' is missing when trying to make the 'cut' (or bisect):

>>> # same 10.0 command (old-style cursor)
>>> rows2 = arcpy.SearchCursor('ovalTendingNEtoSW')
>>> row2 = rows2.next()
>>> feat2 = row2.shape
>>> type(feat2)
<class 'arcpy.arcobjects.geometries.Polygon'>
>>> 
>>> # now try the 'new' geom methods
>>> feat2bound = feat2.boundary()
>>> 
>>> # 'bisect' is the 'cutter':
>>> # cut(cutter)
>>> theBisected = feat2bound.cut(bisect)


Traceback (most recent call last):
  File "<pyshell#43>", line 1, in <module>
    theBisected = feat2bound.cut(bisect)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\arcobjects\arcobjects.py", line 825, in cut
    return convertArcObjectToPythonObject(self._arc_object.Cut(*gp_fixargs((other,))))
RuntimeError: All geometries involved in this operation must have the same spatial reference.

So, this was corrected...the bisect was accomplished - without the need for a higher license (and without the tolerance error).

>>> feat2.spatialReference.PCSName
u'NAD_1983_HARN_StatePlane_Florida_East_FIPS_0901_Feet'


>>> bisect.spatialReference.PCSName
u''


>>> feat2.spatialReference.PCSCode
2881
>>> feat2.spatialReference.factoryCode
2881


>>> sr = arcpy.SpatialReference(2881)
>>> sr.PCSName
u'NAD_1983_HARN_StatePlane_Florida_East_FIPS_0901_Feet'
>>> 
>>> # build 'bisect' with a spatial reference
>>> bisect = arcpy.Polyline(newArray, sr)
>>> bisect.spatialReference.PCSName
u'NAD_1983_HARN_StatePlane_Florida_East_FIPS_0901_Feet'
>>> 
>>> # no fuss now (hopefully)
>>> theBisected = feat2.cut(bisect)
>>> 
>>> type(theBisected)
<type 'list'>


>>> for eachGeom in theBisected:
  print str(type(eachGeom))

<class 'arcpy.arcobjects.geometries.Polygon'>
<class 'arcpy.arcobjects.geometries.Polygon'>


>>> theBisectedLeft = theBisected[0].boundary()
>>> theBisectedRight = theBisected[1].boundary()
>>> 
>>> # copy the results to gdb
>>> arcpy.Merge_management([theBisectedLeft,theBisectedRight],'bisectedOvalBoundary')
<Result '\\\\MC-GISAPPSERVER\\mapdocs\\temp\\testData.gdb\\bisectedOvalBoundary'>
>>> 
>>> # not quite yet, but close...
>>> # need the bisection of the orig boundary poly; not the boundaries of the bisected polys
>>> theBisected = feat2bound.cut(bisect)


>>> arcpy.Merge_management([theBisected[0],theBisected[1]],'bisectedOvalBoundary2')
<Result '\\\\MC-GISAPPSERVER\\mapdocs\\temp\\testData.gdb\\bisectedOvalBoundary2'>
>>> # ah, beautiful!

Re: Bisect Polygon w/ Arcpy

NoahHuntington — Fri, 18 Jul 2014 13:27:16 GMT

Hey Wayne.

Thanks for choosing my example for your experiment. This is a very ingenious solution and a thorough explanation! I will have to spend a little time working thru the details but I am very grateful for this excellent response. I like you am still familiarizing myself with Geonet. But I like what I have seen so far... Thanks a million, I will add to when I have a grasp on this approach. So appreciated. Thank you!!

-Noah