Memisahkan poligon dengan garis menggunakan ArcPy?

Dalam melihat Buffer dengan penghalang fisik menggunakan ArcGIS untuk Desktop? , terpikir oleh saya bahwa saya tidak yakin bagaimana seseorang akan menggunakan alat geoprocessing di ArcGIS untuk membagi poligon dengan garis secara terprogram.

Secara manual, Anda akan menggunakan alat Cut Polygons atau alat Split Polygons pada toolbar Topology , tetapi bagaimana Anda menyelesaikan tugas yang sama menggunakan modelbuilder atau python groprocessing scripting tools?

Langsung dari kelelawar saya memikirkan semua alat di toobox Analisis seperti Union, Identity, dll, tetapi semua itu adalah alat Polygon-Polygon, BUKAN alat Polygon-Line. Bahkan alat Split adalah Polygon-Polygon.

Ada ide?

Menggunakan ET Geowizard Anda dapat mengakses kode untuk alat Split Polygons with Polylines:

Berikut tautan ke skrip .

Atau, Anda dapat menggunakan ArcObjects untuk melakukan ini:

Potong Cuplikan Polygon

Anda juga dapat menggunakan metode buffer satu sisi yang dijelaskan di sini .

Setelah faktanya, saya akhirnya membuat alat ModelBuilder saya sendiri. Saya lupa tentang pertanyaan ini dan memposting solusi saya untuk pertanyaan serupa lainnya . Untuk kelengkapan, ini adalah repost dari jawabannya:

Saya pikir pasti ada cara untuk melakukan ini, jadi saya menciptakan apa yang saya yakini sebagai solusi yang cukup bagus. Saya telah mempostingnya di situs ArcGIS Resources di Community-> Technical-> Analysis & Geoprocessing-> Analysis-> Gallery.

Alat ini disebut Split Polygons With Lines dan memerlukan lisensi ArcInfo karena beberapa alat yang digunakan dalam model. Pada dasarnya apa yang saya lakukan adalah membuat kotak batas minimum untuk poligon dan memperluas garis ke mereka. Jadi menggunakan beberapa voodoo ModelBuilder, saya bisa mengubah garis menjadi poligon, yang kemudian saya gunakan Identity untuk membagi polys asli.

Silakan mengujinya dan lihat apakah itu cocok untuk Anda. Dalam tes saya (terbatas) itu mempertahankan atribut dari poligon asli, dan hanya membelah poligon yang ada.

Jika Anda ingin keluar dari ArcGIS maka gunakan geom.splitpolysbylines .

Secara pribadi saya belum pernah menggunakannya dalam suatu program tetapi saya pikir Anda dapat mengakses commondline ini dengan python, silakan lihat bantuan untuk rincian lebih lanjut.

jika Anda tidak memiliki beberapa masalah akurasi tinggi, Anda akan buffer garis dengan jarak minimum untuk misalnya (0,002 saya pikir itu harus lebih tinggi daripada akurasi kelas fitur Anda), kemudian menerapkan alat hapus ke poligon oleh garis buffered.

Kode arcpy yang diperbarui untuk membagi poligon dalam arah horizontal atau vertikal menggunakan nilai persen

'''-------------------------------------------------------------------------------------------
 Tool Name:   Polygon Bisector in percent parts
 Source Name: polygonbisector.py and splitPolygonsWithLines
 Version:     ArcGIS 10.0
 Author:      ESRI, Inc.
 Required Arguments:
              Input Features (Feature Layer)
              Output Feature Class (Feature Class)
 Optional Arguments:
              Axis (X|Y)
              Group Field(s) (Field)
              Acceptable Error Percent (Double)

 Description: Computes a line that bisects, or divides in Commulative Percent value, a polygon area along a line
              of constant latitude or longitude. each successive input polygon's area will be
              on either side of the bisecting line.
https://www.arcgis.com/home/item.html?id=9aadb577ccb74f0e88b13a0e3643ca4d
Credits (Attribution)
Esri, Inc. [email protected]

http://www.arcgis.com/home/item.html?id=cd6b2d45df654245b7806a896670a431
Split Polygons using Line features
Shapefile by daltongis
Credits (Attribution)
GIS.StackExchange.com

updated by : Sham Davande, GIS Analyst - [email protected]
-------------------------------------------------------------------------------------------'''

# Import system modules
import arcpy
from arcpy import env
import os
import sys
import math
import time
arcpy.env.overwriteOutput = True


# Change input polygon feature class name and path here
split_poly = "D:\\temp\\polygon2.shp"

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
#   Split percent   Commulative Percent value for python code
#   3.7         3.7
#   25.5            29.2
#   20.8            50.0
#   10.5            60.5
#   31.5            92.0
#   8.0 

# its cumulative percent values to split polygon with N number percent parts
list = ["3.7","29.2","50.0","60.5","92.0"]

# Set local variables
out_folder_path = "C:\\" 
out_name = "temp71"
# Set workspace
env.workspace = "C:\\temp71"

# Execute CreateFolder
arcpy.CreateFolder_management(out_folder_path, out_name)

# Set local variables
out_folder = "C:\\temp71" 
out_name = "NewGdb.gdb"
out_name2 = "C:\\temp71\\NewGdb.gdb"
geometry_type = "POLYLINE"
template = ""
has_m = "DISABLED"
has_z = "DISABLED"

# Execute CreateFileGDB
arcpy.CreateFileGDB_management(out_folder, out_name)

percent_lines = "C:\\temp71\\NewGdb.gdb\\line1"
line2 = "C:\\temp71\\NewGdb.gdb\\line2"
line3 = "C:\\temp71\\NewGdb.gdb\\line3"
allLines = "C:\\temp71\\NewGdb.gdb\\all_lines"
outFeatureClass1 = "polygon"
outFeatureClass2 = "C:\\temp71\\NewGdb.gdb\\polygon"
expression = ""

#############################################
# No hard coded input files path mentioned in the code below here after
# and don't want interfere somebodies C:\Temp folder, so C:\temp71 folder automatically created by code below. 
#############################################


file_prj1 = os.path.splitext(split_poly)[0]
file_prj2 = str(file_prj1 + ".prj")
out_coordinate_system = str(file_prj2)

# Execute CreateFeatureclass
arcpy.CreateFeatureclass_management(out_name2,"all_lines", "POLYLINE", "", "DISABLED", "DISABLED", out_coordinate_system)
# Execute FeatureClassToFeatureClass
arcpy.FeatureClassToFeatureClass_conversion(split_poly, out_name2,
                                            outFeatureClass1, expression)

for i in list:
    if i > 0:
        percent_area = float(i)

        ith_part = 100/percent_area

        print "Generating a polyline to split a polygon into two horizontal parts of", percent_area, "% and",100-percent_area, "% areas"
        print ith_part, "ith part of polygon"

        schemaType = "NO_TEST"
        fieldMappings = ""
        subtype = ""

        # Main function, all functions run in GravityModel
        def PolygonBisector(in_features, out_fc, axis="x", groupfields=[], error=0.001):
        # Error if sufficient license is not available
            if arcpy.ProductInfo().lower() not in ['arcinfo']:
                arcpy.AddError("An ArcInfo/Advanced license is required.")
                sys.exit()

            # Set geoprocessing environments
            arcpy.env.overwriteOutput = True
            arcpy.env.qualifiedFieldNames = False

            shapefield = arcpy.Describe(in_features).shapeFieldName
            rounder = GetRounder(in_features)

        # If group fields are specified, dissolve by them
            if groupfields:
                in_features = arcpy.management.Dissolve(in_features, "in_memory/grouped", groupfields)
            else:
                groupfields = [arcpy.Describe(in_features).OIDFieldName]
            fields = [shapefield] + groupfields

        # Create output feature class and set up cursor
            icur = irow = scur = None
            arcpy.management.CreateFeatureclass(os.path.dirname(out_fc), os.path.basename(out_fc), "POLYLINE", "", "", "", arcpy.Describe(in_features).spatialReference)
            arcpy.management.AddField(out_fc, "Group_", "TEXT", "", "", "", "Group: {0}".format(", ".join(groupfields)))
            icur = arcpy.InsertCursor(out_fc)
            scur = arcpy.SearchCursor(in_features, "", "", ";".join(fields))
            count = int(arcpy.management.GetCount(in_features).getOutput(0))
            arcpy.SetProgressor("step", "Processing polygons...", 0, count, 1)
            bigi = 1

            # Begin processing
            try:
                for row in scur:
                    minx = miny = float("inf")
                    maxx = maxy = float("-inf")
                    totalarea = 0
                    feat = row.getValue(shapefield)
                    totalarea = row.getValue(shapefield).area
                    group = []
                    for field in groupfields:
                        group.append(str(row.getValue(field)))
                    partnum = 0
                    # Get the min and max X and Y
                    for part in feat:
                        for point in feat.getPart(partnum):
                            if point:
                                minx = point.X if point.X < minx else minx
                                miny = point.Y if point.Y < miny else miny
                                maxx = point.X if point.X > maxx else maxx
                                maxy = point.Y if point.Y > maxy else maxy
                        partnum += 1

                    # Process the polygon
                    # Some variables
                    conditionmet = False
                    difference = 0
                    lastdifference = float("inf")
                    differences = {}
                    itys = {}
                    i = 1
                    strike = 0

                    # The starting bisector (half the distance from min to max)
                    if axis == "x":
                        ity = (miny + maxy)/2.0
                    else:
                        ity = (minx + maxx)/2.0

                    while not conditionmet:
                        # Construct a line through the middle
                        if axis == "x":
                            line = MakeBisector(minx, maxx, ity, in_features, axis)
                        else:
                            line = MakeBisector(miny, maxy, ity, in_features, axis)
                        # The FeatureToPolygon function does not except a geometry object, so make a temporary feature class
                        templine = arcpy.management.CopyFeatures(line, "in_memory/templine")
                        temppoly = arcpy.management.CopyFeatures(feat, "in_memory/temppoly")
                        # Intersect then Feature To Polygon
                        bisected = arcpy.management.FeatureToPolygon([temppoly, templine], "in_memory/bisected")
                        clip = arcpy.analysis.Clip(bisected, in_features, "in_memory/clip")

                        # Group bisected polygons according to above or below the bisector
                        arcpy.management.AddField(clip, "FLAG", "SHORT")
                        ucur = arcpy.UpdateCursor(clip, "", "")
                        flag = 0
                        try:
                            for urow in ucur:
                                ufeat = urow.getValue(arcpy.Describe(clip).shapeFieldName)
                                partnum = 0
                                for upart in ufeat:
                                    for upoint in ufeat.getPart(partnum):
                                        if upoint:
                                            if axis == "x":
                                                if round(upoint.Y, rounder) > round(ity, rounder):
                                                    flag = 1
                                                    break
                                                elif round(upoint.Y, rounder) < round(ity, rounder):
                                                    flag = -1
                                                    break
                                            else:
                                                if round(upoint.X, rounder) > round(ity, rounder):
                                                    flag = 1
                                                    break
                                                elif round(upoint.X, rounder) < round(ity, rounder):
                                                    flag = -1
                                                    break
                                    partnum += 1
                                urow.setValue("FLAG", flag)
                                ucur.updateRow(urow)
                        except:
                            raise
                        finally:
                            if ucur:
                                del ucur

                        # Check if the areas are halved
                        dissolve = arcpy.management.Dissolve(clip, "in_memory/dissolve", "FLAG")
                        scur2 = arcpy.SearchCursor(dissolve)
                        try:
                            for row2 in scur2:
                                firstarea = row2.getValue(arcpy.Describe(dissolve).shapeFieldName).area
                                firstflag = row2.getValue("FLAG")
                                break
                        except:
                            raise
                        finally:
                            if scur2:
                                del scur2

                        difference = abs(firstarea - (totalarea/ith_part))
                        differences[i] = difference
                        itys[i] = ity
                        print round(100*(difference/(totalarea/ith_part)),5)
                        #arcpy.AddWarning(round(100*(difference/(totalarea/ith_part)),5))
                        # Stop if tolerance is achieved
                        if (difference/(totalarea/ith_part))*100 <= error:
                            conditionmet = True
                            break
                        # Moving the line in the wrong direction? due to coordinate system origins or over-compensation
                        if difference > lastdifference:
                            firstflag = firstflag*-1.0
                        # If we're not improving
                        if abs(difference) > min(differences.values()):
                            strike+=1
                        # Or if the same values keep appearing
                        if differences.values().count(difference) > 3 or strike >=3:
                            arcpy.AddWarning("Tolerance could not be achieved. Output will be the closest possible.")
                            # Reconstruct the best line
                            if axis == "x":
                                line = MakeBisector(minx, maxx, itys[min(differences,key = lambda a: differences.get(a))], in_features, axis)
                            else:
                                line = MakeBisector(miny, maxy, itys[min(differences,key = lambda a: differences.get(a))], in_features, axis)
                            break
                        # Otherwise move the bisector so that the areas will be more evenly split
                        else:
                            if firstflag == 1:
                                if axis == "x":
                                    ity = ((ity-miny)/((totalarea/ith_part)/firstarea)) + miny
                                else:
                                    ity = ((ity-minx)/((totalarea/ith_part)/firstarea)) + minx
                            elif firstflag == -1:
                                if axis == "x":
                                    ity = ((ity-miny)*math.sqrt((totalarea/ith_part)/firstarea)) + miny
                                else:
                                    ity = ((ity-minx)*math.sqrt((totalarea/ith_part)/firstarea)) + minx
                            lastdifference = difference
                        i +=1
                    irow = icur.newRow()
                    irow.setValue(arcpy.Describe(out_fc).shapeFieldName, line)
                    irow.setValue("Group_", ", ".join(group))
                    icur.insertRow(irow)
                    arcpy.SetProgressorPosition()
                    arcpy.AddMessage("{0}/{1}".format(bigi, count))
                    bigi +=1

            except:
                if arcpy.Exists(out_fc):
                    arcpy.management.Delete(out_fc)
                raise
            finally:
                if scur:
                    del scur
                    if icur:
                        del icur
                if irow:
                    del irow
                for data in ["in_memory/grouped", temppoly, templine, clip, bisected, dissolve]:
                    if data:
                        try:
                            arcpy.management.Delete(data)
                        except:
                            ""

        def MakeBisector(min,max,constant, templatefc, axis):
            if axis == "x":
                array = arcpy.Array()
                array.add(arcpy.Point(min, constant))
                array.add(arcpy.Point(max, constant))
            else:
                array = arcpy.Array()
                array.add(arcpy.Point(constant, min))
                array.add(arcpy.Point(constant, max))
            line = arcpy.Polyline(array, arcpy.Describe(templatefc).spatialReference)
            return line

        def GetRounder(in_features):
            try:
                unit = arcpy.Describe(in_features).spatialReference.linearUnitName.lower()
            except:
                unit = "dd"
            if unit.find("foot") > -1:
                rounder = 1
            elif unit.find("kilo") > -1:
                rounder = 3
            elif unit.find("meter") > -1:
                rounder = 1
            elif unit.find("mile") > -1:
                rounder = 3
            elif unit.find("dd") > -1:
                rounder = 5
            else:
                rounder = 3
            return rounder

        # Run the script
        if __name__ == '__main__':
            # Get Parameters
            in_features = str(outFeatureClass2)
            out_fc = percent_lines
            axis = arcpy.GetParameterAsText(2).lower() or "x"
            groupfields = arcpy.GetParameterAsText(3).split(";") if arcpy.GetParameterAsText(3) else []
            error = float(arcpy.GetParameter(4)) if arcpy.GetParameter(4) else 0.001
            out_data = line2
            # Run the main script
            PolygonBisector(in_features, out_fc, axis, groupfields, error)
            arcpy.Copy_management(out_fc, out_data)
            # Use Append tool to move features into single dataset
            arcpy.Append_management(out_data, allLines, schemaType, fieldMappings, subtype)
            print "again generating a polyline to split same polygon"
            print ""

        else:
            print "Error"

print "finished"
print ""
print "Splitting polygon using multiple lines generated for a percent area values provided by you..."
def splitPolygonsWithLines(Poly, Lines, LinesQuery="", outPoly=""):
    inputPoly=Poly
    inputLines=Lines
    query=LinesQuery

    inputPolyName=os.path.basename(inputPoly)
    inputLinesName=os.path.basename(inputLines)

    parDir=os.path.abspath(inputPoly+"\..")
    if outPoly=="":
        outputPolyName=os.path.splitext(inputPolyName)[0]+u"_Split"+os.path.splitext(inputPolyName)[1]
        outputPoly=os.path.join(parDir,outputPolyName)
    else:
        outputPolyName=os.path.basename(outPoly)
        outputPoly=outPoly

    sr=arcpy.Describe(inputPoly).spatialReference
    fieldNameIgnore=["SHAPE_Area", "SHAPE_Length"]
    fieldTypeIgnore=["OID", "Geometry"]

    #############################################################################################################################
    arcpy.CreateFeatureclass_management (parDir, outputPolyName, "POLYGON", "", "", "", sr)

    arcpy.AddField_management(outputPoly, "OLD_OID", "LONG")
    for field in arcpy.ListFields(inputPoly):
        if (field.type not in fieldTypeIgnore and field.name not in fieldNameIgnore):
            arcpy.AddField_management (outputPoly, field.name, field.type)


    arcpy.MakeFeatureLayer_management(inputLines,inputLinesName+"Layer",query)
    arcpy.MakeFeatureLayer_management(inputPoly,inputPolyName+"Layer")

    arcpy.SelectLayerByLocation_management(inputPolyName+"Layer","INTERSECT",inputLinesName+"Layer","","NEW_SELECTION")
    arcpy.SelectLayerByAttribute_management(inputPolyName+"Layer", "SWITCH_SELECTION")

    fieldmappings = arcpy.FieldMappings()
    for field in arcpy.ListFields(inputPoly):
        if (field.type not in fieldTypeIgnore and field.name not in fieldNameIgnore):
            fm=arcpy.FieldMap()
            fm.addInputField(outputPoly, field.name)
            fm.addInputField(inputPolyName+"Layer", field.name)

            fm_name = fm.outputField
            fm_name.name = field.name
            fm.outputField = fm_name

            fieldmappings.addFieldMap (fm)

    fm=arcpy.FieldMap()
    fm.addInputField(outputPoly, "OLD_OID")
    fm.addInputField(inputPolyName+"Layer", "OBJECTID")

    fm_name = fm.outputField
    fm_name.name = "OLD_OID"
    fm.outputField = fm_name

    fieldmappings.addFieldMap (fm)

    arcpy.Append_management(inputPolyName+"Layer", outputPoly, "NO_TEST", fieldmappings)

    polySelect=arcpy.SelectLayerByLocation_management(inputPolyName+"Layer","INTERSECT",inputLinesName+"Layer","","NEW_SELECTION")
    lineSelect=arcpy.SelectLayerByLocation_management(inputLinesName+"Layer","INTERSECT",inputPolyName+"Layer","","NEW_SELECTION")
    #############################################################################################################################

    fields=[f.name for f in arcpy.ListFields(inputPoly) if (f.type not in fieldTypeIgnore and f.name not in fieldNameIgnore)]
    fields.append("SHAPE@")
    totalFeatures=int(arcpy.GetCount_management(polySelect).getOutput(0))

    count=0
    timePrev=time.time()
    with arcpy.da.SearchCursor(polySelect,["OID@"]+fields) as curInput:
        for rowInput in curInput:
            linesTemp=arcpy.SelectLayerByLocation_management(lineSelect,"INTERSECT",rowInput[-1],"","NEW_SELECTION")
            geometry=arcpy.CopyFeatures_management(linesTemp,arcpy.Geometry())
            geometry.append(rowInput[-1].boundary())
            arcpy.FeatureToPolygon_management (geometry, "in_memory\polygons_init")
            arcpy.Clip_analysis ("in_memory\polygons_init", rowInput[-1], "in_memory\polygons_clip")
            with arcpy.da.SearchCursor("in_memory\polygons_clip","SHAPE@") as curPoly:
                newGeom=[]
                for rowP in curPoly:
                    if not rowP[0].disjoint(rowInput[-1]):
                        newGeom.append(rowP[0])
            arcpy.Delete_management("in_memory")

            with arcpy.da.InsertCursor(outputPoly, ["OLD_OID"]+fields) as insCur:
                for geom in newGeom:
                    insertFeature=[r for r in rowInput[:-1]]
                    insertFeature.append(geom)
                    insCur.insertRow(insertFeature)
            count+=1
            if int(time.time()-timePrev)%5==0 and int(time.time()-timePrev)>0:
                timePrev=time.time()
                arcpy.AddMessage("\r{0}% done, {1} features processed".format(int(count*100/totalFeatures),int(count)))

def main():
    arcpy.env.overwriteOutput = True
    arcpy.env.XYTolerance = "0.1 Meters"

    inputPoly = arcpy.GetParameterAsText(0) # required
    inputLines = arcpy.GetParameterAsText(1) # required
    linesQuery = arcpy.GetParameterAsText(2) # optional
    outPoly = arcpy.GetParameterAsText(3) # optional

    if inputPoly=="":
        inputPoly=outFeatureClass2
    if arcpy.Exists(inputPoly):
        arcpy.AddMessage("Input polygons: "+inputPoly)
    else:
        arcpy.AddError("Input polygons layer %s is invalid" % (inputPoly))

    if inputLines=="":
        inputLines=allLines
    if arcpy.Exists(inputLines):
        arcpy.AddMessage("Input lines: "+inputPoly)
    else:
        arcpy.AddError("Input lines layer %s is invalid" % (inputLines))

    if linesQuery=="":
        arcpy.AddMessage("Performing without query")

    if outPoly == "":
        arcpy.AddMessage("Output will be created at the same location as input polygons layer is.")

    splitPolygonsWithLines(inputPoly, inputLines, linesQuery, outPoly)

if __name__ == "__main__":
    main()
print ""
print "Done"