Implement WCS consideration in XYZ <-> ENH conversions.

src/ifcopenshell-python/ifcopenshell/api/georeference/edit_georeferencing.py

@@ -81,28 +81,32 @@ def edit_georeferencing(
         if not (project := file.by_type("IfcProject")):
             return
         project = project[0]
-        if crs := ifcopenshell.util.element.get_pset(project, "ePSet_ProjectedCRS"):
+        if projected_crs:
-            crs = file.by_id(crs["id"])
+            if crs := ifcopenshell.util.element.get_pset(project, "ePSet_ProjectedCRS"):
-            for k, v in projected_crs.items():
+                crs = file.by_id(crs["id"])
-                if k == "Description":
+                for k, v in projected_crs.items():
-                    v = file.createIfcText(v)
+                    if k == "Description":
-                elif k == "Name":
+                        v = file.createIfcText(v)
-                    v = file.createIfcLabel(v)
+                    elif k == "Name":
-                else:
+                        v = file.createIfcLabel(v)
-                    v = file.createIfcIdentifier(v)
+                    else:
-            ifcopenshell.api.pset.edit_pset(file, crs, properties=projected_crs)
+                        v = file.createIfcIdentifier(v)
-        if conversion := ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion"):
+                ifcopenshell.api.pset.edit_pset(file, crs, properties=projected_crs)
-            conversion = file.by_id(conversion["id"])
+        if coordinate_operation:
-            for k, v in coordinate_operation.items():
+            if conversion := ifcopenshell.util.element.get_pset(project, "ePSet_MapConversion"):
-                if k in ("XAxisAbscissa", "XAxisOrdinate", "Scale"):
+                conversion = file.by_id(conversion["id"])
-                    v = file.createIfcReal(v)
+                for k, v in coordinate_operation.items():
-                else:
+                    if k in ("XAxisAbscissa", "XAxisOrdinate", "Scale"):
-                    v = file.createIfcLengthMeasure(v)
+                        v = file.createIfcReal(v)
-            ifcopenshell.api.pset.edit_pset(file, conversion, properties=coordinate_operation)
+                    else:
+                        v = file.createIfcLengthMeasure(v)
+                ifcopenshell.api.pset.edit_pset(file, conversion, properties=coordinate_operation)
         return
-    conversion = file.by_type("IfcCoordinateOperation")[0]
+    if projected_crs:
-    crs = file.by_type("IfcProjectedCRS")[0]
+        crs = file.by_type("IfcProjectedCRS")[0]
-    for name, value in coordinate_operation.items():
+        for name, value in projected_crs.items():
-        setattr(conversion, name, value)
+            setattr(crs, name, value)
-    for name, value in projected_crs.items():
+    if coordinate_operation:
-        setattr(crs, name, value)
+        conversion = file.by_type("IfcCoordinateOperation")[0]
+        for name, value in coordinate_operation.items():
+            setattr(conversion, name, value)

src/ifcopenshell-python/ifcopenshell/api/georeference/edit_true_north.py

@@ -50,26 +50,26 @@ def edit_true_north(file: ifcopenshell.file, true_north: Optional[Union[tuple[fl
         # This unsets true north
         ifcopenshell.api.run("georeference.edit_true_north", model, true_north=None)
     """
-    if not true_north:
+    if isinstance(true_north, (float, int)):
-        return
-
-    if true_north is None:
-        pass
-    elif isinstance(true_north, (float, int)):
         x, y = ifcopenshell.util.geolocation.angle2yaxis(true_north)
-    else:
+    elif true_north is not None:
         x, y = true_north
 
     for context in file.by_type("IfcGeometricRepresentationContext", include_subtypes=False):
+        if context.TrueNorth and true_north is None:
+            old_true_north =  context.TrueNorth
+            context.TrueNorth = None
+            if not file.get_total_inverses(old_true_north):
+                ifcopenshell.util.element.remove_deep2(file, old_true_north)
+            continue
+
         if context.TrueNorth:
             if file.get_total_inverses(context.TrueNorth) != 1:
                 context.TrueNorth = file.create_entity("IfcDirection")
         else:
             context.TrueNorth = file.create_entity("IfcDirection")
         direction = context.TrueNorth
-        if true_north is None:
+        if context.CoordinateSpaceDimension == 2:
-            context.TrueNorth = None
-        elif context.CoordinateSpaceDimension == 2:
             direction.DirectionRatios = (x, y)
         else:
             direction.DirectionRatios = (x, y, 0.0)

src/ifcopenshell-python/ifcopenshell/api/georeference/edit_wcs.py

@@ -33,9 +33,19 @@ def edit_wcs(
 ) -> None:
     """Edits the WCS for all geometric contexts to a translation and rotation
 
-    It's recommended to leave the WCS at 0,0,0. You should generally not be
+    Typically, a project's local engineering origin (0, 0, 0) has a coordinate
-    setting it to any value. Instead, your project's origin should use a
+    operation (e.g. map conversion) to a projected CRS. If a WCS is provided,
-    coordinate operation to convert to the projected CRS.
+    the coordinate operation is relative to the WCS, not the local engineering
+    origin.
+
+    For example, if I have an IfcSite with a placement at (10, 0, 0) and a map
+    conversion of (50, 0, 0), my IfcSite's local XYZ is at (10, 0, 0) with an
+    ENH (Easting, Northing, Height) of (60, 0, 0).
+
+    If I then define by WCS at (15, 0, 0), my IfcSite's local XYZ is still at
+    (10, 0, 0) but its ENH is now at (45, 0, 0).
+
+    It's recommended to leave the WCS at 0,0,0. Please :)
 
     :param x: The X translation of the WCS
     :param y: The Y translation of the WCS

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

@@ -150,6 +150,9 @@ def auto_xyz2enh(
     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
@@ -181,7 +184,11 @@ def auto_enh2xyz(ifc_file, easting, northing, height, is_specified_in_map_units:
         easting *= parameters.scale
         northing *= parameters.scale
         height *= parameters.scale
-    return enh2xyz(easting, northing, height, *parameters)
+    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]:
@@ -423,6 +430,9 @@ def auto_local2global(
     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
@@ -517,7 +527,11 @@ def auto_global2local(
         matrix[0][3] *= parameters.scale
         matrix[1][3] *= parameters.scale
         matrix[2][3] *= parameters.scale
-    return global2local(matrix, *parameters)
+    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:

src/ifcopenshell-python/test/api/georeference/test_edit_georeferencing.py

@@ -39,6 +39,11 @@ class TestEditGeoreferencing(test.bootstrap.IFC4):
         assert conversion.Eastings == 123.45
         assert conversion.Northings == 234.56
 
+        ifcopenshell.api.georeference.edit_georeferencing(self.file, projected_crs={"Name": "EPSG:1234"})
+        assert crs.Name == "EPSG:1234"
+        ifcopenshell.api.georeference.edit_georeferencing(self.file, coordinate_operation={"Eastings": 42})
+        assert conversion.Eastings == 42
+
 
 class TestEditGeoreferencingIFC2X3(test.bootstrap.IFC2X3):
     def test_editing_georeferencing(self):

src/ifcopenshell-python/test/util/test_geolocation.py

@@ -81,6 +81,22 @@ class TestAutoXYZ2ENH(test.bootstrap.IFC4X3):
         )
         assert np.allclose(subject.auto_xyz2enh(self.file, 1000, 1000, 0), (0, 3, 3))
 
+    def test_map_conversion_with_wcs(self):
+        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"},
+            coordinate_operation={"Eastings": 1, "Northings": 2, "OrthogonalHeight": 3},
+        )
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1, y=2, z=3)
+        assert np.allclose(subject.auto_xyz2enh(self.file, 0, 0, 0), (0, 0, 0))
+        assert np.allclose(subject.auto_xyz2enh(self.file, 1, 2, 3), (1, 2, 3))
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1, y=1, z=2)
+        assert np.allclose(subject.auto_xyz2enh(self.file, 0, 0, 0), (0, 1, 1))
+        assert np.allclose(subject.auto_xyz2enh(self.file, 1, 2, 3), (1, 3, 4))
+
     def test_map_conversion_scaled(self):
         ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcProject")
         ifcopenshell.api.context.add_context(self.file, "Model")
@@ -152,6 +168,22 @@ class TestAutoENH2XYZ(test.bootstrap.IFC4X3):
         )
         assert np.allclose(subject.auto_enh2xyz(self.file, 0, 3, 3), (1000, 1000, 0))
 
+    def test_map_conversion_with_wcs(self):
+        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"},
+            coordinate_operation={"Eastings": 1, "Northings": 2, "OrthogonalHeight": 3},
+        )
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1, y=2, z=3)
+        assert np.allclose(subject.auto_enh2xyz(self.file, 0, 0, 0), (0, 0, 0))
+        assert np.allclose(subject.auto_enh2xyz(self.file, 1, 2, 3), (1, 2, 3))
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1, y=1, z=2)
+        assert np.allclose(subject.auto_enh2xyz(self.file, 0, 1, 1), (0, 0, 0))
+        assert np.allclose(subject.auto_enh2xyz(self.file, 1, 3, 4), (1, 2, 3))
+
     def test_map_conversion_scaled(self):
         ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcProject")
         ifcopenshell.api.context.add_context(self.file, "Model")
@@ -324,6 +356,25 @@ class TestAutoLocal2Global(test.bootstrap.IFC4):
         m2[:, 3][0:3] = [0, 3, 3]
         assert np.allclose(subject.auto_local2global(self.file, m), m2)
 
+    def test_map_conversion_with_wcs(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"},
+            coordinate_operation={"Eastings": 1, "Northings": 2, "OrthogonalHeight": 3},
+        )
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1, y=2, z=3)
+        assert np.allclose(subject.auto_local2global(self.file, m), m2)
+        ifcopenshell.api.georeference.edit_georeferencing(self.file, coordinate_operation={"Scale": 0.001})
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1000, y=1000, z=2000)
+        m[:, 3][0:3] = [1000, 2000, 3000]
+        m2[:, 3][0:3] = [1, 3, 4]
+        assert np.allclose(subject.auto_local2global(self.file, m), m2)
+
 
 class TestAutoGlobal2Local(test.bootstrap.IFC4):
     def test_no_georeferencing(self):
@@ -359,6 +410,25 @@ class TestAutoGlobal2Local(test.bootstrap.IFC4):
         m2[:, 3][0:3] = [0, 3, 3]
         assert np.allclose(subject.auto_global2local(self.file, m2), m)
 
+    def test_map_conversion_with_wcs(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"},
+            coordinate_operation={"Eastings": 1, "Northings": 2, "OrthogonalHeight": 3},
+        )
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1, y=2, z=3)
+        assert np.allclose(subject.auto_global2local(self.file, m2), m)
+        ifcopenshell.api.georeference.edit_georeferencing(self.file, coordinate_operation={"Scale": 0.001})
+        ifcopenshell.api.georeference.edit_wcs(self.file, x=1000, y=1000, z=2000)
+        m[:, 3][0:3] = [1000, 2000, 3000]
+        m2[:, 3][0:3] = [1, 3, 4]
+        assert np.allclose(subject.auto_global2local(self.file, m2), m)
+
 
 class TestXAxis2Angle(test.bootstrap.IFC4):
     def test_run(self):