<jupyter_start><jupyter_text>Step 3 - Snap CXB OD Estimates This notebook is the final step of calculating the base OD matrices which underlie all GOSTNets-derived accessibility analysis. It's also the slowest step<jupyter_code>import geopandas as gpd
import pandas as pd
import os, sys, time, importlib

# add to your system path the location of the LoadOSM.py and GOSTnet.py scripts
# sys.path.append("../../../GOSTnets/GOSTnets")

import GOSTnets as gn
import importlib
# import Network_Clean as gnClean
importlib.reload(gn)

import networkx as nx
import osmnx
from shapely.geometry import Point
import numpy as np

import rasterio
from rasterio import features
from shapely.wkt import loads

import ipyparallel as ipp<jupyter_output><empty_output><jupyter_text>Setup paths<jupyter_code>input_pth = r'inputs\\dests'
geo_pth = r'../../../GEO'
fin_pth = 'final'<jupyter_output><empty_output><jupyter_text>pickles<jupyter_code>fin_pckle = r'final_current_G.pickle'
upgr_all = r'final_upgrade_all_G.pickle'
upgr_nosouth = r'final_upgrade_nosouth_G.pickle'
upgr_noferry = r'final_upgrade_noferry_G.pickle'
G_current = nx.read_gpickle(os.path.join(fin_pth, fin_pckle))
G_upgr_all = nx.read_gpickle(os.path.join(fin_pth, upgr_all))
G_upgr_nosouth = nx.read_gpickle(os.path.join(fin_pth, upgr_nosouth))
G_upgr_noferry = nx.read_gpickle(os.path.join(fin_pth, upgr_noferry))<jupyter_output><empty_output><jupyter_text>origins<jupyter_code># grid_name = r'hrsl_2018_cxb_pts_snapped.csv'
grid_name = r'growth_center_origins_snapped.csv'
origins = pd.read_csv(os.path.join(fin_pth,grid_name))
grid = origins<jupyter_output><empty_output><jupyter_text>destinations<jupyter_code>cxb_fil = r'cxb_ctr.shp'
chitt_fil = r'chittagong.shp'
health_fil = r'hc_merge_200324_4326.shp'
primary_fil = r'schools/school_category_primary.gpkg'
secondary_fil = r'schools/school_category_secondary.gpkg'
tertiary_fil = r'schools/school_category_tertiary.gpkg'
matar_fil = r'martarbari.shp'
mkts_fil = r'mkts_merge_4326.shp'
gc_fil = r'cxb_lged_gc_moved_4326.shp'
dests = {"CXB" : cxb_fil, "Chittagong" : chitt_fil, "Health" : health_fil, \
         "Primary_education" : primary_fil, "Secondary_education" : secondary_fil, "Tertiary_education" : tertiary_fil, \
         "Martarbari" : matar_fil, "All_markets" : mkts_fil, "Growth_centers" : gc_fil}

# note you can use smaller / larger dest dictionaries as needed
# This is helpful for going back and re-running only certain destinations, or adding in new ones.

# dests = {"CXB" : cxb_fil}<jupyter_output><empty_output><jupyter_text>scenarios (for looping)<jupyter_code># the dict here allows us to match a label to a pickle

scenarios = {'current' : [G_current,'_current_snapped.csv'],\
             'upgrade_all' : [G_upgr_all,'_ua_snapped.csv'],\
             'upgrade_nosouth' : [G_upgr_nosouth,'_uns_snapped.csv'],\
             'upgrade_noferry' : [G_upgr_noferry,'_unf_snapped.csv']}<jupyter_output><empty_output><jupyter_text>Settings<jupyter_code>walk_speed = 4.5
WGS = {'epsg':'4326'}
measure_crs = {'epsg':'32646'}
# date = 'May2020'<jupyter_output><empty_output><jupyter_text>Troubleshooting<jupyter_code># dests = {"All_markets" : mkts_fil, "Growth_centers" : gc_fil}
# cur = {'current' : [G_current,'_current_snapped.csv']}

# for scen, values in cur.items():
# #     namestr(values, globals())
#     od_routine(values[0],scen,values[1])

# for scen, values in cur.items():
#     add_walking_time(scen,values[1])
# for u, v, data in G_upgr_all.edges(data=True):
#     if data['osm_id'] == '244861071':
#         print(data['infra_type'])
# for u, v, data in G_upgr_nosouth.edges(data=True):
#     if data['osm_id'] == '244861071':
#         print(data['infra_type'])<jupyter_output><empty_output><jupyter_text>Define functions that will generate OD matrices, and then generate direct walking times from origins to destinations The benefit of functions is they allow us to loop the analysis over the dictionary items These functions are big, messy, and adapted from old GOSTNets code. Though functional, the code could likely be streamlined as it's tricky to troubleshoot.<jupyter_code>def od_routine(G_input,scenario='',snap_ending='_snapped.csv'):
    
    origins_pop_nodes = list(set(origins.NN)) # consolidates by shared nodes. 
    
    for dest_type, fpth in dests.items():

        snapfile = dest_type + snap_ending

        print(dest_type)

        dest = pd.read_csv(os.path.join(fin_pth,snapfile))
        dest_nodes = list(set(dest.NN))

        print(len(list(set(dest.NN))))

        od_time = gn.calculate_OD(G_input, origins=origins_pop_nodes, 
                              destinations=dest_nodes, fail_value=99999999, weight='time')

        od_time_df = pd.DataFrame(od_time, index=origins_pop_nodes, columns=dest_nodes)

        print(od_time_df.shape)

        # Add walking time (from origin to NN) for each OD

        # origins_join = origins_pop_snapped.merge(od_time_df, how='left', on='NN')
        origins['NN_dist_seconds'] = ((origins.NN_dist / 1000) / walk_speed) * 60 * 60
        origins_join = origins.join(od_time_df, on='NN', rsuffix="dist_")

    #     print(origins_join.head())

        origins_join.columns[6:len(origins_join.columns)]

        origins_join.to_csv(os.path.join(fin_pth,'origins_walktime_{}_NN_{}.csv'.format(scenario,dest_type)))
        od_time_df.to_csv(os.path.join(fin_pth,'OD_matrix_{}_NN_{}.csv'.format(scenario,dest_type)))
def add_walking_time(scenario='',snap_ending='_snapped.csv'):
    # Main method
    
    print(scenario)

    for dest_type, fpth in dests.items():

        snapfile = dest_type + snap_ending

        print(dest_type)

        OD_name = r'OD_matrix_{}_NN_{}.csv'.format(scenario,dest_type)

        OD = pd.read_csv(os.path.join(fin_pth, OD_name))
        # OD = od_time_df

        OD = OD.rename(columns = {'Unnamed: 0':'O_ID'})
        OD = OD.set_index('O_ID')
        OD = OD.replace([np.inf, -np.inf], np.nan)

    #     # Filtering by only desired destination in an all-destination OD matrix.
    #     # Skipping for now

        od_dest_df = pd.read_csv(os.path.join(fin_pth,snapfile))
        od_dest_df['geometry'] = od_dest_df['geometry'].apply(loads)
#         od_dest_gdf = gpd.GeoDataFrame(od_dest_df, crs = {'init':'epsg:4326'}, geometry = 'geometry')

    #     accepted_facilities = list(set(list(acceptable_df.NN)))
    #     accepted_facilities_str = [str(i) for i in accepted_facilities]

    #     print(accepted_facilities)
    #     print(accepted_facilities_str)
    #     # OD = OD_original[accepted_facilities_str] # not necessary, already done
    #     # acceptable_df.to_csv(os.path.join(basepth,'Output','%s.csv' % subset))

        # Computing walk time from network to destination

        dest = pd.read_csv(os.path.join(fin_pth,snapfile))

        dest_df = dest[['NN','NN_dist']]
        dest_df = dest_df.set_index('NN')

        dest_df['NN_dist'] = dest_df['NN_dist'] / 1000 * 3600 / walk_speed
        dest_df.index = dest_df.index.map(str)

        d_f = OD.transpose()

        for i in d_f.columns:
            dest_df[i] = d_f[i]

        for i in dest_df.columns:
            if i == 'NN_dist':
                pass
            else:
                dest_df[i] = dest_df[i] + dest_df['NN_dist']
                
        dest_df = dest_df.drop('NN_dist', axis = 1)
        dest_df = dest_df.transpose()
        dest_df['min_time'] = dest_df.min(axis = 1)

#         dest_df['geometry'] = dest_df['geometry'].apply(loads)
        dest_gdf = gpd.GeoDataFrame(od_dest_df, geometry = 'geometry', crs = {'init':'epsg:4326'})

        # Add walk time from origin to network

        grid = pd.read_csv(os.path.join(fin_pth, grid_name))
        grid = grid.rename(columns = {'NN':'O_ID','NN_dist':'walk_to_road_net'})
        grid = grid.set_index(grid['O_ID'])

        grid['on_network_time'] = dest_df['min_time']
        grid['walk_to_road_net'] = grid['walk_to_road_net'] / 1000 * 3600 / walk_speed 
        grid['total_time_net'] = grid['on_network_time'] + grid['walk_to_road_net']

#         print(grid.head())

        grid['geometry'] = grid['geometry'].apply(loads)
        o_2_d = gpd.GeoDataFrame(grid, crs = {'init':'epsg:4326'}, geometry = 'geometry')

            # Snapping!

        print('start of snapping: %s\n' % time.ctime())
        o_2_d = gn.pandana_snap_points(o_2_d, 
                                   dest_gdf, # eventually just use dest_gdf
                                   source_crs='epsg:4326',
                                   target_crs='epsg:32646',
                                   add_dist_to_node_col = True)
        print('\nend of snapping: %s' % time.ctime())
        print('\n--- processing complete')

            # Recalculating the resulting walking times into seconds and minutes.
            # Make sure that if walking is faster than on-network travel, it prefers walking

        o_2_d['walk_time_direct'] = o_2_d['idx_dist'] / 1000 * 3600 / walk_speed

        grid['walk_time_direct'] = o_2_d['walk_time_direct']

        # grid['PLOT_TIME_SECS'] = grid[['total_time_net']].min(axis = 1)
        # grid['PLOT_TIME_SECS'] = grid[['walk_to_road_net','total_time_net']].min(axis = 1)

            # The city locations / port location don't have walk_time_direct values so we use if/else logic to work around them.

        if 'walk_time_direct' in grid.columns:
            grid['PLOT_TIME_SECS'] = grid[['walk_time_direct','total_time_net']].min(axis = 1)
        else:
            grid['PLOT_TIME_SECS'] = grid[['total_time_net']]

        grid['PLOT_TIME_MINS'] = grid['PLOT_TIME_SECS'] / 60

        if 'walk_time_direct' in grid.columns:
            def choice(x):
                if x.walk_time_direct < x.total_time_net:
                    return 'walk'
                else:
                    return 'net'

            grid['choice'] = grid.apply(lambda x: choice(x), axis = 1)
            grid['choice'].value_counts()

#         print(grid.head())

        # Export

        grid.to_csv(os.path.join(fin_pth,'final_cxb_{}_od_grid_{}.csv'.format(scenario,dest_type)))<jupyter_output><empty_output><jupyter_text>Load origins and destinations, get unique origin nodes, run OD, export Run the OD routine function on all destinations for the 3 different scenarios<jupyter_code>for scen, values in scenarios.items():
    od_routine(values[0],scen,values[1])<jupyter_output>CXB
1
(34, 1)
CXB
1
(34, 1)
CXB
1
(34, 1)
CXB
1
(34, 1)<jupyter_text>Import completed OD matrix, calcualte walking times from origins to the destinations<jupyter_code>for scen, values in scenarios.items():
    add_walking_time(scen,values[1])<jupyter_output>current
CXB
start of snapping: Fri Jun 12 14:54:38 2020
{'fraction_non_alphanumeric': 0.1234948193783254, 'fraction_numerical': 0.0156818818258191, 'mean_word_length': 3.553336166595835, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '51905', 'license': 'MIT', 'n_tokens_mistral': 4032, 'n_tokens_neox': 3782, 'n_words': 726, 'path': 'Implementations/FY21/ACC_Bangladesh_CXB_Accessibility_and_Growth/Step 3 - CXB OD Estimates.ipynb', 'repo_name': 'worldbank/GOST_PublicGoods', 'size': '10713'}