<jupyter_start><jupyter_text>Step 5 - Gravity Modeling The gravity modeling process can be broken into two primary sections -- one in which the OD matrices are manipulated into shapes compatiable with the gravity modeling tools, and one in which the actual gravity models are developed and exported for origins and destinations. This notebook covers both.<jupyter_code>import os, sys
import time
import pandas as pd
import geopandas as gpd
import numpy as np

import networkx as nx
from scipy import spatial
from functools import partial
import pyproj

import shapely
from shapely.ops import transform
from shapely import wkt
from shapely.geometry import Point, MultiPoint
from shapely.wkt import loads

import GOSTnets as gn
import geopy<jupyter_output><empty_output><jupyter_text>Define stuff Fixed parameters<jupyter_code>simplif_meters = 25
source_epsg = 4326
target_epsg = 3106<jupyter_output><empty_output><jupyter_text>Variable parameters<jupyter_code>scenario = 'Current'
# scenario = 'Padma'
# WorldPop data determinants

constraint_status = 'constrained'
# constraint_status = 'unconstrained'

wp_res = 100
# wp_res = 250
# wp_res = '1k'
# Production date for outputs being used

prod_date = '210329'
origins = 'PopOrigins'
# origins = 'CityOrigins'
# origins = 'All_SEZ_Origins'
# origins = 'Active_SEZ_Origins'<jupyter_output><empty_output><jupyter_text>Path names and file names<jupyter_code>input_pth = r'inputs\\dests'
interm_pth = r'intermediate'
fin_pth = r'final'
res_pth = r'results'

tab_pth = r'../../../Tabular'
geo_pth = r'../../../GEO'
origin_folder = r'..\..\..\GEO\Population'<jupyter_output><empty_output><jupyter_text>Settings and functions<jupyter_code># pd.set_option('max_columns',None)
# the gravity model throws error messages that reduce readibility, let's get rid of them

import warnings

def fxn():
    warnings.warn("deprecated", DeprecationWarning)

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fxn()
# function for sorting alphanumerically

import re

def sorted_nicely( l ): 
    """ Sort the given iterable in the way that humans expect.""" 
    convert = lambda text: int(text) if text.isdigit() else text 
    alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ] 
    return sorted(l, key = alphanum_key)

# funciton for sorting matrices smallest to largest, by origin ID then destination ID

def sort_od_matrix(od_matrix):
    
    # sort by O_IDs, then dest node IDs
    od_matrix = od_matrix.sort_values('Unnamed: 0').reindex(sorted_nicely(od_matrix.columns), axis=1)

#     # reset O_ID column to the front
#     od_matrix = od_matrix[ ['Unnamed: 0'] + [ col for col in od_matrix.columns if col != 'Unnamed: 0' ] ]

    # set the Dest_ID column back to index so the shape is the same as the dWeight shape
    od_matrix.set_index('Unnamed: 0',inplace=True)
# spatial join admin file to origin file based on within location

def georef_origins(orig_fil,adm_fil):
    orig_fil['geometry'] = orig_fil['geometry'].apply(wkt.loads)
    orig_fil = gpd.GeoDataFrame(orig_fil,geometry='geometry')
    orig_fil = gpd.sjoin(orig_fil,adm_fil,how='left',op='within')
# Optimization functions to slim the file sizes of our many large input files and thus speed up processing / export

from typing import List

def optimize_floats(df: pd.DataFrame) -> pd.DataFrame:
    floats = df.select_dtypes(include=['float64']).columns.tolist()
    df[floats] = df[floats].apply(pd.to_numeric, downcast='float')
    return df

def optimize_ints(df: pd.DataFrame) -> pd.DataFrame:
    ints = df.select_dtypes(include=['int64']).columns.tolist()
    df[ints] = df[ints].apply(pd.to_numeric, downcast='integer')
    return df

def optimize(df: pd.DataFrame, datetime_features: List[str] = []):
    return optimize_floats(optimize_ints(df))<jupyter_output><empty_output><jupyter_text>Data Preparation Create a dictionary of destination DFs<jupyter_code># Destination DFs

City_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'All_cities_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
Dhaka_Chitt_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'Dhaka_Chitt_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
Minor_cities_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'Minor_cities_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
Dry_ports_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'Dry_ports_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
River_ports_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'River_ports_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
Deep_sea_ports_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'Deep_sea_ports_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
All_SEZs_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'All_SEZs_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))
Active_SEZs_df = pd.read_csv(os.path.join(fin_pth,prod_date,\
                        f'Active_SEZs_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification_snapped.csv'))

# Dest dictionary

dests = {"All_cities" : City_df, "Dhaka_Chitt" : Dhaka_Chitt_df, "Minor_cities" : Minor_cities_df, \
         "Dry_ports" : Dry_ports_df, "River_ports" : River_ports_df, "Deep_sea_ports" : Deep_sea_ports_df, \
         "All_SEZs" : All_SEZs_df, "Active_SEZs" : Active_SEZs_df}<jupyter_output><empty_output><jupyter_text>Admin data<jupyter_code># load admin spatial data

adm2 = gpd.read_file(os.path.join(geo_pth,'Boundaries/bd_adm_bbs.gpkg'),layer='bd_adm2_wgs84')
adm2.crs = 'epsg:4326'

adm3 = gpd.read_file(os.path.join(geo_pth,'Boundaries/bd_adm_bbs.gpkg'),layer='bd_adm3_wgs84')
adm3.crs = 'epsg:4326'

adm4 = gpd.read_file(os.path.join(geo_pth,'Boundaries/bd_adm_bbs.gpkg'),layer='bd_adm4_wgs84')
adm4.crs = 'epsg:4326'

adm5 = gpd.read_file(os.path.join(geo_pth,'Boundaries/bd_adm_bbs.gpkg'),layer='bd_adm5_wgs84')
adm5.crs = 'epsg:4326'
# Fix types for later joining
adm3['adm3_pcode'] = adm3['adm3_pcode'].astype(str)
adm4['adm4_pcode'] = adm4['adm4_pcode'].astype(str)

# Fix names for later display
adm3['adm3_en'] = adm3['adm3_en'].str.title()<jupyter_output><empty_output><jupyter_text>Data prep set an index on the NN, reindex by an alphanumerically sorted list of IDs, then reset_index to return the NNs to the data frame<jupyter_code>for dest, df in dests.items():
    df[['Unweighted']] = 1
    df = df.sort_values('NN').set_index('NN').reset_index() # sort for correct joining, set/reset index to move NN column to the front for easy checking
    dests.update({dest:df})<jupyter_output><empty_output><jupyter_text>Check results, especially for IDs for use in gravity modeling<jupyter_code># check first key of dests dict
dests[list(dests.keys())[0]].head()<jupyter_output><empty_output><jupyter_text>OD matrix data prep  Prepare OD files<jupyter_code># origin to destination
orig_raw = gpd.read_file(os.path.join(res_pth,prod_date,'spatial',f'{scenario}_origins_{wp_res}m_TTs.gpkg'),layer=f"{origins}",driver='GPKG')
orig_raw = optimize(orig_raw)
orig_raw = orig_raw.sort_values('O_ID')
orig_raw.rename({f'{scenario}_{origins}_DhakaChitt_mins' : f'{scenario}_{origins}_Dhaka_Chitt_mins', \
                 f'{scenario}_{origins}_deep_ports_mins' : f'{scenario}_{origins}_deep_sea_ports_mins'},axis=1,inplace=True)

# Summarize origin populations by origin node, to enable appropriate population weighting
# population assigned to each origin node

if origins == 'PopOrigins':
    raw_orig = pd.read_csv(os.path.join(fin_pth,prod_date,f'bgd_wp_{constraint_status}_origins_{wp_res}m_2020_snapped_25m.csv'))
    raw_orig.rename({'NN':'O_ID','VALUE':'Pop'},axis=1,inplace=True)
    orig_wtd = pd.pivot_table(raw_orig,values='Pop',index='O_ID',aggfunc='sum')
    orig_wtd = optimize(orig_wtd)

elif origins == 'All_SEZ_Origins':
    raw_orig = pd.read_csv(os.path.join(fin_pth,prod_date,f'All_SEZs_{constraint_status}_{wp_res}m_res_25m_simplification_snapped.csv'))
    raw_orig.rename({'NN':'O_ID','VALUE':'Pop'},axis=1,inplace=True)
    orig_wtd = pd.pivot_table(raw_orig,values='All_SEZ_count',index='O_ID',aggfunc='sum')
    orig_wtd = optimize(orig_wtd)
    
elif origins == 'Active_SEZ_Origins':
    raw_orig = pd.read_csv(os.path.join(fin_pth,prod_date,f'Active_SEZs_{constraint_status}_{wp_res}m_res_25m_simplification_snapped.csv'))
    raw_orig.rename({'NN':'O_ID','VALUE':'Pop'},axis=1,inplace=True)
    orig_wtd = pd.pivot_table(raw_orig,values='Active_SEZ_count',index='O_ID',aggfunc='sum')
    orig_wtd = optimize(orig_wtd)
orig_wtd.head()<jupyter_output><empty_output><jupyter_text>Populate a dictionary with dests as keys and a list containing the vehicle ODM and a filtered raw origins file as the value<jupyter_code>orig_raw.columns
odm_dct= {}

for dest, dest_gdf in dests.items():
    
    # filter origins to just relevant information for this destination
    
    if origins == 'PopOrigins':
        dest_orig = orig_raw[['geometry','O_ID','Row_ID','VALUE',f'{scenario}_{origins}_{dest}_mins','adm2_pcode','adm3_pcode','adm4_pcode','adm5_pcode','adm2_pop','adm3_pop','adm4_pop','adm5_pop']]
    else:
        dest_orig = orig_raw
        
    dest_orig = dest_orig.set_index('O_ID')
    
    # Read in raw OD grid
    
    dest_odm = pd.read_csv(os.path.join(fin_pth,prod_date,f'OD_matrix_{scenario}_{origins}_to_{dest}_{constraint_status}_{wp_res}m_res_{simplif_meters}m_simplification.csv'))
    dest_odm = optimize(dest_odm)
    
    # Put the origins and destinations in alphanumeric order by node so that they line up and process correctly within the gravity model. Sort by O_IDs, then dest node IDs
    # Should work as a function but the function is returning unsorted data -- while raw code works. Curious?
    
    dest_odm = dest_odm.sort_values('Unnamed: 0').reindex(sorted_nicely(dest_odm.columns), axis=1)

    # set the O_ID column back to index so the shape is the same as the dWeight shape
    dest_odm.set_index('Unnamed: 0',inplace=True)
    
    # populate a dictionary to use in gravity modeling
    
    odm_dct.update({dest:[dest_orig,dest_odm,dest_gdf]})<jupyter_output><empty_output><jupyter_text>Gravity models<jupyter_code>import GOSTnets.calculate_od_raw as calcOD
# calcOD.calculate_gravity??<jupyter_output><empty_output><jupyter_text>Check out weighting options<jupyter_code>rename_decay_dct = {
    'd_0.01' : '2.9mins_decay',
    'd_0.005' : '5.8mins_decay',
    'd_0.001' : '11.5mins_decay',
    'd_0.0007701635' :'15mins_decay',
    'd_0.0003850818' :'30mins_decay',
    'd_0.0001925409' :'60mins_decay',
    'd_9.62704e-05' : '120mins_decay',
    'd_3.85082e-05' : '300mins_decay',
    'd_1e-05' : '1200mins_decay' }<jupyter_output><empty_output><jupyter_text>Origins  Weights Define your gravity weights<jupyter_code>grav_models_Orig_wtless_dct = {
    'All_cities' : dict.fromkeys(['Pop_all_2011','Pop_all_2021','Unweighted']),
    'Dhaka_Chitt' : dict.fromkeys(['Pop_all_2011','Pop_all_2021','Unweighted']),
    'Minor_cities' : dict.fromkeys(['Pop_all_2011','Pop_all_2021','Unweighted']),
    'Dry_ports' : dict.fromkeys(['Unweighted']),
    'River_ports' : dict.fromkeys(['EXIM','Unweighted']),
    'Deep_sea_ports' : dict.fromkeys(['Annual_capacity_TEU','Annual_capacity_TEU_2025','Unweighted']),
    'All_SEZs' : dict.fromkeys(['All_SEZ_count']),
    'Active_SEZs' : dict.fromkeys(['Active_SEZ_count'])
}<jupyter_output><empty_output><jupyter_text>Gravity calculations Create a function to generate the gravity models and append them to a dictionary<jupyter_code>def orig_grav_loop(orig_file,od_matrix,orig_id_col,orig_weight_col,dest_weight_col):
    
    orig_gravity = calcOD.calculate_gravity(od = np.array(od_matrix), \
                                            oWeight = orig_weight_col,\
                                            dWeight = dest_weight_col)

    orig_gravity['O_ID'] = np.array(orig_id_col) # will assign incorrectly if not converted to an array first
    orig_gravity = orig_gravity.set_index('O_ID')
    
    # merge on indices to speed this step up
    orig_gravity_gdf = pd.merge(orig_file, orig_gravity, how='left',left_index=True,right_index=True)

    return orig_gravity_gdf<jupyter_output><empty_output><jupyter_text>Use the function to create gravity models for all the desired weights and export them to a geopackage<jupyter_code>from pprint import pprint
pprint(grav_models_Orig_wtless_dct)
for dest_key, weight_dct in grav_models_Orig_wtless_dct.items():
    
    for dest, input_list in odm_dct.items():
        
        # match on destination names so as to call the right weights
        
        if dest_key == dest:
            
            print(dest_key)
            
            for wt in weight_dct:

                print(wt)

                # Start timer
                func_start = time.time()

                # we call in the original destinations file here, referencing the weights column calculated above
                df = orig_grav_loop(input_list[0],input_list[1],orig_wtd.index,[1],input_list[2][wt]) 

                # create geodataframe, rename decay dict
                df = gpd.GeoDataFrame(df,geometry='geometry')
                df.rename(columns=rename_decay_dct,inplace=True)

                # export to geopackage
                gpkg_layername = str(dest_key) + '_' + wt # rename as needed
                df.to_file(os.path.join(f'results\\{prod_date}\\spatial\\origins_gravity_{prod_date}_{scenario}_{origins}_no_orig_weight.gpkg'),\
                                           layer=gpkg_layername,\
                                           driver="GPKG",\
                                           ignore_errors=True)

                # populate nested weight dict with gravity model
                grav_models_Orig_wtless_dct[dest_key][wt] = df

                # Report time

                func_end = time.time()
                print('\n time elapsed for function')
                print(str((func_end - func_start) / 60) + ' minutes \n')<jupyter_output>All_cities
Pop_all_2011

 time elapsed for function
22.75905690987905 minutes 

Pop_all_2021

 time elapsed for function
23.530368304252626 minutes 

Unweighted

 time elapsed for function
23.82742187579473 minutes 

Dhaka_Chitt
Pop_all_2011

 time elapsed for function
23.817124489943186 minutes 

Pop_all_2021

 time elapsed for function
24.154891757170358 minutes 

Unweighted

 time elapsed for function
23.197356200218202 minutes 

Minor_cities
Pop_all_2011

 time elapsed for function
22.830280526479086 minutes 

Pop_all_2021

 time elapsed for function
22.74421218633652 minutes 

Unweighted

 time elapsed for function
24.439106325308483 minutes 

Dry_ports
Unweighted

 time elapsed for function
24.71010059515635 minutes 

River_ports
EXIM

 time elapsed for function
23.31399001677831 minutes 

Unweighted

 time elapsed for function
24.555287166436514 minutes 

Deep_sea_ports
Annual_capacity_TEU

 time elapsed for function
24.80245396296183 minutes 

Annual_capacity_TEU_2025

 time el[...]<jupyter_text>Currently the gravity statistics are per origin point. Here we aggregate these statistics at the adm4 level, weighting by origin population<jupyter_code># grav_models_Orig_wtless_dct
# iterate through the resulting dict of gravity models, joining their results to adm3 and adm4 files for export and visualization

adm3_grav = adm3.set_index('adm3_pcode')
adm4_grav = adm4.set_index('adm4_pcode')

# adm3_grav.set_index('adm3_pcode',inplace=True)
# adm4_grav.set_index('adm4_pcode',inplace=True)

for key, weight_dct in grav_models_Orig_wtless_dct.items():
    
    print(key)
    
    for wt, layer in weight_dct.items():

        print(wt)

        # Start timer
        func_start = time.time()

        # Create separate adm3 and adm4 layers

        adm3_layer = layer.filter(regex='mins_decay|adm3_pcode|VALUE|adm3_pop')
        adm4_layer = layer.filter(regex='mins_decay|adm4_pcode|VALUE|adm4_pop')

        # weight index by population
        # Note from previous troubleshooting -- perhaps this WITH the OWeight is causing weird patterns

        for col in adm3_layer.loc[:,'2.9mins_decay':'1200mins_decay']:
            adm3_layer[col] = (adm3_layer[col] * (adm3_layer['VALUE'] / adm3_layer['adm3_pop'])) # weighting by share of overall adm population

        for col in adm4_layer.loc[:,'2.9mins_decay':'1200mins_decay']:
            adm4_layer[col] = (adm4_layer[col] * (adm4_layer['VALUE'] / adm4_layer['adm4_pop'])) # weighting by share of overall adm population

        # aggregate weighted indices by upazila

        upz = adm3_layer.filter(regex='mins_decay|adm3_pcode').groupby('adm3_pcode').apply(lambda x: x.sum()) # new version works with non-sequential columns
        upz.drop(columns='adm3_pcode',inplace=True)
        upz = upz.add_prefix(key + '_' + wt + '_') # if running scenarios, modify by scenario here

        # aggregate weighted indices by union

        union = adm4_layer.filter(regex='mins_decay|adm4_pcode').groupby('adm4_pcode').apply(lambda x: x.sum()) # new version works with non-sequential columns
        union.drop(columns='adm4_pcode',inplace=True)
        union = union.add_prefix(key + '_' + wt + '_')

        # join to adm4 file

        adm3_grav = pd.merge(adm3_grav,upz,how='left',left_index=True,right_index=True)
        adm4_grav = pd.merge(adm4_grav,union,how='left',left_index=True,right_index=True) 
        
        # Report time

        func_end = time.time()
        print('\n time elapsed for function')
        print(str((func_end - func_start) / 60) + ' minutes \n')<jupyter_output>All_cities
Pop_all_2011<jupyter_text>Export the finished adm3 and adm4 layers for data analysis and visualization<jupyter_code># adm3

adm3_grav.drop('geometry',axis=1).to_csv(os.path.join(res_pth,prod_date,f'tables\\adm3_gravity_models_{prod_date}_{origins}_nopop_wt.csv'))
adm3_grav.to_file(os.path.join(res_pth,prod_date,f'spatial\\adm3_gravity_models_{scenario}_{prod_date}_{origins}_nopop_wt.geojson'),driver="GeoJSON")
adm3_grav.to_file(os.path.join(res_pth,prod_date,f'spatial\\adm3_gravity_models_{scenario}_{prod_date}_{origins}_nopop_wt.gpkg'),driver="GPKG", ignore_errors=True)

# adm4
adm4_grav.drop('geometry',axis=1).to_csv(os.path.join(res_pth,prod_date,f'tables\\adm4_gravity_models_{scenario}_{prod_date}_{origins}_nopop_wt.csv'))
adm4_grav.to_file(os.path.join(res_pth,prod_date,f'spatial\\adm4_gravity_models_{scenario}_{prod_date}_{origins}_nopop_wt.geojson'),driver="GeoJSON")
adm4_grav.to_file(os.path.join(res_pth,prod_date,f'spatial\\adm4_gravity_models_{scenario}_{prod_date}_{origins}_nopop_wt.gpkg'),driver="GPKG", ignore_errors=True)<jupyter_output><empty_output><jupyter_text>Destinations<jupyter_code># {scenario}_mkts_odm.transpose()
# .iloc[1:]
# List previously created weights

gm_weights = ['good_jobs_tot_','bad_jobs_tot_','good_jobs_4and5','bad_jobs_4and5',\
              'n_firms_4and5','n_firms_tot','tw_tot','tw_tot_4and5','tw_ser_4and5',\
              'tw_ind_4and5','no_weight']

# Insert weights into dict

grav_models_dict = dict.fromkeys(gm_weights)<jupyter_output><empty_output><jupyter_text>Create a function to generate the gravity models per destination<jupyter_code>def dest_grav_loop(orig_file,od_matrix,grav_param,dest_df,orig_weight_df):
    dest_gravity = calcOD.calculate_gravity(np.array(od_matrix), oWeight=dest_df[grav_param], dWeight = orig_weight_df['Pop'])
    dest_gravity['NN'] = dest_df['NN']
    dest_gravity_gdf = pd.merge(mkts, dest_gravity, how='left',on='NN')

    return dest_gravity_gdf

# OLD

# def dest_grav_loop(orig_file,od_matrix,grav_param,oW_df,dW):
#     dest_gravity = calcOD.calculate_gravity(np.array(od_matrix), oWeight=oW_df[grav_param], dWeight = dW['Pop'])
#     dest_gravity['NN'] = oW_df['NN']
#     dest_gravity_gdf = pd.merge(mkts, dest_gravity, how='left',on='NN')

#     return dest_gravity_gdf<jupyter_output><empty_output><jupyter_text>Apply the function and append results to a dictionary<jupyter_code>for key, v in grav_models_dict.items():
    
    print(key)
    
    df = dest_grav_loop(current_orig_mkts,current_mkts_odm.transpose(),key,mkts,cur_orig_grp)
    
    gpkg_layername = 'mkts_cur_' + key
    df.to_file(os.path.join('results\\spatial\\gravity_models_mkts.gpkg'),layer=gpkg_layername,driver="GPKG", ignore_errors=True)
    
    grav_models_dict[key] = df<jupyter_output><empty_output><jupyter_text>OLD code Sorting process<jupyter_code># # OLD manual process

# # sort by O_IDs, then dest node IDs

# current_mkts_odm = current_mkts_odm.sort_values('Unnamed: 0').reindex(sorted_nicely(current_mkts_odm.columns), axis=1)
# ua_mkts_odm = ua_mkts_odm.sort_values('Unnamed: 0').reindex(sorted_nicely(ua_mkts_odm.columns), axis=1)

# # reset O_ID column to the front
# current_mkts_odm = current_mkts_odm[ ['Unnamed: 0'] + [ col for col in current_mkts_odm.columns if col != 'Unnamed: 0' ] ]
# ua_mkts_odm = ua_mkts_odm[ ['Unnamed: 0'] + [ col for col in ua_mkts_odm.columns if col != 'Unnamed: 0' ] ]

# # set the Dest_ID column back to index so the shape is the same as the dWeight shape
# current_mkts_odm.set_index('Unnamed: 0',inplace=True)
# ua_mkts_odm.set_index('Unnamed: 0',inplace=True)<jupyter_output><empty_output><jupyter_text>Georeference, spatial join files<jupyter_code># # Old manual routine 

# current_orig_mkts['geometry'] = current_orig_mkts['geometry'].apply(wkt.loads)
# current_orig_mkts = gpd.GeoDataFrame(current_orig_mkts,geometry='geometry')
# current_orig_mkts = gpd.sjoin(current_orig_mkts,adm4[['geometry','adm3_en','adm4_en','adm3_pcode','adm4_pcode']],how='left',op='within')

# ua_orig_mkts['geometry'] = ua_orig_mkts['geometry'].apply(wkt.loads)
# ua_orig_mkts = gpd.GeoDataFrame(ua_orig_mkts,geometry='geometry')
# ua_orig_mkts = gpd.sjoin(ua_orig_mkts,adm4[['geometry','adm3_en','adm4_en','adm3_pcode','adm4_pcode']],how='left',op='within')<jupyter_output>C:\Users\wb467985\.conda\envs\geo\lib\site-packages\geopandas\tools\sjoin.py:61: UserWarning: CRS of frames being joined does not match!(None != epsg:4326)
  "(%s != %s)" % (left_df.crs, right_df.crs)<jupyter_text>Scenario based approach to gravity stats Now repeat the process for the all updated roads<jupyter_code># for key, v in grav_models_dict.items():
    
#     print(key)
    
#     df = orig_grav_loop(ua_orig_mkts,ua_mkts_odm,ua_orig_grp,mkts[key])
    
#     gpkg_layername = 'mkts_ua_' + key
#     df.to_file(os.path.join('results\\spatial\\gravity_models.gpkg'),layer=gpkg_layername,driver="GPKG", ignore_errors=True)
    
#     grav_models_dict[key] = df
# # iterate through the dict

# # adm4_grav = adm4_grav

# for key, layer in grav_models_dict.items():
    
#     print(key)
    
#     # weight index by population
    
#     layer['adm4_pop'] = layer.groupby('adm4_pcode')['VALUE'].transform(np.sum)
    
#     for col in layer.loc[:,'d_0.01':'d_1e-05']:
#         layer[col] = (layer[col] * (layer['VALUE'] / layer['adm4_pop']))
    
# #     [ col * (layer['VALUE'].astype(float) / layer['adm4_pop'].astype(float)) for col in layer.loc[:,'d_0.01':'d_1e-05']] # failed atempt at list comprehension
    
#     # aggregate weighted indices by union
    
#     union = layer.loc[:,'adm4_pcode':'d_1e-05'].groupby('adm4_pcode').apply(lambda x: x.sum())
#     union.drop(columns='adm4_pcode',inplace=True)
#     union = union.add_prefix(key + '_ua_')
#     union.reset_index(inplace=True)
#     union['adm4_pcode'] = union['adm4_pcode'].astype(str)

#     # join to adm4 file

#     adm4_grav = pd.merge(adm4_grav,union,how='left',on=['adm4_pcode'])<jupyter_output><empty_output><jupyter_text>If doing a scenario, calculate changes in gravity resulting from that scenario<jupyter_code># adm4_grav['delta60_gj_tot'] = adm4_grav['good_jobs_tot__ua_d_0.0001925409'] - adm4_grav['good_jobs_tot__cur_d_0.0001925409']
# adm4_grav['delta60_bj_tot'] = adm4_grav['bad_jobs_tot__ua_d_0.0001925409'] - adm4_grav['bad_jobs_tot__cur_d_0.0001925409']
# adm4_grav['delta30_bj_tot'] = adm4_grav['bad_jobs_tot__ua_d_0.0003850818'] - adm4_grav['bad_jobs_tot__cur_d_0.0003850818']

# adm4_grav['delta60_gj_4and5'] = adm4_grav['good_jobs_4and5_ua_d_0.0001925409'] - adm4_grav['good_jobs_4and5_cur_d_0.0001925409']
# adm4_grav['delta60_bj_4and5'] = adm4_grav['bad_jobs_4and5_ua_d_0.0001925409'] - adm4_grav['bad_jobs_4and5_cur_d_0.0001925409']

# adm4_grav['delta60_tw_tot'] = adm4_grav['tw_tot_ua_d_0.0001925409'] - adm4_grav['tw_tot_cur_d_0.0001925409']
# adm4_grav['delta60_nf_tot'] = adm4_grav['n_firms_tot_ua_d_0.0001925409'] - adm4_grav['n_firms_tot_cur_d_0.0001925409']

# adm4_grav['delta60_tw_4and5'] = adm4_grav['tw_tot_4and5_ua_d_0.0001925409'] - adm4_grav['tw_tot_4and5_cur_d_0.0001925409']
# adm4_grav['delta60_nf_4and5'] = adm4_grav['n_firms_4and5_ua_d_0.0001925409'] - adm4_grav['n_firms_4and5_cur_d_0.0001925409']

# adm4_grav['delta60_tw_ser4and5'] = adm4_grav['tw_ser_4and5_ua_d_0.0001925409'] - adm4_grav['tw_ser_4and5_cur_d_0.0001925409']
# adm4_grav['delta60_tw_ind4and5'] = adm4_grav['tw_ind_4and5_ua_d_0.0001925409'] - adm4_grav['tw_ind_4and5_cur_d_0.0001925409']

# adm4_grav['delta60_now'] = adm4_grav['no_weight_ua_d_0.0001925409'] - adm4_grav['no_weight_cur_d_0.0001925409']
# adm4_grav['delta30_now'] = adm4_grav['no_weight_ua_d_0.0003850818'] - adm4_grav['no_weight_cur_d_0.0003850818']<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1192435172548255, 'fraction_numerical': 0.0396568531877559, 'mean_word_length': 5.027262044653349, 'pattern_counts': {'":': 0, '<': 98, '<?xml version=': 0, '>': 100, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '126682', 'license': 'MIT', 'n_tokens_mistral': 10208, 'n_tokens_neox': 9330, 'n_words': 1567, 'path': 'Implementations/FY21/ACC_BGD_National/Step 5 - Gravity Modeling.ipynb', 'repo_name': 'worldbank/GOST_PublicGoods', 'size': '25645'}