build_up_geonmetries.py 14.9 KB
Newer Older
Ehlers's avatar
Ehlers committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
from lxml import etree
import copy
import uuid
import os
#os.chdir(r'C:\Users\ge29duf\Documents\02_Forschung\P62\Tool_ne\gml_changes')

class bldg_extrusion:
    def __init__(self, gml_par, file_path, output_file_path, floor_high, residential_code):
        self.gml_par = gml_par
        self.file_path = file_path
        self.output_file_path = output_file_path
        self.floor_high = floor_high
        self.residential_code = residential_code
        
    # Function to extrude buildings
    def bldg_extrusion(self):
        gml_par = self.gml_par
        file_path = self.file_path
        output_file_path = self.output_file_path
        floor_high = self.floor_high
        bldg_function = self.residential_code

        extrusion_height = 1000.0  # Default extrusion height if measuredHeight is not found
        

        # get building IDs from gml_par and number of floors to extrude from gml_par['extruded']
        bldg_IDs = gml_par['id'][gml_par['extruded'] > 0].tolist()
        #nr_of_floors = gml_par['extruded'][gml_par['extruded'] > 0].tolist()
        
        tree = etree.parse(file_path)

        # Define the namespaces
        namespaces = {
            'gml': 'http://www.opengis.net/gml',
            'bldg': 'http://www.opengis.net/citygml/building/1.0',
            'core': 'http://www.opengis.net/citygml/1.0',
            'xlink': 'http://www.w3.org/1999/xlink'
        }
        # Register the namespace
        for prefix, uri in namespaces.items():
            etree.register_namespace(prefix, uri)

        # Loop through each Building element
        buildings = tree.xpath('//bldg:Building', namespaces=namespaces)
        for building in buildings:
            # Check if the building ID is in the list of IDs to extrude
            building_id = building.get('{http://www.opengis.net/gml}id')

            if building_id in bldg_IDs or bldg_IDs == 'all':
                extra_building_hight = float(gml_par[gml_par['id']==building_id]['extruded']) * floor_high

                function = building.xpath('./bldg:function/text()', namespaces=namespaces)[0]
                if function.strip() == bldg_function:
                    # print('function is residential: '+function)
                    # Find the measuredHeight for the building, use extrusion_height if not found
                    measured_height_list = building.xpath('./bldg:measuredHeight/text()', namespaces=namespaces)
                    measured_height = float(measured_height_list[0]) if measured_height_list else extrusion_height

                    # Find all WallSurface elements
                    wall_surfaces = building.xpath('.//bldg:WallSurface', namespaces=namespaces)

                    # Find all GroundSurface elements
                    ground_surfaces = building.xpath('.//bldg:GroundSurface', namespaces=namespaces)

                    # Find all RoofSurface elements
                    roof_surfaces = building.xpath('.//bldg:RoofSurface', namespaces=namespaces)

                    # change roof type to flat <bldg:roofType> with '1000'
                    roof_type = building.xpath('./bldg:roofType/text()', namespaces=namespaces)
                    if roof_type[0] != '1000':
                        roof_type_element = building.xpath('./bldg:roofType', namespaces=namespaces)[0]
                        roof_type_element.text = '1000'
                        #print('roof type changed to flat')
                    
                    # change measured height to + extra_building_hight
                    measured_height_element = building.xpath('./bldg:measuredHeight', namespaces=namespaces)[0]
                    measured_height_element.text = str(measured_height + extra_building_hight)

                    # deleting all <bldg:lod2Solid> elements in building
                    lod2Solid = building.xpath('.//bldg:lod2Solid', namespaces=namespaces)
                    if lod2Solid:
                        building.remove(lod2Solid[0])
                        #print('removed lod2Solid')

                    for ground_surface in ground_surfaces:
                        pos_lists_ground = ground_surface.xpath('.//gml:posList', namespaces=namespaces)
                        for pos_list_ground in pos_lists_ground:
                            coordinates = pos_list_ground.text.split()
                            coord_triplets = [(float(coordinates[j]), float(coordinates[j+1]), float(coordinates[j+2])) for j in range(0, len(coordinates), 3)]
                            extruded_coord_triplets = [(x, y, z + measured_height + extra_building_hight) for x, y, z in coord_triplets]
                            # Assuming one WallSurface for each ground edge for simplicity
                            for i, new_wall in enumerate(coord_triplets[:-1]):
                                
                                if i < len(wall_surfaces):  # Ensure there is a corresponding wall
                                    wall_surface = wall_surfaces[i]
                                    multi_surfaces = wall_surface.xpath('.//gml:MultiSurface', namespaces=namespaces)
                                    if multi_surfaces:
                                        multi_surface = multi_surfaces[0]
                                        pos_lists_wall = multi_surface.xpath('.//gml:posList', namespaces=namespaces)
                                        if pos_lists_wall:
                                            pos_list_wall = pos_lists_wall[0]
                                            # Generate and set new coordinates for the wall
                                            x, y, z = coord_triplets[i]
                                            next_x, next_y, next_z = coord_triplets[(i + 1) % len(coord_triplets)]
                                            extruded_x, extruded_y, extruded_z = extruded_coord_triplets[i]
                                            next_extruded_x, next_extruded_y, next_extruded_z = extruded_coord_triplets[(i + 1) % len(coord_triplets)]
                                            pos_list_wall.text = f"{x} {y} {z} {extruded_x} {extruded_y} {extruded_z} {next_extruded_x} {next_extruded_y} {next_extruded_z} {next_x} {next_y} {next_z} {x} {y} {z}"
                                            #pos_list_wall.text = f"{x} {y} {z} {next_x} {next_y} {next_z} {next_extruded_x} {next_extruded_y} {next_extruded_z} {extruded_x} {extruded_y} {extruded_z} {x} {y} {z}"

                                if i >= len(wall_surfaces)-1:
                                    # Create new wall surfaces for the additional edges
                                    new_wall_surface = copy.deepcopy(wall_surface)

                                    # Create a new <bldg:boundedBy> element with the 'bldg' namespace prefix
                                    bounded_by = etree.Element('{http://www.opengis.net/citygml/building/1.0}boundedBy')
                                    # change the id of the new wall surface to avoid duplicates; <bldg:WallSurface gml:id=
                                    # replace the last character of the id with the current index
                                    new_wall_surface.attrib['{http://www.opengis.net/gml}id'] = new_wall_surface.attrib['{http://www.opengis.net/gml}id'][:-1] + str(i)
                                    # change the id of the new multiSurface <gml:MultiSurface gml:id=
                                    new_wall_surface_multi_surfaces = new_wall_surface.xpath('.//gml:MultiSurface', namespaces=namespaces)
                                    if new_wall_surface_multi_surfaces:
                                        new_wall_surface_multi_surface = new_wall_surface_multi_surfaces[0]
                                        new_wall_surface_multi_surface.attrib['{http://www.opengis.net/gml}id'] = new_wall_surface_multi_surface.attrib['{http://www.opengis.net/gml}id'][:-1] + str(i)
                                        #change the id of the new <gml:Polygon gml:id=
                                        new_wall_surface_polygons = new_wall_surface.xpath('.//gml:Polygon', namespaces=namespaces)
                                        if new_wall_surface_polygons:
                                            new_wall_surface_polygon = new_wall_surface_polygons[0]
                                            new_wall_surface_polygon.attrib['{http://www.opengis.net/gml}id'] = new_wall_surface_polygon.attrib['{http://www.opengis.net/gml}id'][:-1] + str(i)
                                            # change the id of the new <gml:LinearRing gml:id=
                                            new_wall_surface_linear_rings = new_wall_surface.xpath('.//gml:LinearRing', namespaces=namespaces)
                                            if new_wall_surface_linear_rings:
                                                new_wall_surface_linear_ring = new_wall_surface_linear_rings[0]
                                                new_wall_surface_linear_ring.attrib['{http://www.opengis.net/gml}id'] = new_wall_surface_linear_ring.attrib['{http://www.opengis.net/gml}id'][:-1] + str(i)
                                                #print(new_wall_surface_linear_ring.attrib['{http://www.opengis.net/gml}id'])
                                    
                                    # Add the new wall surface as a child to the boundedBy element
                                    bounded_by.append(new_wall_surface)

                                    # Update the posList with the new coordinates
                                    new_wall_surface_pos_lists = new_wall_surface.xpath('.//gml:posList', namespaces=namespaces)
                                    if new_wall_surface_pos_lists:
                                        new_wall_surface_pos_list = new_wall_surface_pos_lists[0]
                                        x, y, z = coord_triplets[i]
                                        next_x, next_y, next_z = coord_triplets[(i + 1) % len(coord_triplets)]
                                        extruded_x, extruded_y, extruded_z = extruded_coord_triplets[i]
                                        next_extruded_x, next_extruded_y, next_extruded_z = extruded_coord_triplets[(i + 1) % len(coord_triplets)]
                                        new_wall_surface_pos_list.text = f"{x} {y} {z} {extruded_x} {extruded_y} {extruded_z} {next_extruded_x} {next_extruded_y} {next_extruded_z} {next_x} {next_y} {next_z} {x} {y} {z}"
                                        #new_wall_surface_pos_list.text = f"{x} {y} {z} {next_x} {next_y} {next_z} {next_extruded_x} {next_extruded_y} {next_extruded_z} {extruded_x} {extruded_y} {extruded_z} {x} {y} {z}"
                                        # Append the new wall surface to the building
                                        #print(new_wall_surface_pos_list.text)
                                        building.append(bounded_by)
                                        #print("Added additional wall surface.")
                                        #else:
                                            # Handle the case where no posList is found
                                            #pass
                                    #else:
                                        # Handle the case where no MultiSurface is found
                                        #pass

                            # replace roof surfaces with extruded ground surface
                            for y, roof_surface in enumerate(roof_surfaces):
                            
                                if y == 0:
                                # If MultiSurface exists, assume the first one is the primary surface
                                    multi_surfaces_roof = roof_surface.xpath('.//gml:MultiSurface', namespaces=namespaces)
                                    if multi_surfaces_roof:
                                        multi_surface_roof = multi_surfaces_roof[0]
                                        pos_lists_roof = multi_surface_roof.xpath('.//gml:posList', namespaces=namespaces)
                                        if pos_lists_roof:
                                            # Directly update the first posList found
                                            pos_list_roof = pos_lists_roof[0]
                                            # Set or update the posList text with the new extruded coordinates
                                            roof_coords_str = " ".join(f"{x} {y} {z}" for x, y, z in reversed(extruded_coord_triplets))
                                            pos_list_roof.text = roof_coords_str
                                            #print("Updated roof surface with new geometry.")

                                if y > 0:
                                    # Assuming 'roof_surface' is the lxml Element representing the RoofSurface you want to remove
                                    boundedBy = roof_surface.getparent()
                                    grandparent = boundedBy.getparent()

                                    # Remove the boundedBy element, which also removes the RoofSurface child
                                    grandparent.remove(boundedBy)
                                    #print("Removed RoofSurface and its parent boundedBy.")

                    # all polygon gml:id from the building surfaces are collected
                    polygon_ids = []
                    for building_surface in building:
                        polygon_ids.extend(building_surface.xpath('.//gml:Polygon/@gml:id', namespaces=namespaces))

                    # Create the root element
                    lod2Solid = etree.Element('{http://www.opengis.net/citygml/building/1.0}lod2Solid')

                    # Create the gml:Solid element
                    solid = etree.SubElement(lod2Solid, '{http://www.opengis.net/gml}Solid', attrib={'{http://www.opengis.net/gml}id': 'UUID_' + str(uuid.uuid4())})

                    # Create the gml:exterior and gml:CompositeSurface elements
                    exterior = etree.SubElement(solid, '{http://www.opengis.net/gml}exterior')
                    composite_surface = etree.SubElement(exterior, '{http://www.opengis.net/gml}CompositeSurface', attrib={'{http://www.opengis.net/gml}id': 'UUID_' + str(uuid.uuid4())})

                    # Create a gml:surfaceMember element for each polygon id
                    for polygon_id in polygon_ids:
                        etree.SubElement(composite_surface, '{http://www.opengis.net/gml}surfaceMember', attrib={'{http://www.w3.org/1999/xlink}href': '#' + polygon_id})
                    
                    # Append the lod2Solid to the building
                    building.append(lod2Solid)
                

                else:
                    print('function is not residential: '+function)
                    pass
                # Write the modified CityGML back to a new file
                tree.write(output_file_path, pretty_print=True, xml_declaration=True, encoding='utf-8')
            if building_id not in bldg_IDs:
                continue
        return print('Building extrusion completed and safed in:\n'+output_file_path)