5

What I'd like to do, using ArcPy: use a line to split a bunch of polygons into separate parts. Here is what the inputs could look like:

enter image description here

In this case the line should split each polygon into two separate polygons, creating 6 separate polygons in the end.

This question has been asked many times, and there have been a lot of answers provided, but none of them seem to mention the polygon's cut method, perhaps because it seems to only have been introduced in 10.2.

Here is what the ArcGIS Help shows (scroll to cut (cutter)):

enter image description here

This seems to be exactly what I need, but I cannot entirely figure out how to implement it.

From the few examples I was able to find, it seems like I cannot use the method directly on a shapefile that I am referencing, but I need to access the underlying geometry/vertices, so here is what I did (please correct me if there is an entirely different way of doing it):

  • read shapefile's vertices using a cursor
  • push each vertex into an array
  • use that array to create a polygon/line

The next step would be to to use the cut() method, as shown in the screenshot above, to actually split the polygons, but I cannot get it to work.

What follows is my code (only the polygon-part uses comments as the polyline is just repetition):

import arcpy

polygon = r"D:\Data\poly.shp"
line = r"D:\Data\line.shp"

# ------------------- POLYGON -------------------
# first read the geometry
# i.e.: read all of the vertex coordinates, and push them into an array
dataset_vertices = []

# loop through each feature
for row in arcpy.da.SearchCursor(polygon, ["OID@", "SHAPE@"]):
    print("Here comes feature {0}:".format(row[0]))
    partnum = 0
    # Step through each part of the feature
    for part in row[1]:
        # Print the part number
        print("Part {0}:".format(partnum))

        tempPart = [] # push each vertex into an array
        # Step through each vertex in the feature
        for pnt in part:
            if pnt:
                # Print x,y coordinates of current point
                print("{0}, {1}".format(pnt.X, pnt.Y))
                tempPart.append([pnt.X,pnt.Y])
            else:
                # If pnt is None, this represents an interior ring
                print("Interior Ring:")
        # push part-array into entire array
        dataset_vertices.append(tempPart)
        partnum += 1


# A list that will hold each the Polygon objects
polygon_features = []

for feature in dataset_vertices:
    # Create a Polygon object based on the array of points
    # Append to the list of Polygon objects
    polygon_features.append(
        arcpy.Polygon(
            arcpy.Array([arcpy.Point(*coords) for coords in feature])))

# ------------------- POLYLINE -------------------
line_vertices = []

for row in arcpy.da.SearchCursor(line, ["OID@", "SHAPE@"]):
    print("Here comes feature {0}:".format(row[0]))

    partnum = 0
    for part in row[1]:
        print("Part {0}:".format(partnum))

        tempPart = []
        for pnt in part:
            if pnt:
                print("{0}, {1}".format(pnt.X, pnt.Y))
                tempPart.append([pnt.X,pnt.Y])
            else:
                print("Interior Ring:")
        line_vertices.append(tempPart)
        partnum += 1


line_features = []

for feature in dataset_vertices:
    line_features.append(
        arcpy.Polygon(
            arcpy.Array([arcpy.Point(*coords) for coords in feature])))

After running this code I can use methods, such as polygon_features[0].boundary(), but I am not sure how to use cut().

  • 1
    You are correct that documentation on the cut feature is lacking. Try this implementation below. Note: You also need a polyline feature to cut the polygon feature with. gis.stackexchange.com/questions/124198/… – Tristan Forward Jul 20 '15 at 11:36
  • I tried this code earlier today but it only cuts one polygon, deleting all the others. Or did I miss something? :-) Furthermore I would like to export the result into a separate file, without overriding the original polygon file. – BritishSteel Jul 20 '15 at 12:01
  • Try using the code at the very start of the post. This is the one you need I believe, as it updates the geom each time. If your still running into issues with it deleting the other features have it create a new featureclass and write the new polygons to that one versus back in on itself, as how this script is doing at the moment. The "polygon[1]" at the end of the script is what is causing you problems. – Tristan Forward Jul 20 '15 at 12:28
  • Creating a copy is a very good idea, and an easy way to avoid that issue, good call! I've use the code from the other topic now, but the issue persists. It always splits only one polygon, and deletes any other polygon the line runs through. If I have polygons that are not crossed by the line, then those polygons are not touched. So it always leaves me with the untouched polygons as well as one split polyon. – BritishSteel Jul 20 '15 at 15:00
  • You mentioned the "polygon[1]"; what is it actually doing, or why is it causing the problem? – BritishSteel Jul 20 '15 at 15:01
3

Feature To Polygon will slice polygons by line feature classes when you use the two feature classes as its input.

Inputs: enter image description here

Results: enter image description here

Depending on line geometry this method can create excess polygons. If so, clip your output feature class by your original input polygon feature class to eliminate the new polygons.

If you only want certain polygons sliced by certain lines, use an SQL expression and feature layers to limit your inputs.

  • I have seen this tool in some answers, but never really understood what it was doing. A big thumbs up for the clarification. This solves my issue, but I will not check the answer for now, as this question specifically addresses the cut() method. If no answer has been provided in a few days, I will mark yours as the correct answer, as it definitely solves the problem as a whole. – BritishSteel Jul 21 '15 at 9:26
  • To anyone reading this: please note that this solution requires an ArcInfo/Advanced license. – BritishSteel Jul 21 '15 at 9:26
  • Just got it to work with this method, but it is only ideal if your shapes are separate. Let's say you have a third shape that slightly overlaps with both the right and the left one, then the empty area between all of these three shapes and the line will be filled with a new polygon. Depending on what your inputs look like it will take a lof of tweaking to get it all to work. I ended up using a Select by Location and some other tools to get it to work. – BritishSteel Jul 21 '15 at 13:51
  • 1
    @BritishSteel Yea it looks like the tool does more than slice polygons, depending on the geometry. To remove the excess features, my method would be: 1. select tool output by location: input polygons, are within the source layer feature. 2. Select by attribute: Switch selection 3. (if selection) Delete features – Emil Brundage Jul 21 '15 at 16:11
  • That is exactly what I ended up doing; glad the logic is correct. It seems very overkill to end up getting split polygons, but at least that is a way to get it to work in ArcGIS/ArcPy! – BritishSteel Jul 21 '15 at 20:35

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.