Fix geo conversion

Author: arneboockmeyer

test/node_test.py

@@ -2,7 +2,8 @@
 
 from pytest import raises
 
-from yaramo.model import Node
+import yaramo.utils.coordinateconversion
+from yaramo.model import DbrefGeoNode, Node, Wgs84GeoNode
 
 from .helper import create_edge, create_node
 
@@ -150,3 +151,34 @@ def test_implausible_anschluss():
 
     with raises(Exception) as exception:
         point.calc_anschluss_of_all_nodes()
+
+
+def test_wgs84_dbref_coordinate_conversion():
+    def _test_distance(_dbref_x, _dbref_y, _wgs84_x, _wgs84_y, factor, smaller_as):
+        wgs84_geo_node = Wgs84GeoNode(wgs84_x, wgs84_y)
+        dbref_geo_node = DbrefGeoNode(dbref_x, dbref_y)
+        assert wgs84_geo_node.get_distance_to_other_geo_node(dbref_geo_node) * factor < smaller_as
+
+    # Example 1:
+    dbref_x = 4563230.251887853
+    dbref_y = 5601992.441701063
+    wgs84_x = 50.55025861737121
+    wgs84_y = 12.890666340087865
+
+    _test_distance(dbref_x, dbref_y, wgs84_x, wgs84_y, 1000 * 1000, 2)
+
+    # Example 2:
+    dbref_x = 4563437.90802
+    dbref_y = 5601946.51808
+    wgs84_x = 50.54982337298292
+    wgs84_y = 12.893588101538324
+
+    _test_distance(dbref_x, dbref_y, wgs84_x, wgs84_y, 1000 * 1000, 2)
+
+    # Example 3: (close-by)
+    dbref_x = 4564720.99586
+    dbref_y = 5601735.46336
+    wgs84_x = 50.5477883
+    wgs84_y = 12.9116547
+
+    _test_distance(dbref_x, dbref_y, wgs84_x, wgs84_y, 1, 0.3)

yaramo/geo_node.py

@@ -1,21 +1,24 @@
 import math
 from abc import ABC, abstractmethod
 
-import pyproj
+from haversine import Unit, haversine
 
 from yaramo.base_element import BaseElement
 
+from .utils.coordinateconversion import transform_dbref_to_wgs84, transform_wgs84_to_dbref
+
 
 class GeoNode(ABC, BaseElement):
     """This is the baseclass of specific GeoNodes that use different coordinate systems.
 
     A GeoNode is characterized by it's x and y coordinates.
     """
 
-    def __init__(self, x, y, **kwargs):
+    def __init__(self, x, y, dbref_crs: str = "ER0", **kwargs):
         super().__init__(**kwargs)
         self.x = x
         self.y = y
+        self.dbref_crs = dbref_crs
 
     @abstractmethod
     def get_distance_to_other_geo_node(self, geo_node_b: "GeoNode"):
@@ -39,66 +42,43 @@ def to_serializable(self):
 
 
 class Wgs84GeoNode(GeoNode):
-    def get_distance_to_other_geo_node(self, geo_node_b: "Wgs84GeoNode"):
-        assert type(self) == type(
-            geo_node_b
-        ), "You cannot calculate the distance between a Wgs84GeoNode and a geo node of a different type!"
-        return self.__haversine_distance(geo_node_b) / 1000
+    def get_distance_to_other_geo_node(self, geo_node_b: "GeoNode"):
+        geo_node_b = geo_node_b.to_wgs84()
+        return self.__haversine_distance(geo_node_b)
 
     def __haversine_distance(self, geo_node_b: "GeoNode"):
-        pi_over_180 = float(math.pi / 180)
-        return (
-            2
-            * 6371000
-            * math.asin(
-                math.pi
-                / 180
-                * math.sqrt(
-                    math.pow(math.sin((pi_over_180 * (geo_node_b.x - self.x)) / 2), 2)
-                    + math.cos(pi_over_180 * self.x)
-                    * math.cos(pi_over_180 * geo_node_b.x)
-                    * math.pow(math.sin((pi_over_180 * (geo_node_b.y - self.y)) / 2), 2)
-                )
-            )
-        )
-
-    def to_wgs84(self):
+        own = (self.x, self.y)
+        other = (geo_node_b.x, geo_node_b.y)
+        return haversine(own, other, unit=Unit.METERS)
+
+    def to_wgs84(self) -> "Wgs84GeoNode":
         return self
 
-    def to_dbref(self):
-        transformer = pyproj.Transformer.from_crs("epsg:4326", "epsg:31468")
-        x, y = transformer.transform(self.y, self.x)
-        return DbrefGeoNode(x, y)
+    def to_dbref(self) -> "DbrefGeoNode":
+        x, y = transform_wgs84_to_dbref(self.x, self.y, self.dbref_crs)
+        return DbrefGeoNode(x, y, self.dbref_crs)
 
     def to_euclidean(self) -> "EuclideanGeoNode":
-        raise NotImplementedError
+        return self.to_dbref().to_euclidean()
 
 
 class DbrefGeoNode(GeoNode):
-    def get_distance_to_other_geo_node(self, geo_node_b: "DbrefGeoNode"):
-        assert type(self) == type(
-            geo_node_b
-        ), "You cannot calculate the distance between a DbrefGeoNode and a geo node of a different type!"
-        return self.__eucldian_distance(geo_node_b)
-
-    def __eucldian_distance(self, geo_node_b: "GeoNode"):
-        min_x = min(self.x, geo_node_b.x)
-        min_y = min(self.y, geo_node_b.y)
-        max_x = max(self.x, geo_node_b.x)
-        max_y = max(self.y, geo_node_b.y)
-        return math.sqrt(math.pow(max_x - min_x, 2) + math.pow(max_y - min_y, 2))
+    def get_distance_to_other_geo_node(self, geo_node_b: "GeoNode"):
+        # Separate DB Ref distance method not implemented yet, therefore use WGS84 distance
+        return self.to_wgs84().get_distance_to_other_geo_node(geo_node_b)
 
-    def to_wgs84(self):
-        raise NotImplementedError
+    def to_wgs84(self) -> "Wgs84GeoNode":
+        x, y = transform_dbref_to_wgs84(self.x, self.y, self.dbref_crs)
+        return Wgs84GeoNode(x, y, self.dbref_crs)
 
-    def to_dbref(self):
+    def to_dbref(self) -> "DbrefGeoNode":
         return self
 
     def to_euclidean(self) -> "EuclideanGeoNode":
         # This transformation is just for testing purposes and not correct, see documentation in EuclideanGeoNode.
         _x_shift = 4533770.0
         _y_shift = 5625780.0
-        return EuclideanGeoNode(self.x - _x_shift, self.y - _y_shift)
+        return EuclideanGeoNode(self.x - _x_shift, self.y - _y_shift, self.dbref_crs)
 
 
 class EuclideanGeoNode(GeoNode):
@@ -109,9 +89,7 @@ class EuclideanGeoNode(GeoNode):
     """
 
     def get_distance_to_other_geo_node(self, geo_node_b: "EuclideanGeoNode"):
-        assert type(self) == type(
-            geo_node_b
-        ), "You cannot calculate the distance between a EuclideanGeoNode and a geo node of a different type!"
+        geo_node_b = geo_node_b.to_euclidean()
         return self.__eucldian_distance(geo_node_b)
 
     def __eucldian_distance(self, geo_node_b: "GeoNode"):
@@ -122,12 +100,12 @@ def __eucldian_distance(self, geo_node_b: "GeoNode"):
         return math.sqrt(math.pow(max_x - min_x, 2) + math.pow(max_y - min_y, 2))
 
     def to_wgs84(self) -> "Wgs84GeoNode":
-        raise NotImplementedError
+        return self.to_dbref().to_wgs84()
 
-    def to_dbref(self):
+    def to_dbref(self) -> "DbrefGeoNode":
         _x_shift = 4533770.0
         _y_shift = 5625780.0
-        return DbrefGeoNode(self.x + _x_shift, self.y + _y_shift)
+        return DbrefGeoNode(self.x + _x_shift, self.y + _y_shift, self.dbref_crs)
 
     def to_euclidean(self) -> "EuclideanGeoNode":
         return self

yaramo/utils/__init__.py

@@ -0,0 +1,65 @@
+"""
+Copyright (c) 2021 DB Netz AG and others.
+
+All rights reserved. This program and the accompanying materials
+are made available under the terms of the Eclipse Public License v2.0
+which accompanies this distribution, and is available at
+http://www.eclipse.org/legal/epl-v20.html
+
+
+-------------
+
+The DBRef and WGS84 conversions are taken from the Eclipse Signalling Engineering Toolbox (https://github.com/eclipse-set/set/tree/main)
+
+See
+https://github.com/eclipse-set/set/blob/main/java/bundles/org.eclipse.set.feature.siteplan/src/org/eclipse/set/feature/siteplan/positionservice/PositionServiceImpl.java
+and
+https://github.com/eclipse-set/set/blob/main/java/bundles/org.eclipse.set.ppmodel.extensions/src/org/eclipse/set/ppmodel/extensions/geometry/CoordinateExtensions.xtend
+for the original sources.
+"""
+import pyproj
+
+EPSG_3857_STRING = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
+
+
+def transform_dbref_to_wgs84(x: float, y: float, coordinate_str: str):
+    proj_dbref_str = __get_proj_dbref_str(coordinate_str)
+
+    # Step 1 DBRef to Pseudo-Mercator:
+    transformer = pyproj.Transformer.from_crs(crs_from=proj_dbref_str, crs_to=EPSG_3857_STRING)
+    x, y = transformer.transform(x, y)
+
+    # Step 2 Pseudo-Mercator to WGS84:
+    transformer_2 = pyproj.Transformer.from_crs(crs_from=3857, crs_to=4326)
+    return transformer_2.transform(x, y)
+
+
+def transform_wgs84_to_dbref(x: float, y: float, coordinate_str: str):
+    proj_dbref_str = __get_proj_dbref_str(coordinate_str)
+
+    # Step 2 WGS84 to Pseudo-Mercator:
+    transformer_2 = pyproj.Transformer.from_crs(crs_from=4326, crs_to=3857)
+    x, y = transformer_2.transform(x, y)
+
+    # Step 1 Pseudo-Mercator to DBRef:
+    transformer = pyproj.Transformer.from_crs(crs_from=EPSG_3857_STRING, crs_to=proj_dbref_str)
+    return transformer.transform(x, y)
+
+
+def __get_proj_dbref_str(coordinate_str: str):
+    CRS_CR0_STRING = "+proj=tmerc +lat_0=0 +lon_0=6  +k=1 +x_0=2500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs"
+    CRS_DR0_STRING = "+proj=tmerc +lat_0=0 +lon_0=9  +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs"
+    CRS_ERO_STRING = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs"
+    CRS_FRO_STRING = "+proj=tmerc +lat_0=0 +lon_0=15 +k=1 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs"
+
+    match coordinate_str:
+        case "CR0":
+            return CRS_CR0_STRING
+        case "DR0":
+            return CRS_DR0_STRING
+        case "ER0":
+            return CRS_ERO_STRING
+        case "FR0":
+            return CRS_FRO_STRING
+        case _:
+            raise ValueError("No valid coordinate_str given. Supported are: CR0, DR0, ER0 and FR0.")