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
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!
...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'> >>>
Membros conectados podem postar, seguir atualizações e mais. Novo aqui? Registre uma conta gratuita.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.
