IfcOpenShell/src/ifcopenshell-python/ifcopenshell/util/geolocation.py

# IfcOpenShell - IFC toolkit and geometry engine
# Copyright (C) 2021 Dion Moult <dion@thinkmoult.com>
#
# This file is part of IfcOpenShell.
#
# IfcOpenShell is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# IfcOpenShell is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.

import math
from decimal import ROUND_HALF_UP, Decimal
from typing import NamedTuple, Optional, Union, Any

import numpy as np

import ifcopenshell
import ifcopenshell.util.element
import ifcopenshell.util.placement

MatrixType = ifcopenshell.util.placement.MatrixType


class HelmertTransformation(NamedTuple):
    e: float
    n: float
    h: float
    xaa: float
    xao: float
    scale: float
    factor_x: float
    factor_y: float
    factor_z: float


def dms2dd(degrees: int, minutes: int, seconds: int, us: int = 0) -> float:
    """Convert degrees, minutes, and (micro)seconds to decimal degrees
    All components must be either positive or negative.

    :param degrees: The degrees component
    :param minutes: The minutes component
    :param seconds: The seconds component
    :param us: The microseconds component
    :return: The angle in decimal degrees.
    """
    all_positive_or_zero = degrees >= 0 and minutes >= 0 and seconds >= 0 and us >= 0
    all_negative_or_zero = degrees <= 0 and minutes <= 0 and seconds <= 0 and us <= 0
    assert all_positive_or_zero or all_negative_or_zero
    return degrees + minutes / 60.0 + seconds / 3600.0 + us / 3600000000.0


def dd2dms(dd: float, use_us: bool = False) -> Union[tuple[int, int, int, int], tuple[int, int, float]]:
    """Convert decimal degrees to degrees, minutes, and (micro)seconds format

    :param dd: The decimal degrees
    :param use_us: True if to include microseconds and false otherwise. Defaults to false.
    :return: The angle in a tuple of either 3 or 4 values,
        4 values: integer number of degrees, integer number of minutes, integer number of seconds and integer number of microseconds
        3 values: integer number of degrees, integer number of minutes, and a float number for seconds
    :note: the tuple follows the format of IfcCompoundPlaneAngleMeasure. Namely all of its components are either positive or negative.
    """
    dd_decimal = Decimal(str(dd))

    degrees = int(dd_decimal)
    degrees_decimal = Decimal(degrees)

    fractional_part = dd_decimal - degrees_decimal

    minutes_decimal = fractional_part * Decimal(60)
    minutes = int(minutes_decimal)
    minutes_decimal_int = Decimal(minutes)

    seconds_decimal = (minutes_decimal - minutes_decimal_int) * Decimal(60)

    if use_us:
        seconds = int(seconds_decimal)
        seconds_decimal_int = Decimal(seconds)
        microseconds_decimal = (seconds_decimal - seconds_decimal_int) * Decimal(1000000)
        microseconds = int(microseconds_decimal.quantize(Decimal(1), rounding=ROUND_HALF_UP))
        return (degrees, minutes, seconds, microseconds)
    else:
        seconds_float = float(seconds_decimal)
        return (degrees, minutes, seconds_float)


def xyz2enh(
    x: float,
    y: float,
    z: float,
    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,
) -> tuple[float, float, float]:
    """Manually convert local XYZ coordinates to map eastings, northings, and height

    This function is for advanced users as it allows you to specify your own
    helmert transformation parameters (i.e. those typically stored in
    IfcMapConversion). This manual approach is useful for tests or in case your
    are setting your helmert transformations in non-standard locations, or if
    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 :func:`auto_xyz2enh` instead.

    :param x: The X local engineering coordinate.
    :param y: The Y local engineering coordinate.
    :param z: The Z local engineering coordinate.
    :param eastings: The eastings offset to apply.
    :param northings: The northings offset to apply.
    :param orthogonal_height: The orthogonal height offset to apply.
    :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.
    :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.
    :param scale: The unit scale such that local ordinate * scale = map
        ordinate. E.g. if your project is in millimeters but your CRS is in
        meters, your scale should be 0.001.
    :param factor_x: The combined scale factor for the X value to convert from
        local coordinates to map coordinates. Your surveyor will typically know
        this number and approximate it as a constant on a small site. Typically
        factor_x and factor_y will be identical, and factor_z will be 1.
    :param factor_y: Same but for the Y value.
    :param factor_z: Same but for the Z value.
    :return: A tuple of three ordinates representing the easting, northing and height.
    """
    theta = math.atan2(x_axis_ordinate, x_axis_abscissa)
    eastings = (scale * factor_x * math.cos(theta) * x) - (scale * factor_y * math.sin(theta) * y) + eastings
    northings = (scale * factor_x * math.sin(theta) * x) + (scale * factor_y * math.cos(theta) * y) + northings
    height = (scale * factor_z * z) + orthogonal_height
    return (eastings, northings, height)


def auto_xyz2enh(
    ifc_file: ifcopenshell.file, x: float, y: float, z: float, should_return_in_map_units: bool = True
) -> tuple[float, float, float]:
    """Convert from local XYZ coordinates to global map coordinate eastings, northings, and heights

    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 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":
    https://www.buildingsmart.org/standards/bsi-standards/standards-library/

    :param ifc_file: The IFC file
    :param x: The X local engineering coordinate provided in project length units.
    :param y: The Y local engineering coordinate provided in project length units.
    :param z: The Z local engineering coordinate provided in project length units.
    :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.
    """
    parameters = get_helmert_transformation_parameters(ifc_file)
    if not parameters:
        return x, y, z
    wcs = get_wcs(ifc_file)
    if wcs is not None:
        x, y, z = (np.linalg.inv(wcs) @ np.array((x, y, z, 1)))[:3]
    enh = xyz2enh(x, y, z, *parameters)
    if should_return_in_map_units:
        return enh
    return enh[0] / parameters.scale, enh[1] / parameters.scale, enh[2] / parameters.scale


def auto_enh2xyz(
    ifc_file: ifcopenshell.file, easting: float, northing: float, height: float, is_specified_in_map_units: bool = True
) -> tuple[float, float, float]:
    """Convert from global map coordinate eastings, northings, and heights to local XYZ coordinates

    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.

    For IFC2X3, the map conversion is detected from the IfcProject's
    ePSet_MapConversion. See the "User Guide for Geo-referencing in IFC":
    https://www.buildingsmart.org/standards/bsi-standards/standards-library/

    :param ifc_file: The IFC file
    :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)
    if not parameters:
        return easting, northing, height
    if not is_specified_in_map_units:
        easting *= parameters.scale
        northing *= parameters.scale
        height *= parameters.scale
    xyz = enh2xyz(easting, northing, height, *parameters)
    wcs = get_wcs(ifc_file)
    if wcs is not None:
        xyz = tuple((wcs @ np.array((*xyz, 1)))[:3])
    return xyz


def get_helmert_transformation_parameters(ifc_file: ifcopenshell.file) -> Optional[HelmertTransformation]:
    """Retrieves the parameters of a helmert transformation that represents a
    coordinate operation

    This coordinate operation is typically what is used to convert between
    local engineering coordinates and map coordinates.

    :param ifc_file: The IFC model, typically containing an
        IfcCoordinateOperation such as an IfcMapConversion.
    :return: The parameters of the transformation.
    """
    if ifc_file.schema == "IFC2X3":
        project = ifc_file.by_type("IfcProject")[0]
        conversion = ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion")
        if not conversion:
            return
        e = conversion.get("Eastings", None) or 0
        n = conversion.get("Northings", None) or 0
        h = conversion.get("OrthogonalHeight", None) or 0
        xaa = conversion.get("XAxisAbscissa", None) or 0
        xao = conversion.get("XAxisOrdinate", None) or 0
        scale = conversion.get("Scale", None) or 1
        factor_x = factor_y = factor_z = 1
    else:
        conversion = ifc_file.by_type("IfcCoordinateOperation")
        if not conversion:
            return
        conversion = conversion[0]

        if conversion.is_a("IfcMapConversion"):
            e = conversion.Eastings or 0
            n = conversion.Northings or 0
            h = conversion.OrthogonalHeight or 0
            xaa = conversion.XAxisAbscissa or 0
            xao = conversion.XAxisOrdinate or 0
            scale = conversion.Scale or 1
            if conversion.is_a() == "IfcMapConversionScaled":
                factor_x = conversion.FactorX
                factor_y = conversion.FactorY
                factor_z = conversion.FactorZ
            else:
                factor_x = factor_y = factor_z = 1
        elif conversion.is_a() == "IfcRigidOperation":
            e = conversion.FirstCoordinate.wrappedValue
            n = conversion.SecondCoordinate.wrappedValue
            h = conversion.Height or 0
            xaa = 1.0
            xao = 0.0
            scale = factor_x = factor_y = factor_z = 1

    if not xaa and not xao:
        xaa = 1.0
        xao = 0.0

    return HelmertTransformation(e, n, h, xaa, xao, scale, factor_x, factor_y, factor_z)


def get_crs(ifc_file: ifcopenshell.file) -> dict[str, Any]:
    """Get CRS information from an IFC file."""
    if ifc_file.schema == "IFC2X3":
        return ifcopenshell.util.element.get_pset(ifc_file.by_type("IfcProject")[0], "ePSet_ProjectedCRS")
    for context in ifc_file.by_type("IfcGeometricRepresentationContext", include_subtypes=False):
        if operation := context.HasCoordinateOperation:
            return operation[0].TargetCRS.get_info()


def auto_z2e(ifc_file: ifcopenshell.file, z: float, should_return_in_map_units: bool = True) -> float:
    """Convert a Z coordinate to an elevation using model georeferencing data

    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.

    For IFC2X3, the map conversion is detected from the IfcProject's
    ePSet_MapConversion. See the "User Guide for Geo-referencing in IFC":
    https://www.buildingsmart.org/standards/bsi-standards/standards-library/

    :param ifc_file: The IFC file
    :param z: The Z local engineering coordinate provided in project length units.
    :return: The elevation in project length units.
    """
    parameters = get_helmert_transformation_parameters(ifc_file)
    if not parameters:
        return z
    e = z2e(z, parameters.h, parameters.scale, parameters.factor_z)
    if should_return_in_map_units:
        return e
    return e / parameters.scale


def z2e(z: float, orthogonal_height: float = 0.0, scale: float = 1.0, factor_z: float = 1.0) -> float:
    """Manually convert a Z coordinate to a map elevation

    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 :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.
    :param scale: The unit scale such that local ordinate * scale = map
        ordinate. E.g. if your project is in millimeters but your CRS is in
        meters, your scale should be 0.001.
    :param factor_x: The combined scale factor for the Z value to convert from
        local coordinates to map coordinates. Your surveyor will typically know
        this number and approximate it as a constant on a small site. This is
        typically just 1.0, as average combined scale factors usually only
        affect the XY axes.
    :return: The elevation in map units.
    """
    return (scale * factor_z * z) + orthogonal_height


def enh2xyz(
    e: float,
    n: float,
    h: float,
    eastings: float = 0.0,
    northings: float = 0.0,
    orthogonal_height: float = 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,
) -> tuple[float, float, float]:
    """Manually convert map eastings, northings, and height to local XYZ coordinates

    This function is for advanced users as it allows you to specify your own
    helmert transformation parameters (i.e. those typically stored in
    IfcMapConversion). This manual approach is useful for tests or in case your
    are setting your helmert transformations in non-standard locations, or if
    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 :func:`auto_enh2xyz` instead.

    :param e: The global easting map coordinate.
    :param n: The global northing map coordinate.
    :param h: The global height map coordinate.
    :param eastings: The eastings offset to apply.
    :param northings: The northings offset to apply.
    :param orthogonal_height: The orthogonal height offset to apply.
    :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.
    :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.
    :param scale: The unit scale such that local ordinate * scale = map
        ordinate. E.g. if your project is in millimeters but your CRS is in
        meters, your scale should be 0.001.
    :param factor_x: The combined scale factor for the X value to convert from
        local coordinates to map coordinates. Your surveyor will typically know
        this number and approximate it as a constant on a small site. Typically
        factor_x and factor_y will be identical, and factor_z will be 1.
    :param factor_y: Same but for the Y value.
    :param factor_z: Same but for the Z value.
    :return: A tuple of three ordinates representing XYZ.
    """
    theta = math.atan2(x_axis_ordinate, x_axis_abscissa)
    sint = math.sin(theta)
    cost = math.cos(theta)
    x = (((e - eastings) * cost) + ((n - northings) * sint)) / (scale * factor_x)
    y = (((eastings - e) * sint) + ((n - northings) * cost)) / (scale * factor_y)
    z = ((h - orthogonal_height) / scale) / factor_z
    return (x, y, z)


def local2global(
    matrix: MatrixType,
    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,
) -> MatrixType:
    """Manually convert a 4x4 matrix from local to global coordinates

    This function is for advanced users as it allows you to specify your own
    helmert transformation parameters (i.e. those typically stored in
    IfcMapConversion). This manual approach is useful for tests or in case your
    are setting your helmert transformations in non-standard locations, or if
    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 :func:`auto_local2global` instead.

    :param matrix: A 4x4 numpy matrix representing local coordinates.
    :param eastings: The eastings offset to apply.
    :param northings: The northings offset to apply.
    :param orthogonal_height: The orthogonal height offset to apply.
    :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.
    :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.
    :param scale: The combined scale factor to convert from local coordinates
        to map coordinates.
    :return: A numpy 4x4 array matrix representing global coordinates.
    """
    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 = rotation_matrix @ scale_and_factor_matrix @ matrix
    result[:3, 0] /= np.linalg.norm(result[:3, 0])
    result[:3, 1] /= np.linalg.norm(result[:3, 1])
    result[:3, 2] /= np.linalg.norm(result[:3, 2])
    result[0, 3] += eastings
    result[1, 3] += northings
    result[2, 3] += orthogonal_height
    return result


def auto_local2global(
    ifc_file: ifcopenshell.file, matrix: MatrixType, should_return_in_map_units: bool = True
) -> MatrixType:
    """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()
    wcs = get_wcs(ifc_file)
    if wcs is not None:
        matrix = np.linalg.inv(wcs) @ matrix
    result = local2global(matrix, *parameters)
    if should_return_in_map_units:
        return result
    result[:3, 3] /= parameters.scale
    return result


def global2local(
    matrix: MatrixType,
    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,
) -> MatrixType:
    """Manually convert a 4x4 matrix from global to local coordinates

    This function is for advanced users as it allows you to specify your own
    helmert transformation parameters (i.e. those typically stored in
    IfcMapConversion). This manual approach is useful for tests or in case your
    are setting your helmert transformations in non-standard locations, or if
    you are applying your own temporary false origin (such as when federating
    models for digital twins of large cities).

    :param matrix: A 4x4 numpy matrix representing global coordinates.
    :param eastings: The eastings offset to apply.
    :param northings: The northings offset to apply.
    :param orthogonal_height: The orthogonal height offset to apply.
    :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.
    :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.
    :param scale: The combined scale factor to convert from local coordinates
        to map coordinates.
    :return: A numpy 4x4 array matrix representing local coordinates.
    """
    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[:3, 0] /= np.linalg.norm(result[:3, 0])
    result[:3, 1] /= np.linalg.norm(result[:3, 1])
    result[:3, 2] /= np.linalg.norm(result[:3, 2])
    return result


def auto_global2local(
    ifc_file: ifcopenshell.file, matrix: MatrixType, is_specified_in_map_units: bool = True
) -> MatrixType:
    """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[:3, 3] *= parameters.scale
    result = global2local(matrix, *parameters)
    wcs = get_wcs(ifc_file)
    if wcs is not None:
        return wcs @ result
    return result


def xaxis2angle(x: float, y: float) -> float:
    """Converts X axis abscissa and ordinates to an angle in decimal degrees

    The X axis abscissa and ordinate is how IFC stores grid north.

    This X axis vector indicates "where is project east, if grid north is up
    the page?". See the diagram on the IfcGeometricRepresentationContext
    documentation for clarification.

    The angle indicates "how do I rotate project east to get to grid east?".
    Alternatively: "how do I rotate project north to get to grid north?".
    Positive angles are anticlockwise.

    :param x: The X axis abscissa
    :param y: The X axis ordinate
    :return: The equivalent angle in decimal degrees from the X axis
    """
    return math.degrees(math.atan2(y, x)) * -1


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.

    This Y axis vector indicates "where is true north, if project north is up
    the page?". See the diagram on the IfcGeometricRepresentationContext
    documentation for clarification.

    The angle indicates "how do I rotate project north to get to true north?".
    Positive angles are anticlockwise.

    :param x: The Y axis abscissa
    :param y: The Y axis ordinate
    :return: The equivalent angle in decimal degrees from the Y axis
    """
    angle = math.degrees(math.atan2(y, x)) - 90
    if angle < -180:
        angle += 360
    elif angle > 180:
        angle -= 360
    return angle


def get_grid_north(ifc_file: ifcopenshell.file) -> float:
    """Get an angle pointing to map grid north

    Anticlockwise is positive.

    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.

    For IFC2X3, the map conversion is detected from the IfcProject's
    ePSet_MapConversion. See the "User Guide for Geo-referencing in IFC":
    https://www.buildingsmart.org/standards/bsi-standards/standards-library/

    :param ifc_file: The IFC file
    :return: An angle to grid north in decimal degrees
    """
    parameters = get_helmert_transformation_parameters(ifc_file)
    if not parameters:
        return 0
    return xaxis2angle(parameters.xaa, parameters.xao)


def get_true_north(ifc_file: ifcopenshell.file) -> float:
    """Get an angle pointing to global true north

    Anticlockwise is positive.

    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 :func:`get_grid_north`
    instead.

    :param ifc_file: The IFC file
    :return: An angle to true north in decimal degrees
    """
    try:
        for context in ifc_file.by_type("IfcGeometricRepresentationContext", include_subtypes=False):
            if context.TrueNorth:
                return yaxis2angle(*context.TrueNorth.DirectionRatios[0:2])
    except:
        return 0
    return 0


def angle2xaxis(angle: float) -> tuple[float, float]:
    """Converts an angle into an X axis abscissa and ordinate

    The inverse of :func:`xaxis2angle`.

    :param angle: The angle in decimal degrees where anticlockwise is positive.
    :return: A tuple of X axis abscissa and ordinate
    """
    angle_rad = math.radians(angle)
    x = math.cos(angle_rad)
    y = -math.sin(angle_rad)
    return x, y


def angle2yaxis(angle: float) -> tuple[float, float]:
    """Converts an angle into an Y axis abscissa and ordinate

    The inverse of :func:`yaxis2angle`.

    :param angle: The angle in decimal degrees where anticlockwise is positive.
    :return: A tuple of Y axis abscissa and ordinate
    """
    angle_rad = math.radians(angle)
    x = -math.sin(angle_rad)
    y = math.cos(angle_rad)
    return x, y


def get_wcs(ifc_file: ifcopenshell.file) -> Optional[MatrixType]:
    """Gets the WCS (prioritising 3D contexts) as a matrix

    :param: The IFC file
    :return: A 4x4 matrix in project units
    """
    wcs = None
    for context in ifc_file.by_type("IfcGeometricRepresentationContext", include_subtypes=False):
        wcs = context.WorldCoordinateSystem
        if context.ContextType == "Model":
            break
    if wcs:
        return ifcopenshell.util.placement.get_axis2placement(wcs)