Initial implementation of scale dependent map conversion for IFC4X3

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

@@ -54,6 +54,27 @@ def xyz2enh(x, y, z, eastings, northings, orthogonal_height, x_axis_abscissa, x_
     return (eastings, northings, height)
 
 
+def xyz2enh_ifc4x3(
+    x,
+    y,
+    z,
+    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,
+):
+    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_z2e(ifc_file, z):
     """Convert a Z coordinate to an elevation using model georeferencing data
 
@@ -78,6 +99,9 @@ def auto_z2e(ifc_file, z):
     h = conversion.OrthogonalHeight
     map_unit = conversion.TargetCRS.MapUnit
     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")
         h = ifcopenshell.util.unit.convert(
             h,
@@ -128,6 +152,46 @@ def local2global(matrix, eastings, northings, orthogonal_height, x_axis_abscissa
     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(
+        [
+            [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[:, 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])
+    result[0][3] += eastings
+    result[1][3] += northings
+    result[2][3] += orthogonal_height
+    return result
+
+
 def global2local(matrix, eastings, northings, orthogonal_height, x_axis_abscissa, x_axis_ordinate, scale=None):
     if scale is None:
         scale = 1.0
@@ -193,12 +257,13 @@ def get_true_north(ifc_file):
 def angle2xaxis(angle):
     angle_rad = math.radians(angle)
     x = math.cos(angle_rad)
-    y = - math.sin(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):
     angle_rad = math.radians(angle)
-    x = - math.sin(angle_rad)
+    x = -math.sin(angle_rad)
     y = math.cos(angle_rad)
     return x, y