Merge branch 'yaramo2' into operations

Author: arneboockmeyer

test/helper.py

@@ -1,8 +1,8 @@
-from yaramo.model import Edge, Node, Wgs84GeoNode
+from yaramo.model import Edge, EuclideanGeoNode, Node
 
 
 def create_geo_node(x, y):
-    return Wgs84GeoNode(x, y)
+    return EuclideanGeoNode(x, y)
 
 
 def create_node(x, y):

test/node_test.py

@@ -2,8 +2,8 @@
 
 from pytest import raises
 
-from yaramo.geo_node import Wgs84GeoNode
-from yaramo.model import Node
+import yaramo.utils.coordinateconversion
+from yaramo.model import DbrefGeoNode, Node, Wgs84GeoNode
 
 from .helper import create_edge, create_node
 
@@ -151,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/edge.py

@@ -58,10 +58,15 @@ def get_other_node(self, node: Node) -> Node:
             return self.node_b
         return self.node_a
 
-    def update_length(self):
+    def update_length(self, force: bool = False):
+        if force:
+            self.length = None
         self.length = self.__get_length()
 
     def __get_length(self) -> float:
+        if self.length is not None:
+            return self.length
+
         if len(self.intermediate_geo_nodes) == 0:
             return self.node_a.geo_node.get_distance_to_other_geo_node(self.node_b.geo_node)
 

yaramo/geo_node.py

@@ -1,22 +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"):
@@ -31,36 +33,63 @@ def to_wgs84(self) -> "Wgs84GeoNode":
     def to_dbref(self) -> "DbrefGeoNode":
         pass
 
+    @abstractmethod
+    def to_euclidean(self) -> "EuclideanGeoNode":
+        pass
+
     def to_serializable(self):
         return self.__dict__, {}
 
 
 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 DbrefGeoNode!"
+    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"):
         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):
+    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":
+        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 Wgs84GeoNode!"
+    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) -> "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) -> "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, self.dbref_crs)
+
+
+class EuclideanGeoNode(GeoNode):
+    """
+    The Euclidean geo node is a simple geo node based on the simple concept of x and y coordinates and
+    the Euclidean distance between two geo nodes. It's mainly for testing purposes and can only be converted to
+    DBRef geo nodes, but the conversion is just a simple coordinate shift and with this, probably incorrect.
+    """
+
+    def get_distance_to_other_geo_node(self, geo_node_b: "EuclideanGeoNode"):
+        geo_node_b = geo_node_b.to_euclidean()
         return self.__eucldian_distance(geo_node_b)
 
     def __eucldian_distance(self, geo_node_b: "GeoNode"):
@@ -70,8 +99,13 @@ def __eucldian_distance(self, geo_node_b: "GeoNode"):
         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 to_wgs84(self):
-        raise NotImplementedError
+    def to_wgs84(self) -> "Wgs84GeoNode":
+        return self.to_dbref().to_wgs84()
+
+    def to_dbref(self) -> "DbrefGeoNode":
+        _x_shift = 4533770.0
+        _y_shift = 5625780.0
+        return DbrefGeoNode(self.x + _x_shift, self.y + _y_shift, self.dbref_crs)
 
-    def to_dbref(self):
+    def to_euclidean(self) -> "EuclideanGeoNode":
         return self

yaramo/model.py

@@ -1,5 +1,5 @@
 from yaramo.edge import Edge  # pylint: noqa
-from yaramo.geo_node import DbrefGeoNode, GeoNode, Wgs84GeoNode  # pylint: noqa
+from yaramo.geo_node import DbrefGeoNode, EuclideanGeoNode, GeoNode, Wgs84GeoNode  # pylint: noqa
 from yaramo.node import EdgeConnectionDirection, Node  # pylint: noqa
 from yaramo.route import Route  # pylint: noqa
 from yaramo.signal import (  # pylint: noqa

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.")