diff --git a/src/ifcopenshell-python/ifcopenshell/util/geolocation.py b/src/ifcopenshell-python/ifcopenshell/util/geolocation.py index 9c2fd9ba9a..c97558cd24 100644 --- a/src/ifcopenshell-python/ifcopenshell/util/geolocation.py +++ b/src/ifcopenshell-python/ifcopenshell/util/geolocation.py @@ -18,6 +18,7 @@ import math import numpy as np +import numpy.typing as npt import ifcopenshell import ifcopenshell.util.unit import ifcopenshell.util.element @@ -94,7 +95,7 @@ def xyz2enh( you are applying your own temporary false origin (such as when federating models for digital twins of large cities). - For most scenarios you should use ``auto_xyz2enh`` instead. + For most scenarios you should use :func:`auto_xyz2enh` instead. :param x: The X local engineering coordinate. :param y: The Y local engineering coordinate. @@ -131,7 +132,7 @@ def auto_xyz2enh( The necessary georeferencing map conversion is automatically detected from the IFC map conversion parameters present in the IFC model. If no map - conversion is present, then the Z coordinate is returned unchanged. + conversion is present, then the coordinates are returned unchanged. For IFC2X3, the map conversion is detected from the IfcProject's ePSet_MapConversion. See the "User Guide for Geo-referencing in IFC": @@ -144,7 +145,6 @@ def auto_xyz2enh( :param should_return_in_map_units: If true, the result is given in map units. If false, the result will be converted back into project units. :return: The global map coordinate eastings, northings, and height. - :rtype: tuple[float] """ parameters = get_helmert_transformation_parameters(ifc_file) if not parameters: @@ -170,6 +170,7 @@ def auto_enh2xyz(ifc_file, easting, northing, height, is_specified_in_map_units: :param easting: The global easting map coordinate provided in map units. :param northing: The global northing map coordinate provided in map units. :param height: The global height map coordinate provided in map units. + :param is_specified_in_map_units: True if the input eastings, northing, and height are in map units. :return: The local engineering XYZ coordinates in project length units. """ parameters = get_helmert_transformation_parameters(ifc_file) @@ -270,7 +271,7 @@ def z2e(z: float, orthogonal_height: float = 0.0, scale: float = 1.0, factor_z: This function is for advanced users as it allows you to specify your own orthogonal height offset and transformation parameters. - For most scenarios you should use ``auto_z2e`` instead. + For most scenarios you should use :func:`auto_z2e` instead. :param z: The Z local engineering coordinate provided in project length units. :param orthogonal_height: The orthogonal height offset to apply. @@ -310,7 +311,7 @@ def enh2xyz( you are applying your own temporary false origin (such as when federating models for digital twins of large cities). - For most scenarios you should use ``auto_enh2xyz`` instead. + For most scenarios you should use :func:`auto_enh2xyz` instead. :param e: The global easting map coordinate. :param n: The global northing map coordinate. @@ -342,7 +343,18 @@ def enh2xyz( return (x, y, z) -def local2global(matrix, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None): +def local2global( + matrix: npt.NDArray[np.float64], + eastings: float = 0.0, + northings: float = 0.0, + orthogonal_height: float = 0.0, + x_axis_abscissa: float = 1.0, + x_axis_ordinate: float = 0.0, + scale: float = 1.0, + factor_x: float = 1.0, + factor_y: float = 1.0, + factor_z: float = 1.0, +) -> npt.NDArray[np.float64]: """Manually convert a 4x4 matrix from local to global coordinates This function is for advanced users as it allows you to specify your own @@ -352,63 +364,20 @@ def local2global(matrix, eastings, northings, orthogonal_height, x_axis_abscissa you are applying your own temporary false origin (such as when federating models for digital twins of large cities). - No unit conversion is performed. + For most scenarios you should use :func:`auto_local2global` instead. :param matrix: A 4x4 numpy matrix representing local coordinates. - :type matrix: np.array :param eastings: The eastings offset to apply. - :type eastings: float :param northings: The northings offset to apply. - :type northings: float :param orthogonal_height: The orthogonal height offset to apply. - :type orthogonal_height: float :param x_axis_abscissa: The X axis abscissa (i.e. first coordinate) of the 2D vector that points to the local X axis when in map coordinates. - :type x_axis_abscissa: float :param x_axis_ordinate: The X axis ordinate (i.e. second coordinate) of the 2D vector that points to the local X axis when in map coordinates. - :type x_axis_ordinate: float :param scale: The combined scale factor to convert from local coordinates to map coordinates. - :type scale: float :return: A numpy 4x4 array matrix representing global coordinates. - :rtype: np.array """ - if scale is None: - scale = 1.0 - x = np.array([x_axis_abscissa, x_axis_ordinate, 0]) - x /= np.linalg.norm(x) - y = np.cross(np.array([0, 0, 1]), x) - intermediate = ( - np.array( - [ - [x[0], y[0], 0, 0], - [x[1], y[1], 0, 0], - [x[2], y[2], 1, 0], - [0, 0, 0, 1], - ] - ) - @ matrix - ) - intermediate[0, 3] = (intermediate[0, 3] * scale) + eastings - intermediate[1, 3] = (intermediate[1, 3] * scale) + northings - intermediate[2, 3] = (intermediate[2, 3] * scale) + orthogonal_height - return intermediate - - -def local2global_ifc4x3( - matrix, - eastings, - northings, - orthogonal_height, - x_axis_abscissa, - x_axis_ordinate, - scale=1.0, - factor_x=1.0, - factor_y=1.0, - factor_z=1.0, -): - # Matrix is a 4x4 matrix typically describing the object placement of an element. theta = math.atan2(x_axis_ordinate, x_axis_abscissa) scale_and_factor_matrix = np.array( [ @@ -436,7 +405,45 @@ def local2global_ifc4x3( return result -def global2local(matrix, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None): +def auto_local2global( + ifc_file: ifcopenshell.file, matrix: npt.NDArray[np.float64], should_return_in_map_units: bool = True +) -> npt.NDArray[np.float64]: + """Convert a local matrix to a global map matrix + + The necessary georeferencing map conversion is automatically detected from + the IFC map conversion parameters present in the IFC model. If no map + conversion is present, then the matrix is returned unchanged. + + :param ifc_file: The IFC file + :param matrix: A 4x4 numpy matrix representing local coordinates. + :param should_return_in_map_units: If true, the result is given in map units. + If false, the result will be converted back into project units. + :return: A numpy 4x4 array matrix representing global coordinates. + """ + parameters = get_helmert_transformation_parameters(ifc_file) + if not parameters: + return matrix.copy() + result = local2global(matrix, *parameters) + if should_return_in_map_units: + return result + result[0][3] /= parameters.scale + result[1][3] /= parameters.scale + result[2][3] /= parameters.scale + return result + + +def global2local( + matrix: npt.NDArray[np.float64], + eastings: float = 0.0, + northings: float = 0.0, + orthogonal_height: float = 0.0, + x_axis_abscissa: float = 1.0, + x_axis_ordinate: float = 0.0, + scale: float = 1.0, + factor_x: float = 1.0, + factor_y: float = 1.0, + factor_z: float = 1.0, +) -> npt.NDArray[np.float64]: """Manually convert a 4x4 matrix from global to local coordinates This function is for advanced users as it allows you to specify your own @@ -446,76 +453,91 @@ def global2local(matrix, eastings, northings, orthogonal_height, x_axis_abscissa you are applying your own temporary false origin (such as when federating models for digital twins of large cities). - No unit conversion is performed. - :param matrix: A 4x4 numpy matrix representing global coordinates. - :type matrix: np.array :param eastings: The eastings offset to apply. - :type eastings: float :param northings: The northings offset to apply. - :type northings: float :param orthogonal_height: The orthogonal height offset to apply. - :type orthogonal_height: float :param x_axis_abscissa: The X axis abscissa (i.e. first coordinate) of the 2D vector that points to the local X axis when in map coordinates. - :type x_axis_abscissa: float :param x_axis_ordinate: The X axis ordinate (i.e. second coordinate) of the 2D vector that points to the local X axis when in map coordinates. - :type x_axis_ordinate: float :param scale: The combined scale factor to convert from local coordinates to map coordinates. - :type scale: float :return: A numpy 4x4 array matrix representing local coordinates. - :rtype: np.array """ - if scale is None: - scale = 1.0 - x = np.array([x_axis_abscissa, x_axis_ordinate, 0]) - x /= np.linalg.norm(x) - y = np.cross(np.array([0, 0, 1]), x) - result = matrix.copy() - result[0, 3] = (result[0, 3] - eastings) / scale - result[1, 3] = (result[1, 3] - northings) / scale - result[2, 3] = (result[2, 3] - orthogonal_height) / scale - return ( - np.linalg.inv( - np.array( - [ - [x[0], y[0], 0, 0], - [x[1], y[1], 0, 0], - [x[2], y[2], 1, 0], - [0, 0, 0, 1], - ] - ) - ) - @ result + theta = math.atan2(x_axis_ordinate, x_axis_abscissa) + scale_and_factor_matrix = np.array( + [ + [scale * factor_x, 0, 0, 0], + [0, scale * factor_y, 0, 0], + [0, 0, scale * factor_z, 0], + [0, 0, 0, 1], + ] ) + rotation_matrix = np.array( + [ + [math.cos(theta), -math.sin(theta), 0, 0], + [math.sin(theta), math.cos(theta), 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + ] + ) + result = matrix.copy() + result[0][3] -= eastings + result[1][3] -= northings + result[2][3] -= orthogonal_height + result = np.linalg.inv(scale_and_factor_matrix) @ np.linalg.inv(rotation_matrix) @ result + result[:, 0][0:3] /= np.linalg.norm(result[:, 0][0:3]) + result[:, 1][0:3] /= np.linalg.norm(result[:, 1][0:3]) + result[:, 2][0:3] /= np.linalg.norm(result[:, 2][0:3]) + return result -def xaxis2angle(x, y): +def auto_global2local( + ifc_file: ifcopenshell.file, matrix: npt.NDArray[np.float64], is_specified_in_map_units: bool = True +) -> npt.NDArray[np.float64]: + """Convert a global map matrix to a local matrix + + The necessary georeferencing map conversion is automatically detected from + the IFC map conversion parameters present in the IFC model. If no map + conversion is present, then the matrix is returned unchanged. + + :param ifc_file: The IFC file + :param matrix: A 4x4 numpy matrix representing local coordinates. + :param should_return_in_map_units: If true, the result is given in map units. + If false, the result will be converted back into project units. + :param is_specified_in_map_units: True if the input matrix is in map units. + :return: A numpy 4x4 array matrix representing global coordinates. + """ + parameters = get_helmert_transformation_parameters(ifc_file) + if not parameters: + return matrix.copy() + if not is_specified_in_map_units: + matrix = matrix.copy() + matrix[0][3] *= parameters.scale + matrix[1][3] *= parameters.scale + matrix[2][3] *= parameters.scale + return global2local(matrix, *parameters) + + +def xaxis2angle(x: float, y: float) -> float: """Converts X axis abscissa and ordinates to an angle in decimal degrees :param x: The X axis abscissa - :type x: float :param y: The X axis ordinate - :type y: float :return: The equivalent angle in decimal degrees from the X axis - :rtype: float """ return math.degrees(math.atan2(y, x)) * -1 -def yaxis2angle(x, y): +def yaxis2angle(x: float, y: float) -> float: """Converts Y axis abscissa and ordinates to an angle in decimal degrees The Y axis abscissa and ordinate is how IFC stores true north. :param x: The Y axis abscissa - :type x: float :param y: The Y axis ordinate - :type y: float :return: The equivalent angle in decimal degrees from the Y axis - :rtype: float """ angle = math.degrees(math.atan2(y, x)) - 90 if angle < -180: @@ -525,7 +547,7 @@ def yaxis2angle(x, y): return angle -def get_grid_north(ifc_file): +def get_grid_north(ifc_file: ifcopenshell.file) -> float: """Get an angle pointing to map grid north Anticlockwise is positive. @@ -539,31 +561,15 @@ def get_grid_north(ifc_file): https://www.buildingsmart.org/standards/bsi-standards/standards-library/ :param ifc_file: The IFC file - :type ifc_file: ifcopenshell.file :return: An angle to grid north in decimal degrees - :rtype: float """ - conversion = None - try: - conversion = ifc_file.by_type("IfcMapConversion")[0] - except: - pass - if conversion: - if not conversion.XAxisAbscissa or not conversion.XAxisOrdinate: - return 0 - xaa = conversion.XAxisAbscissa - xao = conversion.XAxisOrdinate - else: - project = ifc_file.by_type("IfcProject")[0] - conversion = ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion") - if not conversion: - return 0 - xaa = conversion.get("XAxisAbscissa", None) or 0 - xao = conversion.get("XAxisOrdinate", None) or 0 - return xaxis2angle(xaa, xao) + parameters = get_helmert_transformation_parameters(ifc_file) + if not parameters: + return 0 + return xaxis2angle(parameters.xaa, parameters.xao) -def get_true_north(ifc_file): +def get_true_north(ifc_file: ifcopenshell.file) -> float: """Get an angle pointing to global true north Anticlockwise is positive. @@ -571,13 +577,11 @@ def get_true_north(ifc_file): Always remember that true north is not a constant! (Unless you are working in polar coordinates) This true north is only a reference value useful for things like solar analysis on small sites (<1km). If you're after the north - that your surveyor is using, you're probably after ``get_grid_north`` + that your surveyor is using, you're probably after :func:`get_grid_north` instead. :param ifc_file: The IFC file - :type ifc_file: ifcopenshell.file :return: An angle to true north in decimal degrees - :rtype: float """ try: for context in ifc_file.by_type("IfcGeometricRepresentationContext", include_subtypes=False): @@ -588,15 +592,13 @@ def get_true_north(ifc_file): return 0 -def angle2xaxis(angle): +def angle2xaxis(angle: float) -> tuple[float, float]: """Converts an angle into an X axis abscissa and ordinate - The inverse of ``xaxis2angle``. + The inverse of :func:`xaxis2angle`. :param angle: The angle in decimal degrees where anticlockwise is positive. - :type angle: float :return: A tuple of X axis abscissa and ordinate - :rtype: tuple[float] """ angle_rad = math.radians(angle) x = math.cos(angle_rad) @@ -604,15 +606,13 @@ def angle2xaxis(angle): return x, y -def angle2yaxis(angle): +def angle2yaxis(angle: float) -> tuple[float, float]: """Converts an angle into an Y axis abscissa and ordinate - The inverse of ``yaxis2angle``. + The inverse of :func:`yaxis2angle`. :param angle: The angle in decimal degrees where anticlockwise is positive. - :type angle: float :return: A tuple of Y axis abscissa and ordinate - :rtype: tuple[float] """ angle_rad = math.radians(angle) x = -math.sin(angle_rad) diff --git a/src/ifcopenshell-python/test/util/test_geolocation.py b/src/ifcopenshell-python/test/util/test_geolocation.py index db763b7c6e..76b2b339c8 100644 --- a/src/ifcopenshell-python/test/util/test_geolocation.py +++ b/src/ifcopenshell-python/test/util/test_geolocation.py @@ -16,7 +16,6 @@ # You should have received a copy of the GNU Lesser General Public License # along with IfcOpenShell. If not, see . -import pytest import numpy as np import test.bootstrap import ifcopenshell.api.root @@ -138,7 +137,7 @@ class TestAutoZ2E(test.bootstrap.IFC4): class TestLocal2Global(test.bootstrap.IFC4): - def test_converting_from_a_local_matrix_to_a_global_matrix(self): + def test_run(self): m = np.eye(4) m2 = np.eye(4) assert np.allclose(subject.local2global(m, 0, 0, 0, 1.0, 0.0), m2) @@ -155,42 +154,121 @@ class TestLocal2Global(test.bootstrap.IFC4): m2[:, 3][0:3] = [2, 3, 3] assert np.allclose(subject.local2global(m, 1, 2, 3, 1.0, 0.0), m2) - m2[:, 3][0:3] = [3, 4, 3] - assert np.allclose(subject.local2global(m, 1, 2, 3, 1.0, 0.0, 2), m2) + m[:, 3][0:3] = [1000, 1000, 0] + m2[:, 3][0:3] = [2, 3, 3] + assert np.allclose(subject.local2global(m, 1, 2, 3, 1.0, 0.0, 0.001), m2) + m[:, 3][0:3] = [1, 1, 0] m2[:, 0][0:3] = [0, 1, 0] m2[:, 1][0:3] = [-1, 0, 0] m2[:, 3][0:3] = [0, 3, 3] assert np.allclose(subject.local2global(m, 1, 2, 3, 0.0, 1.0), m2) + m[:, 3][0:3] = [1, 1, 1] + m2 = np.eye(4) + m2[:, 3][0:3] = [5, 8, 11] + assert np.allclose(subject.local2global(m, 1, 2, 3, 1.0, 0.0, 2, 2, 3, 4), m2) -class TestLocal2GlobalIfc4X3(test.bootstrap.IFC4): - def test_converting_from_a_local_matrix_to_a_global_matrix(self): + +class TestGlobal2Local(test.bootstrap.IFC4): + def test_run(self): m = np.eye(4) m2 = np.eye(4) - assert np.allclose(subject.local2global_ifc4x3(m, 0, 0, 0, 1.0, 0.0), m2) + assert np.allclose(subject.global2local(m2, 0, 0, 0, 1.0, 0.0), m) m2[:, 3][0:3] = [1, 2, 3] - assert np.allclose(subject.local2global_ifc4x3(m, 1, 2, 3, 1.0, 0.0), m2) + assert np.allclose(subject.global2local(m2, 1, 2, 3, 1.0, 0.0), m) m2[:, 0][0:3] = [0, 1, 0] m2[:, 1][0:3] = [-1, 0, 0] - assert np.allclose(subject.local2global_ifc4x3(m, 1, 2, 3, 0.0, 1.0), m2) + assert np.allclose(subject.global2local(m2, 1, 2, 3, 0.0, 1.0), m) m[:, 3][0:3] = [1, 1, 0] m2 = np.eye(4) m2[:, 3][0:3] = [2, 3, 3] - assert np.allclose(subject.local2global_ifc4x3(m, 1, 2, 3, 1.0, 0.0), m2) + assert np.allclose(subject.global2local(m2, 1, 2, 3, 1.0, 0.0), m) - m2[:, 3][0:3] = [3, 4, 3] - assert np.allclose(subject.local2global_ifc4x3(m, 1, 2, 3, 1.0, 0.0, 2), m2) + m[:, 3][0:3] = [1000, 1000, 0] + m2[:, 3][0:3] = [2, 3, 3] + assert np.allclose(subject.global2local(m2, 1, 2, 3, 1.0, 0.0, 0.001), m) + m[:, 3][0:3] = [1, 1, 0] m2[:, 0][0:3] = [0, 1, 0] m2[:, 1][0:3] = [-1, 0, 0] m2[:, 3][0:3] = [0, 3, 3] - assert np.allclose(subject.local2global_ifc4x3(m, 1, 2, 3, 0.0, 1.0), m2) + assert np.allclose(subject.global2local(m2, 1, 2, 3, 0.0, 1.0), m) m[:, 3][0:3] = [1, 1, 1] m2 = np.eye(4) m2[:, 3][0:3] = [5, 8, 11] - assert np.allclose(subject.local2global_ifc4x3(m, 1, 2, 3, 1.0, 0.0, 2, 2, 3, 4), m2) + assert np.allclose(subject.global2local(m2, 1, 2, 3, 1.0, 0.0, 2, 2, 3, 4), m) + + +class TestAutoLocal2Global(test.bootstrap.IFC4): + def test_no_georeferencing(self): + m = np.eye(4) + ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcProject") + assert np.allclose(subject.auto_local2global(self.file, m), m) + m[:, 3][0:3] = [1, 2, 3] + assert np.allclose(subject.auto_local2global(self.file, m), m) + + def test_map_conversion(self): + m = np.eye(4) + m2 = np.eye(4) + ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcProject") + ifcopenshell.api.context.add_context(self.file, "Model") + ifcopenshell.api.georeference.add_georeferencing(self.file) + ifcopenshell.api.georeference.edit_georeferencing( + self.file, + projected_crs={"Name": "EPSG:7856"}, + map_conversion={"Eastings": 1, "Northings": 2, "OrthogonalHeight": 3}, + ) + m2[:, 3][0:3] = [1, 2, 3] + assert np.allclose(subject.auto_local2global(self.file, m), m2) + ifcopenshell.api.georeference.edit_georeferencing(self.file, map_conversion={"Scale": 0.001}) + assert np.allclose(subject.auto_local2global(self.file, m), m2) + m2[:, 3][0:3] = [1000, 2000, 3000] + assert np.allclose(subject.auto_local2global(self.file, m, should_return_in_map_units=False), m2) + ifcopenshell.api.georeference.edit_georeferencing( + self.file, map_conversion={"XAxisAbscissa": 0, "XAxisOrdinate": 1} + ) + m[:, 3][0:3] = [1000, 1000, 0] + m2[:, 0][0:3] = [0, 1, 0] + m2[:, 1][0:3] = [-1, 0, 0] + m2[:, 3][0:3] = [0, 3, 3] + assert np.allclose(subject.auto_local2global(self.file, m), m2) + + +class TestAutoGlobal2Local(test.bootstrap.IFC4): + def test_no_georeferencing(self): + m = np.eye(4) + ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcProject") + assert np.allclose(subject.auto_global2local(self.file, m), m) + m[:, 3][0:3] = [1, 2, 3] + assert np.allclose(subject.auto_global2local(self.file, m), m) + + def test_map_conversion(self): + m = np.eye(4) + m2 = np.eye(4) + ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcProject") + ifcopenshell.api.context.add_context(self.file, "Model") + ifcopenshell.api.georeference.add_georeferencing(self.file) + ifcopenshell.api.georeference.edit_georeferencing( + self.file, + projected_crs={"Name": "EPSG:7856"}, + map_conversion={"Eastings": 1, "Northings": 2, "OrthogonalHeight": 3}, + ) + m2[:, 3][0:3] = [1, 2, 3] + assert np.allclose(subject.auto_global2local(self.file, m2), m) + ifcopenshell.api.georeference.edit_georeferencing(self.file, map_conversion={"Scale": 0.001}) + assert np.allclose(subject.auto_global2local(self.file, m2), m) + m2[:, 3][0:3] = [1000, 2000, 3000] + assert np.allclose(subject.auto_global2local(self.file, m2, is_specified_in_map_units=False), m) + ifcopenshell.api.georeference.edit_georeferencing( + self.file, map_conversion={"XAxisAbscissa": 0, "XAxisOrdinate": 1} + ) + m[:, 3][0:3] = [1000, 1000, 0] + m2[:, 0][0:3] = [0, 1, 0] + m2[:, 1][0:3] = [-1, 0, 0] + m2[:, 3][0:3] = [0, 3, 3] + assert np.allclose(subject.auto_global2local(self.file, m2), m)