Dzielenie wielokąta linią za pomocą ArcPy?

Patrząc na buforowanie za pomocą fizycznej bariery za pomocą ArcGIS na komputery? , przyszło mi do głowy, że nie jestem pewien, jak można by użyć narzędzi geoprzetwarzania w ArcGIS do programowego podzielenia wielokąta linią.

Ręcznie użyłbyś narzędzia Wytnij wielokąty lub Podziel wieloboki na pasku narzędzi Topologia , ale w jaki sposób możesz wykonać to samo zadanie za pomocą narzędzia do tworzenia modeli lub narzędzi do tworzenia skryptów w Pythonie?

Od razu myślę o wszystkich narzędziach do analizy, takich jak Unia, Tożsamość itp., Ale wszystkie są narzędziami Wielokąt-Wielokąt, a NIE narzędziami Linia Wielokąt. Nawet narzędzie Podziel to Wielokąt-Wielokąt.

Jakieś pomysły?

Korzystając z ET Geowizard, możesz uzyskać dostęp do kodu dla narzędzia Podziel wieloboki za pomocą polilinii:

Oto link do skryptu .

Możesz również użyć ArcObjects, aby to zrobić:

Wytnij fragment wielokąta

Możesz także użyć opisanej tutaj metody buforowania z jednej strony .

Po tym fakcie stworzyłem własne narzędzie ModelBuilder. Zapomniałem o tym pytaniu i zamieściłem moje rozwiązanie na inne podobne pytanie . Dla kompletności jest to odpowiedź na odpowiedź:

Pomyślałem, że musi być na to sposób, więc stworzyłem coś, co uważam za całkiem dobre rozwiązanie. Zamieściłem go na stronie zasobów ArcGIS w społeczności-> Techniczne-> Analiza i geoprzetwarzanie-> Analiza-> Galeria.

Narzędzie nazywa się Podziel wieloboki z liniami i wymaga licencji ArcInfo z powodu niektórych narzędzi używanych w modelu. Zasadniczo to, co zrobiłem, to utworzenie minimalnego obwiedni dla wielokątów i przedłużenie linii do nich. Korzystając z voodoo ModelBuilder, mogłem przekształcić linie w wielokąty, które następnie wykorzystałem Tożsamość do podzielenia oryginalnych polis.

Sprawdź to i sprawdź, czy to działa dla Ciebie. W moich (ograniczonych) testach zachował atrybuty oryginalnych wielokątów i podzielił tylko istniejące wielokąty.

Jeśli chcesz wyjść poza ArcGIS, skorzystaj z geom.splitpolysbylines .

Osobiście nigdy nie korzystałem z niego w programie, ale myślę, że możesz uzyskać dostęp do tej linii komunikacyjnej za pomocą Pythona, więcej informacji znajdziesz w pomocy.

jeśli nie masz żadnych problemów z wysoką dokładnością, buforujesz linię z minimalną odległością np. (0,002, myślę, że powinna być lepsza niż dokładność twojej klasy obiektów), a następnie zastosuj narzędzie do wymazywania wielokąta przez buforowaną linię.

Zaktualizowano kod łuku w celu dzielenia wielokątów w kierunku poziomym lub pionowym przy użyciu wartości procentowych

'''-------------------------------------------------------------------------------------------
 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"