Pandana Network Store from Shapefile

generate_pdna_store_from_shp.py

# This GIST provides a template for conveting shapefiles into pandana networks. 
# Under a MIT Open Source License. https://opensource.org/licenses/MIT
# Researchers cite as Wasserman, D. Geopandas to Pandana Network Converter. (2019) GitHub repository, https://gist.github.com/d-wasserman/2d50671b37ee46b088e155293399a90c

def get_nodes_and_edges(shp_file,rounding=5):
    """Use geopandas to read line shapefile and compile all paths and nodes in a line file based on a rounding tolerance.
    shp_file:path to polyline file with end to end connectivity
    rounding: tolerance parameter for coordinate precision"""
    edges = gpd.read_file(shp_file)
    edges["from_x"]=edges["geometry"].apply(lambda x:round(x.coords[0][0],rounding))
    edges["from_y"]=edges["geometry"].apply(lambda x:round(x.coords[0][1],rounding))
    edges["to_x"]=edges["geometry"].apply(lambda x:round(x.coords[-1][0],rounding))
    edges["to_y"]=edges["geometry"].apply(lambda x:round(x.coords[-1][1],rounding))
    nodes_from = edges[["from_x","from_y"]].rename(index=str,columns={"from_x":"x","from_y":"y"})
    nodes_to = edges[["to_x","to_y"]].rename(index=str,columns={"to_x":"x","to_y":"y"})
    nodes = pd.concat([nodes_from,nodes_to],axis=0)
    nodes["xy"] = list(zip(nodes["x"], nodes["y"]))
    nodes = pd.DataFrame(nodes["xy"].unique(),columns=["xy"])
    nodes["x"] = nodes["xy"].apply(lambda x: x[0])
    nodes["y"] = nodes["xy"].apply(lambda x: x[1])
    nodes = nodes[["x","y"]].copy()
    return nodes , edges

def generate_pandana_store_from_shp(hdf5_path, shp_file, weights=["weight"], oneway=None, overwrite_existing=True,
                                    rounding=6):
    """Generate a pandana ready HDF5 store using geopandas (gdal required) and pandas. Function assumes shape file has
    the appropriate topology where each edge has end points that lead to nodes. If using OSM, pre-process with Osmnx.
     Python 3.5.
        Parameters
        ----------
        hdf5_path: str  
            output path of HDF5 store holding two dataframes ["nodes","edges"]
        shp_file: str 
            input file that geopandas reads to make a graph based on end to end connectivity
        weights : lists
            weights columns transfered to the store edges. Name is maintained.
        oneway: str
            series name where oneway streets (edges) are denoted with a 1, 0 denotes twoway. None, assumes twoway edge.
        overwrite_existing: boolean
            if true, the existing store is overwritten.
        rounding:int  
            the number of digits to round line coordinates to get unique nodes (precision)
        Returns
        -----------
            hdf5_path:basestring - path of output hdf5 pandana store."""
    if os.path.exists(hdf5_path):
        if overwrite_existing:
            print("Overwriting existing store...")
            os.remove(hdf5_path)
        else:
            print("Existing store at path: {0}".format(hdf5_path))
            return hdf5_path
    all_edges_twoway = True
    oneway_field_list = []
    if oneway is not None:
        all_edges_twoway = False
        oneway_field_list.append(oneway)
    print("Reading shapefile with geopandas: {0}...".format(shp_file))
    nodes, edges = get_nodes_and_edges(shp_file, rounding)
    h5store = pd.HDFStore(hdf5_path)
    print("Establishing node store...")
    df_nodes = nodes
    df_nodes["id"] = df_nodes.index.values
    df_nodes.index.rename("id", True)
    h5store['nodes'] = df_nodes
    edge_cnt = len(edges)
    print("Establishing edge store for {0} edges...".format(edge_cnt))
    df_edges = edges[['from_x', 'from_y', 'to_x', 'to_y'] + weights + oneway_field_list].copy()
    print("Transferring nodeids to edges...")
    df_edges = pd.merge(df_edges, df_nodes, how='left', left_on=['from_x', 'from_y'], right_on=['x', 'y'])
    df_edges = pd.merge(df_edges, df_nodes, how='left', left_on=['to_x', 'to_y'], right_on=['x', 'y'],
                        suffixes=('_from', '_to'))
    # nodeids are duplicated on from the joined nodes, joined first to from, suffix to on next set
    df_edges.rename(columns={'id_from': 'from', 'id_to': 'to'}, inplace=True)
    df_edges = df_edges[['from', 'to'] + weights + oneway_field_list]
    if all_edges_twoway:
        print("""Note: Edges are duplicated in this step, do not use the 'twoway' setting in the pandana network if using this 
             function""")
        twoway_edges = df_edges.copy()
        twoway_to = twoway_edges["to"].copy()
        twoway_edges["to"] = twoway_edges["from"]
        twoway_edges["from"] = twoway_to
        df_edges = pd.concat([df_edges, twoway_edges])
    else:
        print("Setting up edges based on oneway field...")
        twoway_edges = df_edges[df_edges[oneway] == 0].copy()
        twoway_to = twoway_edges["to"].copy()
        twoway_edges["to"] = twoway_edges["from"]
        twoway_edges["from"] = twoway_to
        df_edges = pd.concat([df_edges, twoway_edges])
    h5store['edges'] = df_edges
    h5store.close()
    print("Graph store construction complete...")
    return hdf5_path

Thanks for sharing.

This is a code that should be integrated to the Pandana API.

Thanks for sharing.

This is a code that should be integrated to the Pandana API.

I won't argue against that! It seems they have internal tools for it that they have not exposed for some reason.

Hey David, I tested your work, but it seems to be some nodes isolated from the network.

The examples bellow uses the same bbox to create accessibility with 3 schools. The first one uses the pandana API loader, the second uses your code. Am i doing something wrong?

I calculated the distance outside the api, converting using geopandas to Portugal epsg.

networkGeoDataFrame = networkGeoDataFrame.to_crs({'init':'epsg:3763'})
networkGeoDataFrame['dist'] = networkGeoDataFrame['geometry'].length

After the conversion, I used the distance with the raw file from OSM.

Per our LinkedIn conversation, the function assumes you are using a file that has end to end connectivity. This means all edges that must connect, must share an end point. Pandana's native OSM importer uses OSRM to establish this connectivity, where as my function assumes you are working from a file that does this ahead of time.
Some notes:

• Be careful not to planarize segments that are on highways, bridges, ramps etc. You have to check that directly.
• Check that the impedance are similar between similar edges in OSM or your network. It might be edge "cost" related.

David,

I solved this connectivity problem with graph_from_polygon method from OSMnx. It need some tinkering to get the network from every polygon, But it saved me from the connection verification.

I will leave the link here for future reference.

Hi @fillipefeitosa,

I am glad that worked! Can you share back with me how you used it to create a panda network from OSMnx? I want to try to see if I can recreate the nodes in a similar fashion or at least as an option.

Hey David!

Hope everything is all right you with and your family. Things got kinda crazy here with covid pandemic.

The class I refer above is this:

First I get the municipalities from a shapefile using dissolve.
After getting the network using OSMnx, I save the shp file on a temporary folder and convert it to pandanas using your code.
Slicing the entire country in smaller chunks and using multiprocessing, it was possible to create pandana networks for the entire country. For later use, just load the correspondent HDF5

Notebook

def create_network(polygons_dataframe, desired_network, dissolve_column='Concelho'):
    """Create HDF5 with street network graph from desired polygons from a 
    large dataframe from different places. A good use case would be a Geopandas dataframe
    with polygons from a region, or a district.
    Parameters
    ----
    dissolve_column:The column wicth will be used to fecth the results.
    desired_network:The code, name, or whatever is stored in dissolve_column.
    polygons_dataframe:Geopandas Dataframe with polygons data.
    Returns:
    ----
    None
    """
    polygonAsSeries = polygons_dataframe[polygons_dataframe[dissolve_column]==desired_network]
    if(len(polygonAsSeries)==1):
        polygonAsSeries = polygonAsSeries.squeeze()
    elif(len(polygonAsSeries)>1):
        polygons_dataframe = polygons_dataframe.dissolve(by=dissolve_column, as_index=False)
        polygonAsSeries = polygons_dataframe[polygons_dataframe[dissolve_column]==desired_network].squeeze()
        
    latlng_geom, _ = ox.project_geometry(polygonAsSeries.geometry, crs={'init':'epsg:3763'}, to_latlong=True)
    network = ox.graph_from_polygon(latlng_geom, network_type='all_private')
    # Save Shapefile
    try: 
        ox.save_graph_shapefile(network, filename=desired_network, folder='./data/temp/')
    except:
        raise Exception("It was not possible to save %s network as shapefile.")
    else:
        # Load and Convert Shapefile
        save_path = './data/networks/'+desired_network
        load_shp_path = './data/temp/'+desired_network+'/edges/edges.shp'
        path = generate_pandana_store_from_shp(save_path, load_shp_path, weights=['dist'])
        # Removing useless shp Files
        shutil.rmtree('./data/temp/'+desired_network, ignore_errors=True)

This makes sense! Thanks for sharing. I think I need to make a version just to take osmnx edges, but that should be straight forward.