Comparing two sets of polylines intersecting polygons and get a matrix of gained/lost connections

[GuyRosenblum]

I have three feature classes, one is of Polygons and the other two of Polylines, representing a set of areas in a city and a network of transit lines between those areas - one network for current service and another with some changes in it. There's no topological connection between the feature classes.

1. First I'd like to map out the service between all polygons in the city, in pairs - like in an Adjacency Matrix. As in "How many lines are connecting A & B?" Doing so for each set of two polygons. This process will repeat itself for both sets of polylines, in order to build two adjacency matrices, i.e. how many routes connect each pair of polygons.
2. Then I'd like to compare the results of both networks and find how the service between areas of the city has changed - and return a numeric value of lost/gained connections. For example in TABLE 1, there's one new connection between A&C, and zero new connections between C&D (which has changed between lines yet retained one connecting line).
3. If that's too complicated it would be sufficient to get something like in TABLE 2, in which every pair of polygons gets one field, and each line gets a row. In that way I could see which changes occurred and highlight problems at the individual line level.

I thought of doing many INTERSECT commands for each polygon and to calculate pairs manually in Excel, but how can I automate the process for example for 80 lines * 30 Polygons or more? I use both ArcMap and QGIS but I'm far from an expert in those...

How can I solve this, and is this considered complicated?

[JohannesLindner · Accepted solution]

def analyze_connections(area_fc, line_fc, area_id_field, line_id_field):
    """Analyzes how areas are connected by lines.
    
    area_fc: str, layer name or path of the areas (polygons)
    line_fc: str, layer name or path of the lines (polylines)
    area_id_field: str, uniqe identifier of the areas
    line_id_field: str, unique identifier of the lines
    
    Returns a dict of dicts, specifiying the areas and the count of direct connections to other areas

    """
    areas = [row for row in arcpy.da.SearchCursor(area_fc, ["SHAPE@", area_id_field])]
    lines = [row for row in arcpy.da.SearchCursor(line_fc, ["SHAPE@", line_id_field])]
    areas_on_line = {line[1]: [area[1] for area in areas if not line[0].disjoint(area[0])] for line in lines}
    areas = [area[1] for area in areas]
    area_connections = dict()
    for area in areas:
        area_connections[area] = dict()
        for other_area in areas:
            connections = 0
            for l, a in areas_on_line.items():
                if area != other_area:
                    if area in a and other_area in a:
                        connections += 1
            area_connections[area][other_area] = connections
    return area_connections

def compare_connections(old, new):
    """Compares two connection dictionaries (eg output of analyze_connections).

    old: dict of dicts, the old connections
    new: dict of dicts, the new connections
    
    Returns a dict of dicts, specifying the difference between the connection counts of old and new connections.

    """
    areas = sorted(list(old.keys()))
    diff = dict()
    for area in areas:
        diff[area] = dict()
        for other_area  in areas:
            diff[area][other_area] = new[area][other_area] - old[area][other_area]
    return diff

def print_adjacency_matrix(connections):
    """Takes a dict of dicts (eg output of analyze_connections) and prints
    a 2D matrix.

    """"
    areas = sorted(list(connections.keys()))
    rows = [[connections[area][other_area] for other_area in areas] for area in areas]
    for i, row in enumerate(rows):
        rows[i] = "\t".join([areas[i]] + [str(r) for r in row])
    rows.insert(0, "\t".join([""] + areas))
    print("\n".join(rows))

old_connections = analyze_connections("TestPolygons", "TestLines", "TextField", "TextField")
print_adjacency_matrix(old_connections)
#	A	B	C	D
#A	0	1	0	1
#B	1	0	2	2
#C	0	2	0	1
#D	1	2	1	0

new_connections = analyze_connections("TestPolygons", "TestLines", "TextField", "TextField")
print_adjacency_matrix(new_connections)
#	A	B	C	D
#A	0	0	1	0
#B	0	0	1	0
#C	1	1	0	1
#D	0	0	1	0

difference = compare_connections(old_connections, new_connections)
print_adjacency_matrix(difference)
#	A	B	C	D
#A	0	-1	1	-1
#B	-1	0	-1	-2
#C	1	-1	0	0
#D	-1	-2	0	0

[GuyRosenblum]

Thank you so much Johannes!

I will try to check with my peers that are experienced in Python so we could solve this once and for all.