%config InlineBackend.figure_format = "retina"

%load_ext watermark
%watermark

import geopandas
import libpysal
from libpysal import examples
from libpysal.cg import Point, Chain
import matplotlib
import matplotlib_scalebar
from matplotlib_scalebar.scalebar import ScaleBar
import spaghetti

%matplotlib inline
%watermark -w
%watermark -iv

"""
Explanation: If any part of this notebook is used in your research, please cite with the reference found in README.md.

Connected components in a spatial network
Identifying and visualizing the parts of a network
Author: James D. Gaboardi &#106;&#103;&#97;&#98;&#111;&#97;&#114;&#100;&#105;&#64;&#103;&#109;&#97;&#105;&#108;&#46;&#99;&#111;&#109;
This notebook is a walk-through for:

Instantiating a simple network with libpysal.cg.Chain objects
Working with the network components and isolated rings
Visualizing the components and (non)articulation vertices
Longest vs. Largest components
Extracting network components
End of explanation
"""


plus1 = [
    Chain([Point([1, 2]), Point([0, 2])]),
    Chain([Point([1, 2]), Point([1, 1])]),
    Chain([Point([1, 2]), Point([1, 3])]),
]
plus2 = [
    Chain([Point([2, 1]), Point([2, 0])]),
    Chain([Point([2, 1]), Point([3, 1])]),
    Chain([Point([2, 1]), Point([2, 2])]),
]
lines = plus1 + plus2

ntw = spaghetti.Network(in_data=lines)

"""
Explanation: 1. Instantiate a network from two collections of libpysal.cg.Chain objects
End of explanation
"""


ntw.network_fully_connected

"""
Explanation: Here we get a warning because the network we created is not fully connected
End of explanation
"""


ntw.network_n_components

"""
Explanation: It has 2 connected components
End of explanation
"""


ntw.network_component_labels

"""
Explanation: The network components can be inspected through the following attributes
network_component_labels
End of explanation
"""


ntw.network_component2arc

"""
Explanation: network_component2arc
End of explanation
"""


ntw.network_component_lengths

"""
Explanation: network_component_lengths
End of explanation
"""


ntw.network_longest_component

"""
Explanation: network_longest_component
End of explanation
"""


ntw.network_component_vertices

"""
Explanation: network_component_vertices
End of explanation
"""


ntw.network_component_vertex_count

"""
Explanation: network_component_vertex_count
End of explanation
"""


ntw.network_largest_component

"""
Explanation: network_largest_component
End of explanation
"""


ntw.network_component_is_ring

"""
Explanation: network_component_is_ring
End of explanation
"""


ntw.graph_component_labels

"""
Explanation: The same can be performed for graph representations, for example:
graph_component_labels
End of explanation
"""


ntw.graph_component2edge

"""
Explanation: graph_component2edge
End of explanation
"""


# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)

"""
Explanation: Extract the network arc and vertices as geopandas.GeoDataFrame objects
End of explanation
"""


vertices_df

arcs_df

"""
Explanation: Network component labels are found in the "comp_label" field
End of explanation
"""


base = arcs_df.plot(column="comp_label", cmap="Set2", linewidth=5, figsize=(7, 7))
vertices_df.plot(ax=base, color="k", markersize=100, zorder=2);

"""
Explanation: Plot the disconnected network and symbolize the arcs bases on the value of "comp_label"
End of explanation
"""


new_lines = [
    Chain([Point([1, 1]), Point([2, 2])]),
    Chain([Point([0.5, 1]), Point([0.5, 0.5])]),
    Chain([Point([0.5, 0.5]), Point([1, 0.5])]),
    Chain([Point([2, 2.5]), Point([2.5, 2.5])]),
    Chain([Point([2.5, 2.5]), Point([2.5, 2])]),
]
lines += new_lines

ntw = spaghetti.Network(in_data=lines)

"""
Explanation: 2. Add to the network created above
End of explanation
"""


ntw.network_n_components

ntw.network_component2arc

# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)

arcs_df

"""
Explanation: Now there are 3 connected components in the network
End of explanation
"""


ntw.non_articulation_points

"""
Explanation: We can also inspect the non-articulation points in the network. Non-articulation points are vertices in a network that are degree-2. A vertex is degree-2 if, and only if, it is directly connected to only 2 other vertices.
End of explanation
"""


napts = ntw.non_articulation_points
articulation_vertices = vertices_df[~vertices_df["id"].isin(napts)]
non_articulation_vertices = vertices_df[vertices_df["id"].isin(napts)]

"""
Explanation: Slice out the articulation points and non-articulation points
End of explanation
"""


base = arcs_df.plot(column="comp_label", cmap="Set2", linewidth=5, figsize=(7, 7))
articulation_vertices.plot(ax=base, color="k", markersize=100, zorder=2)
non_articulation_vertices.plot(ax=base, marker="s", color="k", markersize=20, zorder=2);

"""
Explanation: Plot the connected components while making a distinction between articulation points and non-articulation points
End of explanation
"""


new_lines = [
    Chain([Point([3, 1]), Point([3.25, 1.25])]),
    Chain([Point([3.25, 1.25]), Point([3.5, 1.25])]),
    Chain([Point([3.5, 1.25]), Point([3.75, 1])]),
    Chain([Point([3.75, 1]), Point([3.5, 0.75])]),
    Chain([Point([3.5, 0.75]), Point([3.25, 0.75])]),
    Chain([Point([3.25, 0.75]), Point([3, 1])]),
]
lines += new_lines

ntw = spaghetti.Network(in_data=lines)

ntw.network_n_components

ntw.network_component2arc

# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)

arcs_df

"""
Explanation: 3. Add a loop of libpysal.cg.Chain objects
End of explanation
"""


ntw.non_articulation_points

"""
Explanation: Here we can see that all the new network vertices are non-articulation point
End of explanation
"""


napts = ntw.non_articulation_points
articulation_vertices = vertices_df[~vertices_df["id"].isin(napts)]
non_articulation_vertices = vertices_df[vertices_df["id"].isin(napts)]

"""
Explanation: Slice out the articulation points and non-articulation points
End of explanation
"""


base = arcs_df.plot(column="comp_label", cmap="Set2", linewidth=5, figsize=(7, 7))
articulation_vertices.plot(ax=base, color="k", markersize=100, zorder=2)
non_articulation_vertices.plot(ax=base, marker="s", color="k", markersize=20, zorder=2);

"""
Explanation: The new network vertices are non-articulation points because they form a closed ring
End of explanation
"""


cross = [
    Chain([Point([0, 5]), Point([5, 5]), Point([5, 10])]),
    Chain([Point([5, 0]), Point([5, 5]), Point([10, 5])]),
]
hexagon = [
    Chain(
        [
            Point([12, 5]),
            Point([13, 6]),
            Point([14, 6]),
            Point([15, 5]),
            Point([14, 4]),
            Point([13, 4]),
            Point([12, 5]),
        ]
    ),
]
lines = cross + hexagon

ntw = spaghetti.Network(in_data=lines)

# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)

base = arcs_df.plot(column="comp_label", cmap="Set2", linewidth=5, figsize=(7, 7))
vertices_df.plot(ax=base, color="k", markersize=100, zorder=2);

"""
Explanation: 4. Longest vs. largest components — cross vs. hexagon
End of explanation
"""


clongest = ntw.network_longest_component
clength = round(ntw.network_component_lengths[clongest], 5)
clargest = ntw.network_largest_component
cverts = ntw.network_component_vertex_count[clargest]
print("The longest component is %s at %s units of distance." % (clongest, clength))
print("The largest component is %s with %s vertices." % (clargest, cverts))

"""
Explanation: The longest component is not necessarily the largest
This is because in spaghetti the largest compnent equates to the most vertices
End of explanation
"""


longest = spaghetti.extract_component(ntw, ntw.network_longest_component)

# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(longest, vertices=True, arcs=True)

vertices_df

base = arcs_df.plot(column="comp_label", cmap="Set2", linewidth=5, figsize=(7, 7))
vertices_df.plot(ax=base, color="k", markersize=100, zorder=2);

"""
Explanation: 5. Extracting components
Extract the longest component
End of explanation
"""


largest = spaghetti.extract_component(ntw, ntw.network_largest_component)
# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(largest, vertices=True, arcs=True)
base = arcs_df.plot(column="comp_label", cmap="Set2", linewidth=5, figsize=(7, 7))
vertices_df.plot(ax=base, color="k", markersize=100, zorder=2);

"""
Explanation: Extract the largest component and plot
End of explanation
"""


libpysal.examples.available()
newhaven = libpysal.examples.get_path("newhaven_nework.shp")
ntw = spaghetti.Network(in_data=newhaven, extractgraph=False)

"""
Explanation: Empirical Example — New Haven, Connecticut
End of explanation
"""


longest = spaghetti.extract_component(ntw, ntw.network_longest_component)

# network vertices and arcs
vertices_df, arcs_df = spaghetti.element_as_gdf(ntw, vertices=True, arcs=True)
arcs_df.crs = "epsg:4269"
arcs_df = arcs_df.to_crs("epsg:6433")

# longest vertices and arcs
lc_vertices, lc_arcs = spaghetti.element_as_gdf(longest, vertices=True, arcs=True)
lc_arcs.crs = "epsg:4269"
lc_arcs = lc_arcs.to_crs("epsg:6433")

"""
Explanation: Extract the longest component
End of explanation
"""


nlc = ntw.network_longest_component
arcs_df = arcs_df[arcs_df.comp_label != nlc]
ocomp = list(set(ntw.network_component_labels))
ocomp.remove(nlc)

"""
Explanation: Filter non-longest component arcs
End of explanation
"""


def legend(objects):
    """Add a legend to a plot"""
    patches = make_patches(*objects)
    kws = {"fancybox": True, "framealpha": 0.85, "fontsize": "x-large"}
    kws.update({"loc": "lower left", "labelspacing": 2.0, "borderpad": 2.0})
    legend = matplotlib.pyplot.legend(handles=patches, **kws)
    legend.get_frame().set_facecolor("white")

def make_patches(comp_type, in_comp, oc):
    """Create patches for legend"""
    labels_colors_alpha = [
        ["%s component: %s" % (comp_type.capitalize(), in_comp), "k", 0.5],
        ["Other components: %s-%s" % (oc[0], oc[1]), "r", 1],
    ]
    patches = []
    for l, c, a in labels_colors_alpha:
        p = matplotlib.lines.Line2D([], [], lw=2, label=l, c=c, alpha=a)
        patches.append(p)
    return patches

base = arcs_df.plot(color="r", alpha=1, linewidth=3, figsize=(10, 10))
lc_arcs.plot(ax=base, color="k", linewidth=2, alpha=0.5, zorder=2)
# add legend
legend(("longest", nlc, (ocomp[0], ocomp[-1])))
# add scale bar
scalebar = ScaleBar(3, units="m", location="lower right")
base.add_artist(scalebar)
base.set(xticklabels=[], xticks=[], yticklabels=[], yticks=[]);

"""
Explanation: Plot network arcs
End of explanation
"""


