Make all coordinates 4D

src/transformo/_typing.py

@@ -12,8 +12,10 @@
 import numpy as np
 import numpy.typing as npt
 
-Vector = Annotated[npt.NDArray[np.floating], Literal[3, 1]]
-CoordinateMatrix = Annotated[npt.NDArray[np.floating], Literal["N", 3]]
+Vector = Annotated[npt.NDArray[np.floating], Literal["N", 1]]
+Matrix = Annotated[npt.NDArray[np.floating], Literal["M", "N"]]
+CoordinateVector = Annotated[npt.NDArray[np.floating], Literal[4, 1]]
+CoordinateMatrix = Annotated[npt.NDArray[np.floating], Literal["N", 4]]
 
 ParameterValue = str | float | None
 
@@ -46,7 +48,7 @@ def from_str(  # pylint: disable=too-many-arguments
     ) -> CoordinateLike: ...
 
     @property
-    def vector(self) -> np.typing.ArrayLike: ...
+    def vector(self) -> CoordinateVector: ...
 
     @property
     def stddev(self) -> np.typing.ArrayLike: ...

src/transformo/core.py

@@ -81,7 +81,7 @@ class DataSource(pydantic.BaseModel):
 
     # station-based overrides
     overrides: dict[str, CoordinateOverrides] | None = pydantic.Field(
-        default_factory=dict
+        default_factory=dict.__call__
     )
 
     # coordinates are not included in pipeline serialization
@@ -213,6 +213,8 @@ def __post_init__(self) -> None:
                 if (t := self.overrides[key].t) is not None:
                     c.t = t
 
+            self.coordinates = sorted(self.coordinates, key=lambda c: c.station)
+
     def __add__(self, other: DataSource) -> DataSource:
         """
         Add two `DataSource`s.
@@ -331,6 +333,22 @@ def update_coordinates(self, coordinates: CoordinateMatrix) -> DataSource:
 
         return DataSource(coordinates=new_coordinates)
 
+    def station_union(self, other: DataSource) -> list[str]:
+        return list(set(self.stations) & set(other.stations))
+
+    def limit_to_stations(self, stations: list[str]) -> None:
+        """
+        Limit DataSource to stations given in supplied list.
+        """
+        stations_to_remove = []
+
+        for c in self.coordinates:
+            if c.station not in stations:
+                stations_to_remove.append(c)
+
+        for c in stations_to_remove:
+            self.coordinates.remove(c)
+
 
 class CombinedDataSource(DataSource):
     """Combination of two or more DataSource's."""

src/transformo/datatypes.py

@@ -8,7 +8,7 @@
 import pydantic
 from pydantic.dataclasses import dataclass
 
-from transformo._typing import ParameterValue
+from transformo._typing import CoordinateVector, ParameterValue
 from transformo.transformer import Transformer
 
 
@@ -64,9 +64,9 @@ def from_str(  # pylint: disable=too-many-arguments
         )
 
     @property
-    def vector(self) -> np.typing.ArrayLike:
+    def vector(self) -> CoordinateVector:
         """Coordinate given as Numpy vector (1D array)."""
-        return np.array([self.x, self.y, self.z])
+        return np.array([self.x, self.y, self.z, self.t])
 
     @property
     def stddev(self) -> np.typing.ArrayLike:

src/transformo/operators/helmert.py

@@ -12,7 +12,7 @@
 import numpy as np
 from numpy import cos, sin
 
-from transformo._typing import CoordinateMatrix, Vector
+from transformo._typing import CoordinateMatrix, Matrix, Vector
 from transformo.core import Operator
 from transformo.datatypes import Parameter
 
@@ -130,13 +130,17 @@ def forward(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         """
         Forward method of the 3 parameter Helmert.
         """
-        return coordinates + self.T
+        coords = coordinates.copy()
+        coords[:, 0:3] = coords[:, 0:3] + self.T
+        return coords
 
     def inverse(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         """
         Inverse method of the 3 parameter Helmert.
         """
-        return coordinates - self.T
+        coords = coordinates.copy()
+        coords[:, 0:3] = coords[:, 0:3] - self.T
+        return coords
 
     def estimate(
         self,
@@ -152,8 +156,12 @@ def estimate(
         this method is called.
         """
 
-        avg_source = np.average(source_coordinates, axis=0, weights=source_weights)
-        avg_target = np.average(target_coordinates, axis=0, weights=target_weights)
+        avg_source = np.average(
+            source_coordinates[:, 0:3], axis=0, weights=source_weights
+        )
+        avg_target = np.average(
+            target_coordinates[:, 0:3], axis=0, weights=target_weights
+        )
 
         mean_translation = avg_target - avg_source
 
@@ -172,7 +180,7 @@ class RotationConvention(enum.Enum):
 
 
 # fmt: off
-def R1(rx: float) -> CoordinateMatrix:
+def R1(rx: float) -> Matrix:
     """Rotation matrix for rotating about the x-axis."""
     return np.array(
         [
@@ -182,7 +190,7 @@ def R1(rx: float) -> CoordinateMatrix:
         ]
     )
 
-def R2(ry: float) -> CoordinateMatrix:
+def R2(ry: float) -> Matrix:
     """Rotation matrix for rotating about the y-axis."""
     return np.array(
         [
@@ -192,7 +200,7 @@ def R2(ry: float) -> CoordinateMatrix:
         ]
     )
 
-def R3(rz: float) -> CoordinateMatrix:
+def R3(rz: float) -> Matrix:
     """Rotation matrix for rotating about the z-axis."""
     return np.array(
         [
@@ -296,7 +304,7 @@ def _parameter_list(self) -> list[Parameter]:
         return params
 
     @cached_property
-    def R(self) -> CoordinateMatrix:  # pylint: disable=invalid-name
+    def R(self) -> Matrix:  # pylint: disable=invalid-name
         """
         Rotation matrix.
         """
@@ -342,13 +350,18 @@ def forward(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         # from the single-point Helmert formulation of B = T + s * R*A.
         # By transposing the rotation matrix we get the same results
         # when instead doing B = T + s * A*R^T.
-        return self.T + self.scale * coordinates @ self.R.T
+        coords = coordinates.copy()
+        coords[:, 0:3] = self.T + self.scale * coords[:, 0:3] @ self.R.T
+        return coords
 
     def inverse(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         """
         Inverse method of the 7 parameter Helmert.
         """
-        return -self.T + 1 / self.scale * coordinates @ self.R
+        coords = coordinates.copy()
+        coords[:, 0:3] = -self.T + 1 / self.scale * coords[:, 0:3] @ self.R
+
+        return coords
 
     def estimate(
         self,
@@ -382,7 +395,7 @@ def rad2arcsec(rad) -> float:
 
             A[i * 3 : i * 3 + 3, :] = np.column_stack([np.eye(3), c[0:3], R])
 
-        b = target_coordinates.flatten()
+        b = target_coordinates[:, 0:3].flatten()
 
         coeffs, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
 

src/transformo/operators/proj.py

@@ -50,22 +50,24 @@ def forward(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         """
         Forward method of the Transformation.
         """
-        x, y, z = self._transformer.transform(
+        x, y, z, t = self._transformer.transform(
             coordinates[:, 0],
             coordinates[:, 1],
             coordinates[:, 2],
+            coordinates[:, 3],
             direction=pyproj.enums.TransformDirection.FORWARD,
         )
-        return np.array([x, y, z]).transpose()
+        return np.array([x, y, z, t]).transpose()
 
     def inverse(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         """
         Inverse method of the Transformation.
         """
-        x, y, z = self._transformer.transform(
+        x, y, z, t = self._transformer.transform(
             coordinates[:, 0],
             coordinates[:, 1],
             coordinates[:, 2],
+            coordinates[:, 3],
             direction=pyproj.enums.TransformDirection.INVERSE,
         )
-        return np.array([x, y, z]).transpose()
+        return np.array([x, y, z, t]).transpose()

src/transformo/pipeline.py

@@ -81,6 +81,29 @@ def __init__(
         self._combined_source_data = sum(self.source_data, DataSource(None))
         self._combined_target_data = sum(self.target_data, DataSource(None))
 
+        all_source_stations = self._combined_source_data.stations
+        all_target_stations = self._combined_target_data.stations
+
+        station_union = self._combined_source_data.station_union(
+            self._combined_target_data
+        )
+
+        self._combined_source_data.limit_to_stations(station_union)
+        self._combined_target_data.limit_to_stations(station_union)
+
+        limited_source_stations = self._combined_source_data.stations
+        limited_target_stations = self._combined_target_data.stations
+
+        if len(limited_source_stations) != len(all_source_stations):
+            n = len(all_source_stations) - len(limited_source_stations)
+            msg = f"{n} stations removed from source data (not in intersection of source and target data)"
+            logger.warning(msg)
+
+        if len(limited_target_stations) != len(all_target_stations):
+            n = len(all_target_stations) - len(limited_target_stations)
+            msg = f"{n} stations removed from target data (not in intersection of source and target data)"
+            logger.warning(msg)
+
         self._intermediate_results: list[DataSource] = []
 
     @classmethod

src/transformo/presenters/coordinates.py

@@ -212,7 +212,7 @@ def evaluate(
         target = target_data.coordinate_matrix
         model = results[-1].coordinate_matrix
 
-        residuals = np.subtract(target, model)
+        residuals = np.subtract(target[:, 0:3], model[:, 0:3])
 
         # convert to milimeters
         residuals *= 1000

src/transformo/transformer.py

@@ -7,15 +7,15 @@
 import numpy as np
 import pyproj
 
-from transformo._typing import CoordinateMatrix, Vector
+from transformo._typing import CoordinateMatrix, CoordinateVector
 
 
 class Transformer:
     """
     Transform coordinates using PROJ.
 
     Interface for pyproj that works using Transformo's datatypes such
-    as CoordinateMatrix and Vector.
+    as CoordinateMatrix and CoordinateVector.
     """
 
     def __init__(self, transformer: pyproj.Transformer | None = None):
@@ -50,7 +50,7 @@ def transform_many(self, coordinates: CoordinateMatrix) -> CoordinateMatrix:
         results = self.transformer.itransform(coordinates)
         return np.array(list(results))
 
-    def transform_one(self, coordinate: Vector) -> Vector:
+    def transform_one(self, coordinate: CoordinateVector) -> CoordinateVector:
         """
         Transform a single coordinate.
         """

test/datasources/test_datasource.py

@@ -23,23 +23,23 @@ def test_datasource(coordinate_factory: Callable) -> None:
 
     # Is the coordinate matrix property functioning as expected?
     assert isinstance(ds1.coordinate_matrix, np.ndarray)
-    assert ds1.coordinate_matrix.shape == (n, 3)
+    assert ds1.coordinate_matrix.shape == (n, 4)
     assert (ds1.coordinate_matrix[0, :] == ds1.coordinates[0].vector).all()
 
     # Can we add two datasources?
     ds2 = DataSource(coordinates=[coordinate_factory() for _ in range(n)])
 
     ds_combined = ds1 + ds2
     assert len(ds_combined.coordinates) == 2 * n
-    assert ds_combined.coordinate_matrix.shape == (2 * n, 3)
+    assert ds_combined.coordinate_matrix.shape == (2 * n, 4)
     assert isinstance(ds_combined, DataSource)
 
     ds_combined2 = ds1
     for ds in [ds1, ds2]:
         ds_combined2 = ds_combined2 + ds
 
     assert len(ds_combined2.coordinates) == 3 * n
-    assert ds_combined2.coordinate_matrix.shape == (3 * n, 3)
+    assert ds_combined2.coordinate_matrix.shape == (3 * n, 4)
     assert isinstance(ds_combined2, DataSource)
 
 
@@ -238,3 +238,45 @@ def test_station_overrides(files) -> None:
         if c.station == "SULD":
             assert c.sz == 0.42
             assert c.station == "MULD"
+
+
+def test_station_union(files) -> None:
+    """
+    Test DataSource.station_union().
+
+    In the file `dk_cors_etrs89.csv` we have the following stations:
+
+        BUDP, ESBC, FER5, FYHA, GESR, HABY, HIRS, SMID, SULD, TEJH
+
+    By renaming some of them in the import we can simulate to different
+    DataSources that includes different stations.
+    """
+
+    ds1 = CsvDataSource(filename=files["dk_cors_etrs89.csv"])
+    ds2 = CsvDataSource(
+        filename=files["dk_cors_etrs89.csv"],
+        overrides={
+            "BUDP": CoordinateOverrides(station="BIDP"),
+            "FYHA": CoordinateOverrides(station="PYHA"),
+            "HIRS": CoordinateOverrides(station="HALS"),
+            "SULD": CoordinateOverrides(station="MULD"),
+        },
+    )
+
+    union = ds1.station_union(ds2)
+    expected_union = ["ESBC", "FER5", "GESR", "HABY", "SMID", "TEJH"]
+
+    assert set(union) == set(expected_union)
+
+
+def test_limit_to_stations(files) -> None:
+    """Test DataSource.limit_to_stations()."""
+    ds = CsvDataSource(
+        filename=files["dk_cors_etrs89.csv"],
+    )
+
+    stations = ["BUDP", "ESBC"]
+    ds.limit_to_stations(stations)
+    limited_stations = [c.station for c in ds.coordinates]
+
+    assert stations == limited_stations

test/operators/test_helmert.py

@@ -102,7 +102,7 @@ def test_helmerttranslation_as_operator():
     op = HelmertTranslation(x=3, y=5, z=10)
 
     # Can we roundtrip the `forward` and `inverse` methods
-    source_coordinates = np.zeros(shape=(10, 3))
+    source_coordinates = np.zeros(shape=(10, 4))
     roundtripped_coordinates = op.inverse(op.forward(source_coordinates))
     print(roundtripped_coordinates)
     assert np.all(source_coordinates == roundtripped_coordinates)
@@ -218,8 +218,12 @@ def test_helmert7param_transformation(source_coordinates):
     print(projstring)
     T = Transformer.from_projstring(projstring)
 
+    print(source_coordinates[0, :])
     transformo_coords = h7.forward(source_coordinates)
     proj_coords = T.transform_many(source_coordinates)
+    print(source_coordinates[0, :])
+    print(transformo_coords[0, :])
+    print(proj_coords[0, :])
 
     # We are not going to get an exact match but it's very close and good enough
     assert np.allclose(proj_coords, transformo_coords)