Render Shapely multipolygon as a 3D prism and export in gltf2

shapely_to_3d_mesh.py

from shapely import wkb
from shapely.ops import triangulate
from shapely.geometry import LineString
import trimesh

# Hol-y building https://www.openstreetmap.org/relation/3868363
wkb_holes = bytes.fromhex(
    """
01060000000100000001030000000500000036000000BBA067EAA3D836414EC6A5B2124A5A415E4
E319AAED83641F70DA6AD114A5A4151BE9F9CB9D83641BE8D6992104A5A416011309ED4D8364162
1E6FF10D4A5A4193B19785E8D8364167FFCCF00B4A5A41C773A7A303D93641835F4A45094A5A414
C30164A2AD93641DE895A76054A5A415AFD1FE12AD93641F8F9DBE0054A5A41E3D577392BD93641
DCE9A4D1054A5A418D56936732D936419C1644430A4A5A4104797DAB31D936411FADCC4D0A4A5A4
1525BB0C436D93641C7EB1E8F0D4A5A4181CDD93D40D93641B0BFAE8E134A5A41CAB096DA40D936
41FB252684134A5A41AE5D770D4CD9364122937CA91A4A5A41753D48404BD93641066BB9BF1A4A5
A4194EAC97952D93641BB4E49831F4A5A41778656115AD936418B937A06244A5A41EAB77A805AD9
364156B49AF9234A5A419C177E6461D93641461F6465284A5A41EECFA6EF60D9364124C0EC6F284
A5A41C2813F1D61D9364108E3DD91284A5A41F1B6F0FC5ED936417E2706C3284A5A41884BFE505C
D936413312882D294A5A417F653EC759D93641D2BBF882294A5A41AF9AEFA657D93641A161F5B22
94A5A41033A623F55D9364173C79AE0294A5A41B175A37C52D936419EAB34002A4A5A41470F6F34
50D936415574683F2A4A5A4189040AF04FD93641C2EFA21E2A4A5A419B39AFFD36D936419EFF259
32C4A5A417E7DF3982ED93641CD6F7B632D4A5A41C0774CB82ED93641687435762D4A5A411215A5
E029D93641F49091E62D4A5A41368424FD29D93641F8D8FF042E4A5A41DF5E2BFA11D93641B78E4
05C304A5A41B40B06C411D93641F2C58042304A5A41CD61C98404D936414DFE3D91314A5A4111C9
7B7EEBD83641E2C0E1F8334A5A41E47A14ACEBD836412ED2551E344A5A418AD76D33E9D8364118C
E065A344A5A4104E65FE8E6D83641405A7DB6344A5A41588A90E4E4D83641CA1A34F9344A5A4180
DB9BAAE2D8364172F50A2F354A5A4136F6C49DDFD83641D9529069354A5A411569D12ADDD836411
D61A785354A5A41C98D76E5DAD836412BBF2CC0354A5A4181AF77ACDAD836411910EAA2354A5A41
5D3B3A2CDAD83641B4F6C9AF354A5A41E5CF4117CCD83641F1225DA22C4A5A410AA5423AB0D8364
17A28BFC61A4A5A411A269A0CA4D8364112202D24134A5A41C86D7181A4D83641F2474D17134A5A
41BBA067EAA3D836414EC6A5B2124A5A411400000019C01806C0D8364154B513D6134A5A4169A46
58FC7D83641CEAF6583184A5A413B915CE6D0D83641F6A5548C1E4A5A41048E17ECD6D83641120D
7FEF1D4A5A417FA89FD8D7D83641314AC67A1E4A5A41C897F2ACDBD836418FE275181E4A5A41DF3
4AACEDAD8364130A72E8D1D4A5A411886BBFDE3D8364136A81FAA1C4A5A418224AA7BF1D83641E3
BA7D471B4A5A418CCDF16BE8D836419EFE7452154A5A41244271F9E0D836411ECCC094104A5A41C
4470C08CAD836416570C2D5124A5A412CA7687CC9D83641AF66267F124A5A41A0D88CEBC9D83641
61CD9D74124A5A4185433721C9D83641701BE5FA114A5A41E5A8DB7CC4D83641F4055280124A5A4
1A89CA43EC5D836418C9D30F4124A5A415EC821CDC5D83641AFE7CDE3124A5A414042062EC6D836
41A514E736134A5A4119C01806C0D8364154B513D6134A5A410D00000083A5FBE4D9D8364106EE4
D6C244A5A41AB337F44E3D836412B56484C2A4A5A416366D4A6E6D83641505B407F2C4A5A414BB8
C8E9E3D83641482BD1C72C4A5A4132258D8FE5D83641CB5297ED2D4A5A4154A8043BE7D836415AE
94EC92D4A5A419E814510E7D8364193A1E0AA2D4A5A41C842C541FBD8364165339FBA2B4A5A4174
E0CA4209D936418B57A5552A4A5A417D9F243204D9364177EB2F21274A5A41A5CE87C1FAD836411
164191E214A5A41A7785935EDD83641FE96016E224A5A4183A5FBE4D9D8364106EE4D6C244A5A41
1800000061C8986AF7D836418A2F42570E4A5A41C858D740FFD836416BAEC40C134A5A41A8085C6
708D936416D4EB81C194A5A419E7FA23D16D93641A42942BB174A5A412D7240751FD936410CD884
D3164A5A4151E67DF51FD936413E707225174A5A414AFD19781FD93641310BFB2F174A5A418A01E
F5E20D93641424BBFB7174A5A41AE8B3EDE24D93641E423323F174A5A4124AE282224D936413668
BFB2164A5A41B07C04B323D936416AA273BE164A5A411EC0064123D936413067B774164A5A41E7B
CC14629D936413B7590DC154A5A410029970620D93641B604E0E90F4A5A41CCCF238018D93641A0
7A1B170B4A5A4186C25C5512D93641211C68A90B4A5A41634E1FD511D93641A9034F560B4A5A41B
AF4694112D93641C66CC64B0B4A5A416B2D078B11D936413F7639D30A4A5A413ADB6BD80CD93641
B9E874500B4A5A41FDCE349A0DD936412827FCC10B4A5A41B3FAB1280ED9364126521CB50B4A5A4
169262FB70ED93641A46B35080C4A5A4161C8986AF7D836418A2F42570E4A5A41170000005C6A79
BB11D93641F731C5E11E4A5A4154A1A30416D936414730B2A4214A5A41DA815DB414D93641A54BA
3C6214A5A41E944EB8314D93641729086A3214A5A41E7478F7712D936413D0863E0214A5A41E260
456A14D93641C0E6020C234A5A416A4D955116D93641850C52D0224A5A417810232116D9364144A
F66B5224A5A41B919B66B17D936415D13C78E224A5A41B92A0A071BD9364152C92CE0244A5A41F9
44F1EC1FD93641D7181F11284A5A413E791EB22DD93641680868C9264A5A41B40963F140D936413
5A9AFE3244A5A41E057CAC340D936414ECA15C4244A5A4103D6830B42D93641702BD3A6244A5A41
F412F63B42D93641C5E94CD3244A5A4107CE21B443D9364126EBD8AD244A5A412F24EBDD41D9364
13C1C1388234A5A419386EB453FD93641A2D746C7234A5A41C4AC62FA3BD93641314EE9B3214A5A
41D04B2E2332D936412D9EC29B1B4A5A419ABA00671FD93641D267D78A1D4A5A415C6A79BB11D93
641F731C5E11E4A5A41
""".replace(
        "\n", ""
    )
)

# wkb_holes = bytes.fromhex('0106000000010000000103000000020000000500000061CD67C002DB364196BAB7A7D5495A41E226F95343DB3641D8685634CF495A4114B60CA053DB364162272C1AD9495A417DB5472313DB3641407DD9E8DF495A4161CD67C002DB364196BAB7A7D5495A4105000000CCB922731FDB3641B31F738FD6495A4150D53A4C27DB36410FA64F55DB495A41446C0FE939DB3641957DED6AD9495A41BF50F70F32DB36418F4311A5D4495A41CCB922731FDB3641B31F738FD6495A41')
s = wkb.loads(wkb_holes)

# apply polygon triangulation to have triangles
# note that if there are holes they are covered too
internal_triangles = [tri for tri in triangulate(s) if tri.within(s)]


# render the internal and external triangles to show the structure
import shapely.geometry

mul = shapely.geometry.MultiPolygon(internal_triangles)
with open("check_internal.svg", "w") as fw:
    fw.write(mul._repr_svg_())

mul_holes = shapely.geometry.MultiPolygon(
    [tri for tri in triangulate(s) if not tri.within(s)]
)

with open("check_external.svg", "w") as fw:
    fw.write(mul_holes._repr_svg_())

with open("check_all_triangles.svg", "w") as fw:
    fw.write(shapely.geometry.MultiPolygon(triangulate(s))._repr_svg_())


# validate it with trimesh
# it will take care of merging overlapping vertices
vertices = []
faces = []
HEIGHT = 100.0

dbg_int_segments = []
dbg_ext_segments = []


for tri_num, tri in enumerate(internal_triangles):
    # it's 4 coordinates, the last is identical to the first to form a circuit
    for x, y in tri.exterior.coords[:-1]:
        vertices.append([x, y, 0.0])
    for x, y in tri.exterior.coords[:-1]:
        vertices.append([x, y, HEIGHT])
    # now 6 points have been added
    idx = tri_num * 6
    # bottom and top faces of the prism, use as is
    faces.append([idx + 2, idx + 1, idx])
    faces.append([idx + 3, idx + 4, idx + 5])

    # now for the 3 lateral faces check that they do not fall inside the base
    # before adding them, or we'd have an invalid topology
    for i1, i2 in ((0, 1), (1, 2), (2, 0)):
        if LineString([tri.exterior.coords[i1], tri.exterior.coords[i2]]).within(s):
            dbg_int_segments.append(
                LineString([tri.exterior.coords[i1], tri.exterior.coords[i2]])
            )
            continue
        dbg_ext_segments.append(
            LineString([tri.exterior.coords[i1], tri.exterior.coords[i2]])
        )
        # add the two triangles making up the exterior face
        # there are 2 equivalent ways to split it, just use one
        faces.append([idx + i1, idx + i2, idx + i1 + 3])
        faces.append([idx + i2, idx + i2 + 3, idx + i1 + 3])


with open("check_external_segments.svg", "w") as fw:
    fw.write(shapely.geometry.MultiLineString(dbg_ext_segments)._repr_svg_())

with open("check_internal_segments.svg", "w") as fw:
    fw.write(shapely.geometry.MultiLineString(dbg_int_segments)._repr_svg_())


mesh = trimesh.Trimesh(
    vertices=vertices,
    faces=faces,
)
# requires networkx installed
# trimesh.repair.fix_inversion(mesh)
# trimesh.repair.fix_winding(mesh)
# trimesh.repair.fix_normals(mesh, multibody=True)

# mesh.fill_holes()
mesh.export("mesh_topo.ply")
mesh.export("mesh_topo.glb")
print(mesh.euler_number)
# requires pyglet installed
mesh.show()