Add river network topology determination

rivertopo/topology.py

@@ -0,0 +1,41 @@
+from osgeo import ogr
+import numpy as np
+
+ogr.UseExceptions()
+
+class ConnectedPoints:
+    def __init__(self, upstream=frozenset(), downstream=frozenset()):
+        self.upstream = set(upstream)
+        self.downstream = set(downstream)
+
+    def __repr__(self):
+        return f'ConnectedPoints(upstream={self.upstream}, downstream={self.downstream})'
+
+    def __eq__(self, other):
+        return self.upstream == other.upstream and self.downstream == other.downstream
+
+def get_layer_topology(layer):
+    # dictionary of ConnectedPoint, indexed by coordinates
+    point_connections = {}
+
+    for feature in layer:
+        geometry = feature.GetGeometryRef()
+        # geometry = feature.GetGeometryRef().GetGeometryRef(0)
+
+        # TODO flip array if faldretning == "Modsat"?
+        geometry_xy = np.array(geometry.GetPoints())[:,0:2]
+        geometry_xy_tuples = [tuple(xy) for xy in geometry_xy]
+
+        # initialize
+        for xy in geometry_xy_tuples:
+            if not xy in point_connections:
+                point_connections[xy] = ConnectedPoints()
+
+        # Store upstream and downstream points
+        for i in range(1, len(geometry_xy_tuples)):
+            point_connections[geometry_xy_tuples[i]].upstream.add(geometry_xy_tuples[i-1])
+
+        for i in range(0, len(geometry_xy_tuples)-1):
+            point_connections[geometry_xy_tuples[i]].downstream.add(geometry_xy_tuples[i+1])
+
+    return point_connections

tests/conftest.py

@@ -15,21 +15,30 @@ def driver():
     return ogr.GetDriverByName('Memory')
 
 @pytest.fixture()
-def datasrc(driver):
+def centerline_datasrc(driver):
+    return driver.CreateDataSource('river_centerlines')
+
+@pytest.fixture()
+def profile_datasrc(driver):
     return driver.CreateDataSource('vandloebsdata')
 
 @pytest.fixture()
-def regulativprofilsimpel_layer(datasrc, srs):
-    layer = datasrc.CreateLayer('regulativprofilsimpel_nohist', srs=srs, geom_type=ogr.wkbPoint)
+def centerline_layer(centerline_datasrc, srs):
+    layer = centerline_datasrc.CreateLayer('vandloebsmidte', srs=srs, geom_type=ogr.wkbLineString25D)
+    return layer
+
+@pytest.fixture()
+def regulativprofilsimpel_layer(profile_datasrc, srs):
+    layer = profile_datasrc.CreateLayer('regulativprofilsimpel_nohist', srs=srs, geom_type=ogr.wkbPoint)
     layer.CreateField(ogr.FieldDefn('anlaeghoejre', ogr.OFTReal))
     layer.CreateField(ogr.FieldDefn('anlaegvenstre', ogr.OFTReal))
     layer.CreateField(ogr.FieldDefn('bundbredde', ogr.OFTReal))
     layer.CreateField(ogr.FieldDefn('bundkote', ogr.OFTReal))
     return layer
 
 @pytest.fixture()
-def regulativprofilsammensat_layer(datasrc, srs):
-    layer = datasrc.CreateLayer('regulativprofilsammens_nohist', srs=srs, geom_type=ogr.wkbPoint)
+def regulativprofilsammensat_layer(profile_datasrc, srs):
+    layer = profile_datasrc.CreateLayer('regulativprofilsammens_nohist', srs=srs, geom_type=ogr.wkbPoint)
     layer.CreateField(ogr.FieldDefn('afsatbanketbreddehoejre', ogr.OFTReal))
     layer.CreateField(ogr.FieldDefn('afsatbanketbreddevenstre', ogr.OFTReal))
     layer.CreateField(ogr.FieldDefn('afsatsanlaeghoejre', ogr.OFTReal))
@@ -42,6 +51,6 @@ def regulativprofilsammensat_layer(datasrc, srs):
     return layer
 
 @pytest.fixture()
-def opmaaltprofil_layer(datasrc, srs):
-    layer = datasrc.CreateLayer('opmaaltprofil_nohist', srs=srs, geom_type=ogr.wkbLineString25D)
+def opmaaltprofil_layer(profile_datasrc, srs):
+    layer = profile_datasrc.CreateLayer('opmaaltprofil_nohist', srs=srs, geom_type=ogr.wkbLineString25D)
     return layer

tests/test_topology.py

@@ -0,0 +1,54 @@
+from rivertopo import topology
+from rivertopo.topology import ConnectedPoints, get_layer_topology
+
+from osgeo import ogr
+import numpy as np
+
+ogr.UseExceptions()
+
+def test_get_layer_topology(centerline_layer):
+    layer_defn = centerline_layer.GetLayerDefn()
+
+    river_points = [
+        [
+            (0.0, 0.0, 10.0),
+            (4.0, 3.0, 9.5),
+            (5.0, 3.0, 9.0),
+            (7.0, 2.0, 8.5),
+            (8.0, 0.0, 8.0),
+        ],
+        [
+            (1.0, 4.0, 10.0),
+            (4.0, 3.0, 9.5),
+        ],
+        [
+            (7.0, 2.0, 8.5),
+            (9.0, 3.0, 8.0),
+            (10.0, 3.0, 7.5),
+        ],
+    ]
+
+    expected_connections = {
+        river_points[0][0][:2]: ConnectedPoints(downstream=set([river_points[0][1][:2]])),
+        river_points[0][1][:2]: ConnectedPoints(upstream=set([river_points[0][0][:2], river_points[1][0][:2]]), downstream=set([river_points[0][2][:2]])),
+        river_points[0][2][:2]: ConnectedPoints(upstream=set([river_points[0][1][:2]]), downstream=set([river_points[0][3][:2]])),
+        river_points[0][3][:2]: ConnectedPoints(upstream=set([river_points[0][2][:2]]), downstream=set([river_points[0][4][:2], river_points[2][1][:2]])),
+        river_points[0][4][:2]: ConnectedPoints(upstream=set([river_points[0][3][:2]])),
+        river_points[1][0][:2]: ConnectedPoints(downstream=set([river_points[0][1][:2]])),
+        river_points[2][1][:2]: ConnectedPoints(upstream=set([river_points[0][3][:2]]), downstream=set([river_points[2][2][:2]])),
+        river_points[2][2][:2]: ConnectedPoints(upstream=set([river_points[2][1][:2]])),
+    }
+
+    for coord_array in river_points:
+        geometry = ogr.Geometry(ogr.wkbLineString25D)
+        for coords in coord_array:
+            geometry.AddPoint(coords[0], coords[1], coords[2])
+
+        feature = ogr.Feature(layer_defn)
+        feature.SetGeometry(geometry)
+
+        centerline_layer.CreateFeature(feature)
+
+    topology = get_layer_topology(centerline_layer)
+
+    assert topology == expected_connections