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
            if building_id not in bldg_IDs:
                continue
        tree.write(output_file_path, pretty_print=True, xml_declaration=True, encoding='utf-8')
        return print('Building extrusion completed and saved in:\n'+output_file_path)