Missing indices when extracting vertices and edges from polygon

Question

Given a polygon with holes

from shapely.geometry import Polygon
import geopandas as gpd

# Polygon exterior:
p = [[20,767],[54,744],[107,707],
 [190,654],[265,609],[363,548],
 [462,484],[514,447],[603,389],
 [682,337],[726,310],[757,290],
 [786,277],[820,259],[843,249],
 [881,231],[921,215],[975,197],
 [1048,174],[1089,163],[1141,152],
 [1212,137],[1270,121],[1271,64],
 [1207,78],[1163,89],[1096,103],
 [1048,115],[1001,129],[949,144],
 [905,157],[874,170],[830,187],
 [781,208],[730,236],[696,255],
 [652,282],[606,306],[561,340],
 [512,370],[478,393],[436,418],
 [385,453],[330,490],[285,521],
 [229,566],[183,603],[123,652],
 [70,698],[13,749]]

# Define interior "holes":
interiors = {}
interiors[0] = [[290,543],[301,560],[393,501],[377,482]]
interiors[1] = [[507,392],[549,363],[553,367],[572,352],[588,372],[522,415]]
interiors[2] = [[599,340],[636,316],[648,334],[612,357]]
interiors[3] = [[714,262],[727,284],[821,238],[811,215]]
interiors[4] = [[850,218],[935,185],[937,187],[850,221]]
interiors[5] = [[959,159],[1066,129],[1071,146],[966,177]]
interiors[6] = [[1119,133],[1175,122],[1178,123],[1119,134]]
interiors[7] = [[1211,102],[1266,91],[1267,97],[1212,108],[1211,102]]

i_p = {k: Polygon(v) for k, v in interiors.items()}


zone = Polygon(p, [zone.exterior.coords for zone in i_p.values() \
                    if zone.within(Polygon(p)) is True])

polygon = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[zone])

Extract the vertices and edges as a i) coordinate list and ii) list of indices that define the edge (end points of the line segment).

l = 0    # this will be a number ~ len(gdf)
dps = 2

all_coords = []
idx = []

for i, row in polygon.iterrows():
    exterior = row.geometry.exterior
    ring = list(exterior.coords)
    ring.pop() # pop out duplicate vertice
    if exterior.is_ccw == False:
        #-- to get proper orientation of the normals
        ring.reverse()

    for j, v in enumerate(ring):
        idx.append([l+j, l+j+1])
        rounded_x = v[0]# round(float(v[0]), dps)
        rounded_y = v[1] #round(float(v[1]), dps)
        all_coords.append([rounded_x, rounded_y])
        
    idx.append([l+j+1, l])
    l = len(all_coords) + l

    for each, inte in enumerate(row.geometry.interiors):
        ring = list(inte.coords)
        ring.pop()
        if inte.is_ccw == True:
            #-- to get proper orientation of the normals
            ring.reverse()
        
        for p in ring:
            idx.append([l+each, l+each+1])
            rounded_x = p[0] #round(float(v[0]), 2)
            rounded_y = p[1] #round(float(v[1]), 2)
            all_coords.append([rounded_x, rounded_y])#, rounded_z]) 
    
        idx.append([l+each+1, l])
        l = len(all_coords) + l
      
print(all_coords)
print('')
print(idx)

[[13.0, 749.0], [70.0, 698.0], [123.0, 652.0], [183.0, 603.0], [229.0, 566.0], [285.0, 521.0], [330.0, 490.0], [385.0, 453.0], [436.0, 418.0], [478.0, 393.0], [512.0, 370.0], [561.0, 340.0], [606.0, 306.0], [652.0, 282.0], [696.0, 255.0], [730.0, 236.0], [781.0, 208.0], [830.0, 187.0], [874.0, 170.0], [905.0, 157.0], ... [966.0, 177.0], [1071.0, 146.0], [1066.0, 129.0], [959.0, 159.0], [1119.0, 134.0], [1178.0, 123.0], [1175.0, 122.0], [1119.0, 133.0], [1212.0, 108.0], [1267.0, 97.0], [1266.0, 91.0], [1211.0, 102.0]]

[[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14], [14, 15], [15, 16], [16, 17], [17, 18], [18, 19], [19, 20], [20, 21], [21, 22], [22, 23], [23, 24], [24, 25], [25, 26], [26, 27], [27, 28], [28, 29], [29, 30], [30, 31], [31, 32], [32, 33], [33, 34], [34, 35], [35, 36], [36, 37], [37, 38], [38, 39], [39, 40], [40, 41], [41, 42], [42, 43], [43, 44], [44, 45], [45, 46], [46, 47], [47, 48], [48, 49], [49, 50], [50, 0], [50, 51], [51, 52], [52, 53], [53, 54], [54, 50], [104, 105], [105, 106], [106, 107], [107, 108], [108, 109], [109, 110], [110, 104], [164, 165], [165, 166], [166, 167], [167, 168], [168, 164], [228, 229], [229, 230], [230, 231], [231, 232], [232, 228], [296, 297], [297, 298], [298, 299], [299, 300], [300, 296], [368, 369], [369, 370], [370, 371], [371, 372], [372, 368], [444, 445], [445, 446], [446, 447], [447, 448], [448, 444], [524, 525], [525, 526], [526, 527], [527, 528], [528, 524]]

Notice print(idx) goes ....[50, 51], [51, 52], [52, 53], [53, 54], [54, 50], [104, 105], [105, 106], [106, 107], [107, 108], [108, 109], [109, 110], [110, 104]... and then again ...[167, 168], [168, 164], [228, 229], [229, 230], [230, 231], [231, 232], [232, 228], [296, 297]...

Where are the missing vertices? Why is the loop skipping?

Answer 1

You have a problem with your for loops.

As there is only one polygon:

exterior = polygon.geometry[0].exterior
interior = polygon.geometry[0].interiors
ext_coords = list(exterior.coords)[:-1] # eliminate last value

all_coords = []
idx = []
# exterior
for i,coord  in enumerate(ext_coords):
    idx.append([i,i+1])
    all_coords.append(coord)
# interior
l = len(idx)
for int in interior:
   int_coords=list(int.coords)[:-1]
   for i in int_coords:
     idx.append([l,l+1])
     l = l+1
     all_coords.append(i)
print(idx)
[[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14], [14, 15], [15, 16], [16, 17], [17, 18], [18, 19], [19, 20], [20, 21], [21, 22], [22, 23], [23, 24], [24, 25], [25, 26], [26, 27], [27, 28], [28, 29], [29, 30], [30, 31], [31, 32], [32, 33], [33, 34], [34, 35], [35, 36], [36, 37], [37, 38], [38, 39], [39, 40], [40, 41], [41, 42], [42, 43], [43, 44], [44, 45], [45, 46], [46, 47], [47, 48], [48, 49], [49, 50], [50, 51], [51, 52], [52, 53], [53, 54], [54, 55], [55, 56], [56, 57], [57, 58], [58, 59], [59, 60], [60, 61], [61, 62], [62, 63], [63, 64], [64, 65], [65, 66], [66, 67], [67, 68], [68, 69], [69, 70], [70, 71], [71, 72], [72, 73], [73, 74], [74, 75], [75, 76], [76, 77], [77, 78], [78, 79], [79, 80], [80, 81], [81, 82], [82, 83], [83, 84]]