How to do angular closeness centrality from a shapefile

[felipetavares80]

Hello everyone.

I'm trying to do a angular weighted closeness centrality from a segmented road map.

I'm importing it from osmnx, simplifying the segments by Douglas-Peucker on Qgis, extracting the points and creating a graph for momepy closeness centrality analysis.

I did the graph but it seems that the results do not match the manual calculations when I set the angular option for closeness centrality. It matches for topological but not for angular.

I don't know if I'm doing something wrong on manual calculations or in the coding, but 'll appreciate a help on this.

Thanks in advance!

All the Best.

The Codes, the manual calculations document and the shapefiles is attached.
angularclosenesspack.zip

[martinfleis]

Hi @felipetavares80,

First question - which version of momepy do you have? There was a bug in creation of angular graph recently (fixed in the latest release).

Second - what makes you think it matches topological calculation? See

import geopandas as gpd
import osmnx as ox
import networkx as nx
import momepy

# load Shapefiles as node/edge GeoDataFrames indexed as described in OSMnx docs
gdf_nodes = gpd.read_file('./teste/nodes.shp').set_index('osmid')
gdf_edges = gpd.read_file('./teste/edges.shp').set_index(['u', 'v', 'key'])

# Extract columns x e y coordinates from points
gdf_nodes['x'] = gdf_nodes.geometry.x
gdf_nodes['y'] = gdf_nodes.geometry.y

# convert the node/edge GeoDataFrames to a MultiDiGraph
graph_attrs = {'crs': 'epsg:32625', 'simplified': True}

#G2 = ox.graph_from_gdfs(gdf_nodes, gdf_edges)
G2 = ox.graph_from_gdfs(gdf_nodes, gdf_edges, graph_attrs)



# saving a gdf and letting momepy do the conversion back to the graph.
edges = ox.graph_to_gdfs(G2, nodes=False, edges=True,
                                   node_geometry=False, fill_edge_geometry=True)

dual = momepy.gdf_to_nx(edges, approach='dual')

# measure with angle as a weight
dual = momepy.closeness_centrality(dual, weight='angle')
# measure without weight - topological
dual = momepy.closeness_centrality(dual, name="topological")

dual_gdf = momepy.nx_to_gdf(dual, points=False)

print(dual_gdf[["closeness", "topological"]].head())

You can see that the values are different

   closeness  topological
0   0.004002     0.077749
1   0.005351     0.083963
2   0.004401     0.076200
3   0.005510     0.089063
4   0.006022     0.086391

Which is also shown on plots:

dual_gdf.plot("closeness")

dual_gdf.plot("topological")

It seems to work as expected on my side, so I may need a bit more info to see the issue.

I did the graph but it seems that the results do not match the manual calculations when I set the angular option for closeness centrality. It matches for topological but not for angular.

How do you do the manual calculation?

[martinfleis]

It may actually be the best to depend on cityseer directly in this. There's no need to re-implement the algorithm and maintain it at two places. I'm just reluctant to expose yet another object (cityseer.structures.NetworkStructure) momepy users need to deal with... It would certainly break the logic we have now, of having a pure networkx graph for network analysis and GeoDataFrame for geometry analysis.

@songololo can you think about a shortcut converting a networkx object created by momepy.nx_from_gdf to a structure needed by city seer and retrieving the values from it so we can attach them as node/edge attributes? That way the signature of momepy's centrality measures does not have to change and cityseer object would not be exposed here.

[songololo]

Will try to look at this during the next week or so.

[songololo]

For discussion, here is how this would roughly look via cityseer. I'm not using osmnx for conversion to a GeoDataframe because it seems to use directed networks, which are not compatible with cityseer. Once computed the columns could be copied back and forth between dataframes?

import geopandas as gpd
import networkx as nx
from cityseer.metrics import networks
from cityseer.tools import graphs, plot

# cityseer uses undirected networks - assumes pedestrians
G1 = nx.MultiGraph()
# add nodes
gdf_nodes = gpd.read_file("./temp/nodes.shp").set_index("osmid")
for idx, row in gdf_nodes.iterrows():
    G1.add_node(idx, x=row["geometry"].x, y=row["geometry"].y)
# add edges
gdf_edges = gpd.read_file("./temp/edges.shp").set_index(["u", "v", "key"])
for (idx_a, idx_b, edge_idx), row in gdf_edges.iterrows():
    G1.add_edge(idx_a, idx_b, edge_idx, geom=row["geometry"])
# visual check
plot.plot_nx(G1, labels=True)

# calculating angular centrality over the primal representation
nodes_gdf, network_structure = graphs.network_structure_from_nx(G1, crs=32625)
# network structure contains information for jitted functions / numpy like
nodes_gdf = networks.node_centrality(
    measures=["node_harmonic_angular"],
    network_structure=network_structure,
    nodes_gdf=nodes_gdf,
    distances=[100, 400],
    angular=True,
)
print(nodes_gdf.head())

# using the dual representation
G1_dual = graphs.nx_to_dual(G1)
# use plot_geoms to show geometric representation
plot.plot_nx(G1_dual, labels=True, plot_geoms=True)
nodes_gdf_dual, network_structure_dual = graphs.network_structure_from_nx(
    G1_dual, crs=32625
)
nodes_gdf_dual = networks.node_centrality(
    measures=["node_harmonic_angular"],
    network_structure=network_structure_dual,
    nodes_gdf=nodes_gdf_dual,
    distances=[100, 400],
    angular=True,
)
print(nodes_gdf_dual.head())

[songololo]

It would be possible to create a convenience method for going directly from momepy GeoDataframes to cityseer data structures and back again. i.e. to skip the networkX step.

[songololo]

For info - the network_structure in cityseer is basically a data structure containing numpy arrays which are used by the jitted functions (so that it can scale to larger graphs). EDIT: This also pre-computes information such as in-out bearings for segments and segment impedances (distance and angular change), etc.