Document geolocation util and add support for ePSet IFC2X3 fallbacks.

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

@@ -23,11 +23,34 @@ import ifcopenshell.util.unit
 
 
 def dms2dd(degrees, minutes, seconds, ms=0):
+    """Convert degrees, minutes, and (milli)seconds to decimal degrees
+
+    :param degrees: The degrees component
+    :type degrees: int
+    :param minutes: The minutes component
+    :type minutes: int
+    :param seconds: The seconds component
+    :type seconds: int
+    :param ms: The milliseconds component
+    :type ms: int
+    :return: The angle in decimal degrees.
+    :rtype: float
+    """
     dd = float(degrees) + float(minutes) / 60.0 + float(seconds) / (3600.0) + float(ms / 3600000000.0)
     return dd
 
 
 def dd2dms(dd, use_ms=False):
+    """Convert decimal degrees to degrees, minutes, and (milli)seconds format
+
+    :param dd: The decimal degrees
+    :type dd: float
+    :param use_ms: True if to include milliseconds and false otherwise. Defaults to false.
+    :type use_ms: bool
+    :return: The angle in a tuple of either 3 or 4 values, being degrees,
+        minutes, seconds, and optionally milliseconds.
+    :rtype: tuple[float]
+    """
     dd = float(dd)
     sign = 1 if dd >= 0 else -1
     dd = abs(dd)
@@ -43,6 +66,43 @@ def dd2dms(dd, use_ms=False):
 
 
 def xyz2enh(x, y, z, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None):
+    """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).
+
+    No unit conversion is performed.
+
+    For most scenarios you should use ``auto_xyz2enh`` instead.
+
+    :param x: The X local engineering coordinate.
+    :type x: float
+    :param y: The Y local engineering coordinate.
+    :type y: float
+    :param z: The Z local engineering coordinate.
+    :type z: float
+    :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 tuple of three ordinates representing the easting, northing and height.
+    :rtype: tuple[float]
+    """
     if scale is None:
         scale = 1.0
     rotation = math.atan2(x_axis_ordinate, x_axis_abscissa)
@@ -76,20 +136,60 @@ def xyz2enh_ifc4x3(
 
 
 def auto_xyz2enh(ifc_file, x, y, z):
+    """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 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
+    :type ifc_file: ifcopenshell.file.file
+    :param x: The X local engineering coordinate provided in project length units.
+    :type x: float
+    :param y: The Y local engineering coordinate provided in project length units.
+    :type y: float
+    :param z: The Z local engineering coordinate provided in project length units.
+    :type z: float
+    :return: The global map coordinate eastings, northings, and height in map units.
+    :rtype: tuple[float]
+    """
+    conversion = None
     try:
         conversion = ifc_file.by_type("IfcMapConversion")
     except:
-        return (x, y, z)
+        pass
-    if not conversion:
+
-        return (x, y, z)
+    if conversion:
-    conversion = conversion[0]
+        conversion = conversion[0]
-    e = conversion.Eastings or 0
+        e = conversion.Eastings or 0
-    n = conversion.Northings or 0
+        n = conversion.Northings or 0
-    h = conversion.OrthogonalHeight or 0
+        h = conversion.OrthogonalHeight or 0
-    xaa = conversion.XAxisAbscissa or 0
+        xaa = conversion.XAxisAbscissa or 0
-    xao = conversion.XAxisOrdinate or 0
+        xao = conversion.XAxisOrdinate or 0
-    scale = conversion.Scale or 0
+        scale = conversion.Scale or 1
-    map_unit = conversion.TargetCRS.MapUnit
+        map_unit = conversion.TargetCRS.MapUnit
+    else:
+        project = ifc_file.by_type("IfcProject")[0]
+        conversion = ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion")
+        if not conversion:
+            return (x, y, z)
+
+        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
+        map_unit = None
+
+    if not xaa and not xao:
+        xaa = 1.0
+        xao = 0.0
+
     if map_unit:
         # Warning! This definition has changed in IFC4X3 such that map_unit no
         # longer affects unit conversion, only the Scale attribute affects unit
@@ -103,6 +203,74 @@ def auto_xyz2enh(ifc_file, x, y, z):
     return xyz2enh(x, y, z, e, n, h, xaa, xao, scale)
 
 
+def auto_enh2xyz(ifc_file, easting, northing, height):
+    """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
+    :type ifc_file: ifcopenshell.file.file
+    :param easting: The global easting map coordinate provided in map units.
+    :type easting: float
+    :param northing: The global northing map coordinate provided in map units.
+    :type northing: float
+    :param height: The global height map coordinate provided in map units.
+    :type height: float
+    :return: The local engineering XYZ coordinates in project length units.
+    :rtype: tuple[float]
+    """
+    conversion = None
+    try:
+        conversion = ifc_file.by_type("IfcMapConversion")
+    except:
+        pass
+
+    if conversion:
+        conversion = conversion[0]
+        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
+        map_unit = conversion.TargetCRS.MapUnit
+    else:
+        project = ifc_file.by_type("IfcProject")[0]
+        conversion = ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion")
+        if not conversion:
+            return (easting, northing, height)
+
+        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
+        map_unit = None
+
+    if not xaa and not xao:
+        xaa = 1.0
+        xao = 0.0
+
+    if map_unit:
+        # Warning! This definition has changed in IFC4X3 such that map_unit no
+        # longer affects unit conversion, only the Scale attribute affects unit
+        # conversion. TODO: consolidate once IFC4X3 confirmed.
+        project_unit = ifcopenshell.util.unit.get_project_unit(ifc_file, "LENGTHUNIT")
+        map_prefix = getattr(map_unit, "Prefix", None)
+        project_prefix = getattr(project_unit, "Prefix", None)
+        e = ifcopenshell.util.unit.convert(e, map_prefix, map_unit.Name, project_prefix, project_unit.Name)
+        n = ifcopenshell.util.unit.convert(n, map_prefix, map_unit.Name, project_prefix, project_unit.Name)
+        h = ifcopenshell.util.unit.convert(h, map_prefix, map_unit.Name, project_prefix, project_unit.Name)
+    return enh2xyz(easting, northing, height, e, n, h, xaa, xao, scale)
+
+
 def auto_z2e(ifc_file, z):
     """Convert a Z coordinate to an elevation using model georeferencing data
 
@@ -110,6 +278,10 @@ def auto_z2e(ifc_file, z):
     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
     :type ifc_file: ifcopenshell.file.file
     :param z: The Z local engineering coordinate provided in project length units.
@@ -117,15 +289,25 @@ def auto_z2e(ifc_file, z):
     :return: The elevation in project length units.
     :rtype: float
     """
+    conversion = None
     try:
         conversion = ifc_file.by_type("IfcMapConversion")
     except:
-        return z
+        pass
-    if not conversion or not conversion[0].OrthogonalHeight:
+
-        return z
+    if conversion and not conversion[0].OrthogonalHeight:
-    conversion = conversion[0]
+        conversion = conversion[0]
-    h = conversion.OrthogonalHeight
+        h = conversion.OrthogonalHeight
-    map_unit = conversion.TargetCRS.MapUnit
+        map_unit = conversion.TargetCRS.MapUnit
+    else:
+        project = ifc_file.by_type("IfcProject")[0]
+        conversion = ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion")
+        if not conversion:
+            return z
+
+        h = conversion.get("OrthogonalHeight", None) or 0
+        map_unit = None
+
     if map_unit:
         # Warning! This definition has changed in IFC4X3 such that map_unit no
         # longer affects unit conversion, only the Scale attribute affects unit
@@ -142,10 +324,61 @@ def auto_z2e(ifc_file, z):
 
 
 def z2e(z, h):
+    """Manually convert a Z coordinate to an elevation
+
+    This function is for advanced users as it allows you to specify your own
+    orthogonal height offset.
+
+    For most scenarios you should use ``auto_z2e`` instead.
+
+    :param z: The Z local engineering coordinate provided in project length units.
+    :type z: float
+    :param h: The orthogonal height offset in project length units.
+    :type h: float
+    :return: The elevation in project length units.
+    :rtype: float
+    """
     return z + h
 
 
 def enh2xyz(e, n, h, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None):
+    """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).
+
+    No unit conversion is performed.
+
+    For most scenarios you should use ``auto_enh2xyz`` instead.
+
+    :param e: The global easting map coordinate.
+    :type e: float
+    :param n: The global northing map coordinate.
+    :type n: float
+    :param h: The global height map coordinate.
+    :type h: float
+    :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 tuple of three ordinates representing XYZ.
+    :rtype: tuple[float]
+    """
     if scale is None:
         scale = 1.0
     rotation = math.atan2(x_axis_ordinate, x_axis_abscissa)
@@ -158,6 +391,37 @@ def enh2xyz(e, n, h, eastings, northings, orthogonal_height, x_axis_abscissa, x_
 
 
 def local2global(matrix, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None):
+    """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).
+
+    No unit conversion is performed.
+
+    :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])
@@ -221,6 +485,37 @@ def local2global_ifc4x3(
 
 
 def global2local(matrix, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None):
+    """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).
+
+    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])
@@ -245,13 +540,31 @@ def global2local(matrix, eastings, northings, orthogonal_height, x_axis_abscissa
     )
 
 
-# Used for converting the X and Y vectors of the X Axis in IFC grid north geolocation
 def xaxis2angle(x, y):
+    """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
 
 
-# Used for converting the X and Y vectors of the Y Axis in IFC true north geolocation
 def yaxis2angle(x, y):
+    """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:
         angle += 360
@@ -261,36 +574,94 @@ def yaxis2angle(x, y):
 
 
 def get_grid_north(ifc_file):
+    """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
+    :type ifc_file: ifcopenshell.file.file
+    :return: An angle to grid north in decimal degrees
+    :rtype: float
+    """
+    conversion = None
     try:
         conversion = ifc_file.by_type("IfcMapConversion")[0]
     except:
-        return 0
+        pass
-    if not conversion.XAxisAbscissa or not conversion.XAxisOrdinate:
+    if conversion:
-        return 0
+        if not conversion.XAxisAbscissa or not conversion.XAxisOrdinate:
-    return xaxis2angle(conversion.XAxisAbscissa, 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)
 
 
 def get_true_north(ifc_file):
+    """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 ``get_grid_north``
+    instead.
+
+    :param ifc_file: The IFC file
+    :type ifc_file: ifcopenshell.file.file
+    :return: An angle to true north in decimal degrees
+    :rtype: float
+    """
     try:
         for context in ifc_file.by_type("IfcGeometricRepresentationContext", include_subtypes=False):
-            if not context.TrueNorth:
+            if context.TrueNorth:
-                continue
+                return yaxis2angle(*context.TrueNorth.DirectionRatios[0:2])
-            return yaxis2angle(*context.TrueNorth.DirectionRatios[0:2])
     except:
         return 0
     return 0
 
 
-# Used  for converting an angle in degrees to return the X and Y vectors of the X Axis in IFC grid north geolocation:
 def angle2xaxis(angle):
+    """Converts an angle into an X axis abscissa and ordinate
+
+    The inverse of ``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)
     y = -math.sin(angle_rad)
     return x, y
 
 
-# Used for converting True North angle as seen in CAD (relative to +Y)
 def angle2yaxis(angle):
+    """Converts an angle into an Y axis abscissa and ordinate
+
+    The inverse of ``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)
     y = math.cos(angle_rad)