summer_course_2024/imports/geometry/citygml.py
2021-06-01 18:31:50 -04:00

284 lines
13 KiB
Python

"""
CityGml module parses citygml files and import the geometry into the city model structure
SPDX - License - Identifier: LGPL - 3.0 - or -later
Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
"""
import numpy as np
import xmltodict
from city_model_structure.city import City
from city_model_structure.building import Building
from city_model_structure.attributes.surface import Surface
from helpers.geometry_helper import GeometryHelper
from city_model_structure.attributes.polygon import Polygon
class CityGml:
"""
CityGml class
"""
def __init__(self, path):
self._city = None
with open(path) as gml:
# Clean the namespaces is an important task to prevent wrong ns:field due poor citygml implementations
self._gml = xmltodict.parse(gml.read(), process_namespaces=True, xml_attribs=True, namespaces={
'http://www.opengis.net/gml': None,
'http://www.opengis.net/citygml/1.0': None,
'http://www.opengis.net/citygml/building/1.0': None,
'http://schemas.opengis.net/citygml/building/1.0/building.xsd': None,
'http://www.opengis.net/citygml/appearance/1.0': None,
'http://schemas.opengis.net/citygml/appearance/1.0/appearance.xsd': None,
'http://www.opengis.net/citygml/relief/1.0': None,
'http://schemas.opengis.net/citygml/relief/1.0/relief.xsd': None,
'http://www.opengis.net/citygml/generics/1.0': None,
'http://www.w3.org/2001/XMLSchema-instance': None,
'urn:oasis:names:tc:ciq:xsdschema:xAL:2.0': None,
'http://www.w3.org/1999/xlink': None,
'http://www.opengis.net/citygml/relief/2.0': None,
'http://www.opengis.net/citygml/building/2.0': None,
'http://www.opengis.net/citygml/building/2.0 http://schemas.opengis.net/citygml/building/2.0/building.xsd '
'http://www.opengis.net/citygml/relief/2.0 http://schemas.opengis.net/citygml/relief/2.0/relief.xsd" '
'xmlns="http://www.opengis.net/citygml/2.0': None,
'http://www.opengis.net/citygml/2.0': None
}, force_list=('cityObjectMember', 'curveMember', 'boundedBy', 'surfaceMember'))
self._city_objects = None
self._geometry = GeometryHelper()
for bound in self._gml['CityModel']['boundedBy']:
envelope = bound['Envelope']
if '#text' in envelope['lowerCorner']:
self._lower_corner = np.fromstring(envelope['lowerCorner']['#text'], dtype=float, sep=' ')
self._upper_corner = np.fromstring(envelope['upperCorner']['#text'], dtype=float, sep=' ')
else:
self._lower_corner = np.fromstring(envelope['lowerCorner'], dtype=float, sep=' ')
self._upper_corner = np.fromstring(envelope['upperCorner'], dtype=float, sep=' ')
if '@srsName' in envelope:
self._srs_name = envelope['@srsName']
@property
def content(self):
"""
CityGml raw content
:return: str
"""
return self._gml
@property
def city(self) -> City:
"""
City model structure enriched with the geometry information
:return: City
"""
if self._city is None:
# todo: refactor this method to clearly choose the gml type
self._city = City(self._lower_corner, self._upper_corner, self._srs_name)
i = 0
building_part = None
for o in self._gml['CityModel']['cityObjectMember']:
i += 1
lod = 0
surfaces = []
if 'lod1Solid' in o['Building']:
lod += 1
surfaces = CityGml._lod1_solid(o)
elif 'lod1MultiSurface' in o['Building']:
lod += 1
surfaces = CityGml._lod1_multi_surface(o)
elif 'lod2MultiSurface' in o['Building']:
# todo: check if this is a real case or a miss-formed citygml
lod = 2
surfaces = surfaces + CityGml._lod2_solid_multi_surface(o)
else:
for bound in o['Building']['boundedBy']:
surface_type = next(iter(bound))
if 'lod2MultiSurface' in bound[surface_type]:
lod = 2
surfaces = surfaces + CityGml._lod2(bound)
if 'lod3Solid' in o['Building']:
lod += 4
if 'lod4Solid' in o['Building']:
lod += 8
lod_terrain_str = 'lod' + str(lod) + 'TerrainIntersection'
terrains = []
if lod_terrain_str in o['Building']:
terrains = self._terrains(o, lod_terrain_str)
function = None
year_of_construction = None
if 'consistsOfBuildingPart' in o['Building']:
if 'BuildingPart' in o['Building']['consistsOfBuildingPart']:
name = o['Building']['consistsOfBuildingPart']['BuildingPart']['name']
if 'yearOfConstruction' in o['Building']['consistsOfBuildingPart']['BuildingPart']:
year_of_construction = o['Building']['consistsOfBuildingPart']['BuildingPart']['yearOfConstruction']
if 'function' in o['Building']['consistsOfBuildingPart']['BuildingPart']:
function = o['Building']['consistsOfBuildingPart']['BuildingPart']['function']
else:
name = o['Building']['@id']
if 'yearOfConstruction' in o['Building']:
year_of_construction = o['Building']['yearOfConstruction']
if 'function' in o['Building']:
function = o['Building']['function']
self._city.add_city_object(Building(name, lod, surfaces, year_of_construction, function, self._lower_corner,
terrains))
return self._city
def _terrains(self, city_object, lod_terrain_str):
try:
curves = [c['LineString']['posList']['#text']
for c in city_object['Building'][lod_terrain_str]['MultiCurve']['curveMember']]
except TypeError:
curves = [c['LineString']['posList']
for c in city_object['Building'][lod_terrain_str]['MultiCurve']['curveMember']]
terrains = []
for curve in curves:
curve_points = np.fromstring(curve, dtype=float, sep=' ')
curve_points = self._geometry.to_points_matrix(curve_points)
terrains.append(curve_points)
return terrains
@staticmethod
def _lod1_solid(o):
try:
solid_points = [CityGml._solid_points(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']['#text']))
for s in o['Building']['lod1Solid']['Solid']['exterior']['CompositeSurface']['surfaceMember']]
except TypeError:
solid_points = [CityGml._solid_points(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']))
for s in o['Building']['lod1Solid']['Solid']['exterior']['CompositeSurface']['surfaceMember']]
return [Surface(Polygon(sp),Polygon(sp)) for sp in solid_points]
@staticmethod
def _lod1_multi_surface(o):
solid_points = [CityGml._solid_points(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']))
for s in o['Building']['lod1MultiSurface']['MultiSurface']['surfaceMember']]
return [Surface(Polygon(sp),Polygon(sp)) for sp in solid_points]
@staticmethod
def _lod2_solid_multi_surface(o):
if 'boundedBy' in o['Building']['consistsOfBuildingPart']['BuildingPart']:
if 'RoofSurface' in o['Building']['consistsOfBuildingPart']['BuildingPart']['boundedBy']:
if o['Building']['consistsOfBuildingPart']['BuildingPart']['boundedBy']['RoofSurface']['lod2MultiSurface'] != 'None':
polygons = [Polygon(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']))
for s in o['Building']['consistsOfBuildingPart']['BuildingPart']['boundedBy']['RoofSurface']['lod2MultiSurface']['MultiSurface']['surfaceMember']]
elif 'WallSurface' in o['Building']['consistsOfBuildingPart']['BuildingPart']['boundedBy']:
if o['Building']['consistsOfBuildingPart']['BuildingPart']['boundedBy']['WallSurface']['lod2MultiSurface'] != 'None':
polygons = [Polygon(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']))
for s in o['Building']['consistsOfBuildingPart']['BuildingPart']['boundedBy']['WallSurface']['lod2MultiSurface']['MultiSurface']['surfaceMember']]
else:
polygons = [Polygon(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']))
for s in o['Building']['lod2MultiSurface']['MultiSurface']['surfaceMember']]
return [Surface(p,p) for p in polygons]
@staticmethod
def _lod2_composite_surface(s):
solid_points = [CityGml._solid_points((CityGml._remove_last_point(sm['Polygon']['exterior']['LinearRing']['posList'])))
for sm in s['CompositeSurface']['surfaceMember']]
return [Surface(Polygon(sp),Polygon(sp)) for sp in solid_points]
@staticmethod
def _lod2_multi_surface(s, surface_type):
# todo: this need to be changed into surface bounded?
try:
solid_points = [CityGml._solid_points(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList']
['#text']))]
except TypeError:
solid_points = [CityGml._solid_points(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']
['posList']))]
return [Surface(Polygon(sp),Polygon(sp), surface_type=surface_type) for sp in solid_points]
@staticmethod
def _lod2(bound):
surfaces = []
for surface_type in iter(bound):
for s in bound[surface_type]['lod2MultiSurface']['MultiSurface']['surfaceMember']:
if 'CompositeSurface' in s:
surfaces = surfaces + CityGml._lod2_composite_surface(s)
else:
surfaces = surfaces + CityGml._lod2_multi_surface(s, surface_type)
return surfaces
@staticmethod
def _remove_last_point(points):
array = points.split(' ')
res = " "
return res.join(array[0:len(array) - 3])
@staticmethod
def _solid_points(coordinates) -> np.ndarray:
"""
Solid surface point matrix [[x, y, z],[x, y, z],...]
:parameter coordinates: string from file
:return: np.ndarray
"""
solid_points = np.fromstring(coordinates, dtype=float, sep=' ')
solid_points = GeometryHelper.to_points_matrix(solid_points)
return solid_points
@staticmethod
def _holes_points(holes_coordinates) -> [np.ndarray]:
"""
Holes surfaces point matrices [[[x, y, z],[x, y, z],...]]
:parameter holes_coordinates: strings from file
:return: [np.ndarray]
"""
holes_points = []
for hole_coordinates in holes_coordinates:
hole_points = np.fromstring(hole_coordinates, dtype=float, sep=' ')
hole_points = GeometryHelper.to_points_matrix(hole_points)
holes_points.append(hole_points)
return holes_points
@staticmethod
def _perimeter_points(coordinates, solid_points, holes_points) -> np.ndarray:
"""
Matrix of points of the perimeter in the same order as in coordinates [[x, y, z],[x, y, z],...]
:parameter coordinates: string from file
:parameter solid_points: np.ndarray points that define the solid part of a surface
:parameter holes_points: [np.ndarray] points that define each the holes in a surface
:return: np.ndarray
"""
if holes_points is None:
perimeter_points = solid_points
else:
_perimeter_coordinates = coordinates
for hole_points in holes_points:
_hole = np.append(hole_points, hole_points[0])
_closed_hole = ' '.join(str(e) for e in [*_hole[:]])
# add a mark 'M' to ensure that the recombination of points does not provoke errors in finding holes
_perimeter_coordinates = _perimeter_coordinates.replace(_closed_hole, 'M')
_perimeter_coordinates = _perimeter_coordinates.replace('M', '')
perimeter_points = np.fromstring(_perimeter_coordinates, dtype=float, sep=' ')
perimeter_points = GeometryHelper.to_points_matrix(perimeter_points)
# remove duplicated points
pv = np.array([perimeter_points[0]])
for point in perimeter_points:
duplicated_point = False
for p in pv:
if GeometryHelper().almost_equal(0.0, p, point):
duplicated_point = True
if not duplicated_point:
pv = np.append(pv, [point], axis=0)
perimeter_points = pv
return perimeter_points
@staticmethod
def _holes_polygons(holes_points) -> [Polygon]:
"""
hole surfaces, a list of hole polygons found in a surface
:parameter holes_points: [np.ndarray] that define each of the holes
:return: None, [] or [Polygon]
None -> not known whether holes exist in reality or not due to low level of detail of input data
[] -> no holes in the surface
[Polygon] -> one or more holes in the surface
"""
if holes_points is None:
holes_polygons = None
else:
holes_polygons = []
for hole_points in holes_points:
holes_polygons.append(Polygon(hole_points))
return holes_polygons