diff --git a/gml_extrusion/build_up_geonmetries.py b/gml_extrusion/build_up_geonmetries.py new file mode 100644 index 0000000000000000000000000000000000000000..ee956678ce150de7bf2fb5154bebbbb328ebccba --- /dev/null +++ b/gml_extrusion/build_up_geonmetries.py @@ -0,0 +1,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) \ No newline at end of file