Я хочу разделить многоугольник слева направо, основываясь на срезанных областях многоугольника, используя arcpy в ArcGIS

Question

Мне дан многоугольник с известными координатами. Пользователь может ввести области, которые он хотел бы вырезать из этого многоугольника. Например, если площадь составляет 440084,489 кв. М, пользователь может выбрать обрезку полигона при 440000, 84 и 0,489 кв. М или другие требуемые значения. (В более сложной версии пользователь также может выбрать направление (ie. Запад-Восток, Север-Юг и др.)). Затем скрипт проанализирует многоугольник и найдет координаты, в которых будут находиться указанные c области многоугольника. Я анализирую многоугольник, рассматривая направление Запад-Восток, и по оси X я подумал найти значение X, где вертикальная линия, которая будет разрезать многоугольник, будет l ie. Я пытаюсь найти расположение значения X, чтобы, если я обрежу там многоугольник, левая сторона имела бы первое значение в списке. Я должен достичь 5 мм кв. М ошибки, что я сейчас думаю, слишком малая величина. Я нахожу место, где его можно разрезать, непрерывно разрезая половину линии Запад-Восток. Если область слева выше требуемой, я снова разрезал линию пополам. Если область слева становится меньше необходимой площади, я делю линию пополам и принимаю во внимание правую сторону линии разреза. Я считаю, что это работает быстрее, чем увеличение значения X. Вот мой код:

    import arcpy
import os
import time

time_start = int(round(time.time() * 1000))
poly = r'C:\Users\Paul-PC\Documents\ArcGIS\New File Geodatabase.gdb\parcela'
direction = "E-V"
poly_nou = r'C:\Users\Paul-PC\Documents\ArcGIS\New File Geodatabase.gdb\tabel_pardon'
out_path = os.path.split(poly_nou)[0]
name =  os.path.split(poly_nou)[1]
spatial_ref = arcpy.Describe(poly).spatialReference
if arcpy.Exists(poly_nou):
    print("Verified {} exists".format(poly_nou))
else:
    arcpy.management.CreateFeatureclass(out_path, name, geometry_type="POLYGON", template=poly, spatial_reference=spatial_ref)
input_areas = "300000 100000 40000 800 40.089"

areas = []
b = input_areas.split(" ")
areas = []
for item in b:
    areas.append(float(item))

with arcpy.da.SearchCursor(poly, ["SHAPE@"]) as pcursor:
    for prow in pcursor:
        polygon = prow[0]  # polygon to cut
        polyg_area = float("%.5f" %polygon.area)
        #print(polyg_area)
        e = polygon.extent  # bounding extent of polygon
        #print(float("%.5f" %e.XMin), float("%.5f" %e.YMin), float("%.5f" %e.XMax), float("%.5f" %e.YMax))
del pcursor
print("Min x of polygon: ", e.XMin, " Max x of Polygon: ", e.XMax)
print(polyg_area)


cutpoly = polygon
icursor = arcpy.da.InsertCursor(poly_nou, ["SHAPE@"]) #arcpy.GetParameterAsText(4)

for i in areas[:len(areas)-1]:
    print(i)

    # define the extreme coordinates of the remaining polygon
    cutpoly_e = cutpoly.extent
    minX = cutpoly_e.XMin
    maxX = cutpoly_e.XMax
    ymax = cutpoly_e.YMax
    ymin = cutpoly_e.YMin

    pntarray = arcpy.Array()
    #print(float("%.5f" %(maxX)) - float("%.5f" %(minX)))
    midX = (float("%.5f" %(minX)) + ((float(float("%.5f" %(maxX)) - float("%.5f" %(minX)))) / 2))
    print(type(maxX), type(minX))
    jumatate = ((float("%.5f" %maxX)) - float("%.5f" %(minX))) / 2
    print("Mid point: ", midX, "jumatate: ", jumatate)
    pntarray.add(arcpy.Point(midX, ymin))
    pntarray.add(arcpy.Point(midX, ymax))
    pline = arcpy.Polyline(pntarray, arcpy.SpatialReference(31700))
    # cut polygon and get split-parts
    # part 0 is on the right side and part 1 is on the left side of the cut
    cutlist = cutpoly.cut(pline)
    temp_polyg = float("%.5f" %cutlist[0].area)
    print("Temp POlygon area: ", float("%.5f" %temp_polyg))
    left_polyg = cutlist[0]
    print("Left polyg XMin: ", float("%.5f" %left_polyg.extent.XMin), "Left Polyg XMAx: ", float("%.5f" %left_polyg.extent.XMax))
    print("Left polyg area: ", float("%.5f" %left_polyg.area))
    right_polyg = cutlist[1]
    print("Right polyg XMin: ", right_polyg.extent.XMin, "RIght Polyg XMAx: ", right_polyg.extent.XMax)

    new_cut_X = 0

    while abs(temp_polyg-i)>0.005 :

        # construct NS line
        while (temp_polyg < i-0.005):
            print(" temp polyg ", temp_polyg, " i area: ", i)
            print("Searched area is larger than temp polygon")
            right_polyg_e = right_polyg.extent
            r_poly_xMin = float("%.5f" %right_polyg_e.XMin)
            #r_poly_xMax = float("%.5f" %right_polyg_e.XMax)
            #r_poly_yMax = float("%.5f" %right_polyg_e.YMax)
            #r_poly_yMin = float("%.5f" %right_polyg_e.YMin)
            if (jumatate <= 0.001):
                jumatate = 0.0001
            else:
                jumatate = jumatate / 2
            print("Right poly x_min: ", r_poly_xMin)
            new_cut_X = r_poly_xMin + jumatate
            print(" new cut X: ", float("%.5f" %new_cut_X) , " jumatate: ", jumatate)

            pntarray = arcpy.Array()
            #newX = float("%.5f" % (minX + (r_poly_xMin - r_poly_xMin) / 2))
            #print(midX)
            pntarray.add(arcpy.Point(new_cut_X, ymin))
            pntarray.add(arcpy.Point(new_cut_X, ymax))

            pline = arcpy.Polyline(pntarray, arcpy.SpatialReference(31700))
            # cut polygon and get split-parts
            # part 0 is on the right side and part 1 is on the left side of the cut
            cutlist = cutpoly.cut(pline)
            temp_polyg = float("%.5f" %cutlist[0].area) #aici trebuie sa preia doar jumatatea stanga minus aria deja cercetata
            print("new Temp POlygon area: ", float("%.5f" %temp_polyg))
            left_polyg = cutlist[0]
            print("new left polygon area: ", float("%.5f" %left_polyg.area))
            right_polyg = cutlist[1]

        while (temp_polyg > i+0.005):
            print("Searched area is smaller than temp polygon")
            l_poly_e = left_polyg.extent
            #print(" new cut X: ", new_cut_X , " jumatate: ", jumatate)
            print("Left polygon area: ", left_polyg.area)
            # l_poly_xMin = float("%.5f" %l_poly_e.XMin)
            l_poly_XMax = float("%.5f" %l_poly_e.XMax)
            # l_poly_yMax = float("%.5f" %l_poly_e.YMax)
            # l_poly_yMin = float("%.5f" %l_poly_e.YMin)
            if(jumatate<=0.001):
                jumatate = 0.00001
            else:
                jumatate = jumatate/2
            new_cut_X = l_poly_XMax - jumatate
            print(" new cut X: ", float("%.5f" %new_cut_X), " jumatate: ", jumatate)

            pntarray = arcpy.Array()
            #midX = float("%.5f" % (cutpoly.extent. + (maxX - minX) / 2))
            pntarray.add(arcpy.Point(new_cut_X, ymin))
            pntarray.add(arcpy.Point(new_cut_X, ymax))

            pline = arcpy.Polyline(pntarray, arcpy.SpatialReference(31700))
            # cut polygon and get split-parts
            # part 0 is on the right side and part 1 is on the left side of the cut
            cutlist = cutpoly.cut(pline) #aici trebuie sa preia doar jumatatea dreapta minus aria deja cercetata
            temp_polyg = float("%.5f" %cutlist[0].area)
            print("Temp POlygon area: ", float("%.5f" %temp_polyg))
            left_polyg = cutlist[0]
            print("left polygon area: ", float("%.5f" %left_polyg.area))
            right_polyg = cutlist[1]

        if (jumatate<=0.0001):
            break
    print("found point cut for area: ", i, "left polygon xMax: ", float("%.5f" %left_polyg.extent.XMax), " left polygon area: ", left_polyg.area, " right polygon area: ", right_polyg.area)
    cutpoly = right_polyg
    print(right_polyg.area)
    #with arcpy.da.InsertCursor(poly_nou, ['SHAPE@XY']) as cursor:
    #icursor.insertRow([left_polyg])

# insert last cut remainder
# icursor.insertRow([right_polyg])
del icursor

time_end = int(round(time.time() * 1000))
print(time_end - time_start)

Он продолжает сокращать "jumatate" и добавлять или вычитать из оставшихся полигонов горизонтальную линию. Однако в какой-то момент он выходит за пределы значения 0,006, разрезая «джумат» пополам и переходит в бесконечный l oop, где он показывает те же значения.