Upload New File

gml_extrusion/build_up_geonmetries.py

+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)
\ No newline at end of file