wip: interpret network from vector features

sample_data/use_cases/new_york_energy/nysdp.py

@@ -1,9 +1,12 @@
 import json
 import geopandas
+import numpy
+import pandas
 import requests
+import shapely
 
-from django.contrib.gis.geos import GEOSGeometry
-from uvdat.core.models import Dataset, VectorMapLayer, VectorFeature
+from django.contrib.gis.geos import GEOSGeometry, Point, LineString
+from uvdat.core.models import Dataset, Network, NetworkNode, NetworkEdge, VectorMapLayer, VectorFeature
 
 
 NYDSP_URL = 'https://systemdataportal.nationalgrid.com/arcgis/rest/services/NYSDP'
@@ -55,53 +58,219 @@ def create_vector_features(dataset, service_name=None):
     VectorFeature.objects.bulk_create(vector_features)
 
 
+def vector_features_from_network(network):
+    VectorMapLayer.objects.filter(dataset=network.dataset).delete()
+    map_layer = VectorMapLayer.objects.create(dataset=network.dataset)
+    VectorFeature.objects.bulk_create([
+        VectorFeature(
+            map_layer=map_layer,
+            geometry=node.location,
+            properties=node.metadata or {},
+        )
+        for node in network.nodes.all()
+    ])
+    VectorFeature.objects.bulk_create([
+        VectorFeature(
+            map_layer=map_layer,
+            geometry=edge.line_geometry,
+            properties=edge.metadata or {},
+        )
+        for edge in network.edges.all()
+    ])
+
+
 def create_consolidated_network(dataset):
-    VectorMapLayer.objects.filter(dataset=dataset).delete()
-    map_layer = VectorMapLayer.objects.create(dataset=dataset, index=0)
-    include_services = [
-        "Substations",
-        "DistAssetsOverview",
-        "Electrification_Data",
-        "EV_Load_Serving_Capacity",
-        "Hosting_Capacity_Data",
-        "LSRV",
-        "NY_SubT_SDP"
-    ]
+    # VectorMapLayer.objects.filter(dataset=dataset).delete()
+    # map_layer = VectorMapLayer.objects.create(dataset=dataset, index=0)
+    # include_services = [
+    #     "Substations",
+    #     "DistAssetsOverview",
+    #     "Electrification_Data",
+    #     "EV_Load_Serving_Capacity",
+    #     "Hosting_Capacity_Data",
+    #     "LSRV",
+    #     "NY_SubT_SDP"
+    # ]
 
-    features = []
-    for service_name in include_services:
-        try:
-            dataset = Dataset.objects.get(name=f'National Grid {service_name}')
-            print('\t', f'Using saved vector features for {service_name}...')
-            features += [
-                dict(
-                    type='Feature',
-                    geometry=json.loads(feature.geometry.json),
-                    properties=feature.properties
-                )
-                for feature in VectorFeature.objects.filter(map_layer__dataset=dataset)
-            ]
-        except Dataset.DoesNotExist:
-            print('\t', f'Querying {service_name}...')
-            feature_sets = fetch_vector_features(service_name=service_name)
-            for layer_id, feature_set in feature_sets.items():
-                features += feature_set
-    print('\t', len(features), 'features found. Combining like geometries...')
+    # features = []
+    # for service_name in include_services:
+    #     try:
+    #         d = Dataset.objects.get(name=f'National Grid {service_name}')
+    #         print('\t', f'Using saved vector features for {service_name}...')
+    #         features += [
+    #             dict(
+    #                 type='Feature',
+    #                 geometry=json.loads(feature.geometry.json),
+    #                 properties=feature.properties
+    #             )
+    #             for feature in VectorFeature.objects.filter(map_layer__dataset=d)
+    #         ]
+    #     except Dataset.DoesNotExist:
+    #         print('\t', f'Querying {service_name}...')
+    #         feature_sets = fetch_vector_features(service_name=service_name)
+    #         for layer_id, feature_set in feature_sets.items():
+    #             features += feature_set
+    # print('\t', len(features), 'features found. Combining like geometries...')
+    # print('features found:', len(features))
 
     # Use geopandas to merge properties of duplicate features
+    # gdf = geopandas.GeoDataFrame.from_features(features)
+    # gdf["geometry"] = gdf.normalize()
+
+    # grouped = gdf.groupby(gdf['geometry'].to_wkt()).first()
+    # print('\t', len(grouped), 'consolidated features found. Saving to database...')
+
+    # TODO: color by properties
+
+    # geojson = json.loads(grouped.to_json())
+    # VectorFeature.objects.bulk_create([
+    #     VectorFeature(
+    #         map_layer=map_layer,
+    #         geometry=GEOSGeometry(json.dumps(feature['geometry'])),
+    #         properties=feature['properties'],
+    #     )
+    #     for feature in geojson.get('features')
+    # ], batch_size=10000)
+
+    from django.contrib.gis.geos import Polygon
+
+
+    test_dataset, created = Dataset.objects.get_or_create(
+        name="Test Network",
+        category="energy",
+        dataset_type="VECTOR"
+    )
+    Network.objects.filter(dataset=test_dataset).delete()
+    network = Network.objects.create(dataset=test_dataset)
+
+    test_bbox = dict(xmin=-74, ymin=42.8, xmax=-73.9, ymax=42.9)
+    test_area = Polygon.from_bbox(
+        [test_bbox['xmin'], test_bbox['ymin'], test_bbox['xmax'], test_bbox['ymax']]
+    )
+    features = [
+        dict(
+            type='Feature',
+            geometry=json.loads(feature.geometry.json),
+            properties=feature.properties
+        )
+        for feature in VectorFeature.objects.filter(
+            map_layer__dataset=dataset,
+            geometry__coveredby=test_area,
+        )
+    ]
+    show_columns = [
+        'geometry',
+    ]
+    pandas.set_option('display.max_columns', None)
+    print('test features:', len(features))
     gdf = geopandas.GeoDataFrame.from_features(features)
     gdf["geometry"] = gdf.normalize()
+    spatial_index = gdf.sindex
 
-    grouped = gdf.groupby(gdf['geometry'].to_wkt()).first()
-    intersects = grouped.sjoin(grouped, how='left', predicate='intersects')
-    print('\t', len(grouped), 'consolidated features found. Saving to database...')
+    # construct adjacencies list
+    adjacencies = {}
+    for feat_index, feature in gdf.iterrows():
+        adjacencies[feat_index] = {}
+        geom = feature['geometry']
+        if geom.geom_type == 'MultiLineString':
+            # attempt to merge contiguous portions
+            geom = shapely.ops.linemerge(geom)
+        geom_list = [geom]
+        if geom.geom_type == 'MultiLineString':
+            # if still multi, contains non-contiguous portions
+            geom_list = geom.geoms
 
-    geojson = json.loads(grouped.to_json())
-    VectorFeature.objects.bulk_create([
-        VectorFeature(
-            map_layer=map_layer,
-            geometry=GEOSGeometry(json.dumps(feature['geometry'])),
-            properties=feature['properties'],
-        )
-        for feature in geojson.get('features')
-    ])
+        for geom_index, geom in enumerate(geom_list):
+            adjacencies[feat_index][geom_index] = {}
+            # first and last points (endpoints of current geom)
+            # assumes no intersections will occur in the middle of a line
+            for coord_index, coord_name in [(0, 'start'), (-1, 'end')]:
+                coord = shapely.geometry.Point(geom.coords[coord_index])
+                nearest_indexes, distances = spatial_index.nearest(coord, max_distance=10, return_distance=True)
+                geom_inds, tree_inds = nearest_indexes
+                # exclude self from nearest results
+                tree_inds = [int(i) for i in tree_inds if i != feat_index]
+                node = None
+                if len(tree_inds) != 1:
+                    node = NetworkNode.objects.create(
+                        name=f'{feat_index}_{geom_index}_{coord_name}',
+                        network=network,
+                        location=Point(*geom.coords[coord_index])
+                    )
+                adjacencies[feat_index][geom_index][coord_name] = dict(
+                    node=node,
+                    adjs=tree_inds,
+                )
+
+    visited = []
+    def follow_adjacencies(feat_index, current_line=None, current_from_node=None):
+        if feat_index not in visited:
+            visited.append(feat_index)
+            feat = gdf.iloc[feat_index]
+            for geom_index, geom_record in adjacencies[feat_index].items():
+                start_record = geom_record.get('start', {})
+                start_node = start_record.get('node')
+                start_adjs = start_record.get('adjs')
+
+                end_record = geom_record.get('end', {})
+                end_node = end_record.get('node')
+                end_adjs = end_record.get('adjs')
+                coords = []
+                if feat['geometry'].geom_type == 'MultiLineString':
+                    for geom in feat['geometry'].geoms:
+                        coords += geom.coords
+                else:
+                    coords = feat['geometry'].coords
+
+                if start_node and end_node:
+                    # print('edge', start_node, end_node, 'len', len(coords))
+                    NetworkEdge.objects.create(
+                        name='edge',
+                        network=network,
+                        from_node=start_node,
+                        to_node=end_node,
+                        metadata={},  # TODO: figure out capturing metadata
+                        line_geometry=LineString(*[
+                            Point(*p) for p in coords
+                        ])
+                    )
+
+                if current_line is None:
+                    current_line = list(coords)
+                else:
+                    current_line += list(coords)
+
+                follow_start_results = None
+                follow_end_results = None
+                if len(start_adjs) > 0:
+                    for adj in start_adjs:
+                        follow_start_results = follow_adjacencies(
+                            adj, current_line=current_line, current_from_node=current_from_node or end_node
+                        )
+                if len(end_adjs) > 0:
+                    for adj in end_adjs:
+                        follow_end_results = follow_adjacencies(
+                            adj, current_line=current_line, current_from_node=current_from_node or start_node
+                        )
+
+                start_line, start_from_node, end_line, end_from_node = None, None, None, None
+                if follow_start_results:
+                    start_line, start_from_node = follow_start_results
+                if follow_end_results:
+                    end_line, end_from_node = follow_end_results
+                if current_from_node is None:
+                    if start_line is not None and end_line is not None:
+                        print(feat_index, 'combine lines', len(start_line), '+', len(end_lines), 'between', start_from_node, end_from_node)
+                    if start_line is not None:
+                        print(feat_index, 'start line!', start_from_node, end_from_node, current_from_node, start_node, end_node)
+                    if end_line is not None:
+                        print(feat_index, 'end line!', start_from_node, end_from_node, current_from_node, start_node, end_node)
+
+            return current_line, current_from_node
+
+    for feat_index in list(adjacencies.keys())[:5]:
+        follow_adjacencies(feat_index)
+    print('total nodes', NetworkNode.objects.filter(network=network).count())
+    print('total edges', NetworkEdge.objects.filter(network=network).count())
+
+    vector_features_from_network(network)