<jupyter_start><jupyter_text>Proposed Roads Processing and Raster Export This notebook processes a collection of line vectors into raster files with the 'speed' attribute used as the raster pixel value.<jupyter_code>import os, sys
from datetime import date

import common_rasterio_ops as rast_ops

import re
import numpy as np
from numpy import random

import rasterio
from rasterio import features, transform
from rasterio.transform import Affine
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.io import MemoryFile

import pandas as pd
import geopandas as gpd

import shapely
from shapely.geometry import shape, Point, LineString, MultiLineString

sys.path.append('../../src/')

from gostrocks.src.GOSTRocks.misc import tPrint<jupyter_output><empty_output><jupyter_text>Setup File paths<jupyter_code>geo_pth = r'P:\PAK\GEO'
data_dir = r'../../data'

rast_dir = r'rast_inputs'
vect_in_dir = r'vect_inputs'
vect_out_dir = r'vect_out'

rd_dir = r'roads'
dest_dir = r'destinations'
speed_dir = r'speed'
fric_dir = r'friction'
acc_dir = r'access'<jupyter_output><empty_output><jupyter_text>Date<jupyter_code>today = date.today().strftime("%y%m%d")<jupyter_output><empty_output><jupyter_text>Projections<jupyter_code># change this to whatever the desired metric output projection is -- this notebook assumes you want outputs in meters
dest_crs = 'EPSG:32642'<jupyter_output><empty_output><jupyter_text>Load in KP as clipping object<jupyter_code>kp = gpd.read_file(os.path.join(geo_pth,r'Boundaries/OCHA/pak_admbnda_adm1_ocha_pco_gaul_20181218.shp'))
kp = kp[kp['ADM1_EN'] == 'Khyber Pakhtunkhwa']
kp = kp.to_crs(dest_crs)

# Buffer the polygon by 20km so we take in nearby markets and roads that may be used
kp.geometry = kp.buffer(20000)<jupyter_output><empty_output><jupyter_text>Important file names<jupyter_code>prop_tab_raw = r'Proposed_raw//needed_proposed_roads_data_211017.csv'
pgr = r'Proposed_raw//KP_RoadforUpgrade_Correct_211014'
prop_geo_raw = pgr +'.gpkg'
prop_geo_clean = pgr + '_clean.shp'<jupyter_output><empty_output><jupyter_text>Pre-processing Tabular<jupyter_code>prop_tab = pd.read_csv(os.path.join(data_dir,rd_dir,prop_tab_raw),na_values='')
prop_tab[::20]
na_fill_dict = {
    'Current_Road_Cond' : 'Poor',
    'Current_Surface' : 'Unpaved',
    'Current_Road_Class' : 'Access Road',
    'Upgrade_Road_Cond' : 'Very Good',
    'Upgrade_Surface' : 'Asphaltic',
    'Upgrade_Road_Class' : 'District Road'
}
prop_tab = prop_tab.fillna(na_fill_dict)
prop_tab.to_csv(os.path.join(data_dir,rd_dir,r'Proposed_raw//needed_proposed_roads_data_211004_nafill.csv'))
prop_tab.columns
prop_tab = prop_tab[['Road_ID','Road_Name','District','Current_Road_Cond','Current_Surface','Current_Road_Class','Upgrade_Road_Cond','Upgrade_Surface','Upgrade_Road_Class']]
prop_tab.rename({'Road_ID' : 'SN'},axis=1,inplace=True)
# prop_tab.to_csv(os.path.join(data_dir,rd_dir,r'Proposed_raw//Clean_proposed_roads_attributes_211017.csv'))<jupyter_output><empty_output><jupyter_text>Geo<jupyter_code>pre_prop = gpd.read_file(os.path.join(data_dir,rd_dir,prop_geo_raw),driver="GPKG").to_crs(dest_crs)
pre_prop = pre_prop[pre_prop['SN'].isna() == False]
pre_prop = pre_prop[['SN','geometry']]
# pre_prop.to_file(os.path.join(data_dir,rd_dir,prop_geo_clean))<jupyter_output><empty_output><jupyter_text>**Stop**At this stage go to ArcGIS Pro and "Add Surface Information" from a DEM to the cleaned roads. Add Z_Mean, S_Length, and Avg_Slope to the road lines.Hacky yes, but this is very painful to do in Python and very simple to do there. Load the data back in below. Once you are done continue on in this notebook.  Load in and process roads vectors Note: need to add elevation data and slope to proposed roads prior to this step using ArcGIS's Add Surface Info<jupyter_code>prop_tab = pd.read_csv(os.path.join(data_dir,rd_dir,r'Proposed_raw//Clean_proposed_roads_attributes_211017.csv'))
# load back in your roads data after joining surface information in ArcGIS
prop_raw = gpd.read_file(os.path.join(data_dir,rd_dir,prop_geo_clean)).to_crs(dest_crs)
prop_raw.head(2)<jupyter_output><empty_output><jupyter_text>Each road has multiple segments which are of different lengths and slopes. To get the most accurate value for each road we must weight segment's values by their length when combining Calculate length-weighted slope and elevation by road<jupyter_code># calculate road segment weight (based on length as part of road's total length

rds_length = prop_raw.groupby('SN')['SLength'].sum().rename('total_length',axis=1)
prop_raw = pd.merge(prop_raw,rds_length,how='left',on='SN')
prop_raw['length_wt'] = prop_raw['SLength'] / prop_raw['total_length']
# calculate weighted elev and slope per road segment

prop_raw['Z_Mean'] = prop_raw.Z_Mean * prop_raw.length_wt
prop_raw['Avg_Slope'] = prop_raw.Avg_Slope * prop_raw.length_wt
# calculate weighted average elevation and slope per road, slim down to one data frame, join in other tabular info

prop_slim = prop_raw.groupby('SN')[['Z_Mean','Avg_Slope']].sum()
prop_slim = pd.merge(prop_slim,prop_tab,how='left',on='SN')
prop_slim.SN = prop_slim.SN.astype(int)<jupyter_output><empty_output><jupyter_text>Dissolve the results by SN<jupyter_code># the agg function just works on the road lengths
prop_d = prop_raw.dissolve(by='SN',aggfunc='sum').reset_index()<jupyter_output><empty_output><jupyter_text>Consolidate the road lines on their common IDs and then join the resulting geometries to the slimmed down proposed roads DF<jupyter_code>prop_raw = gpd.GeoDataFrame(pd.merge(prop_slim,prop_d[['SN','geometry']],how='left',on='SN'))
prop_raw[::15]<jupyter_output><empty_output><jupyter_text>Populate speeds for the roads based on their type, status, pavement, and locationFirst align the data model with the existing roads data model<jupyter_code># create a temporary copy to iterate over
prop_rds = prop_raw.copy()
prop_rds.dtypes<jupyter_output><empty_output><jupyter_text>Terrain<jupyter_code># narrowing range to 8/16 based on satellite imagery observation
prop_rds['Terrain'] = pd.cut(prop_rds['Avg_Slope'], [-np.inf, 8, 16, np.inf], 
                           labels = ['Plains', 'Hills', 'Mountains']) # change labels here<jupyter_output><empty_output><jupyter_text>Conditions -- leftover<jupyter_code>prop_rds.Current_Road_Class = prop_rds.Current_Road_Class.replace('Village road','Collector Road')
prop_rds.Current_Road_Class = prop_rds.Current_Road_Class.fillna('Collector Road')
road_condition_filter = [prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('Collector'),
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('Collector'),
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('Collector')]
prop_rds.columns<jupyter_output><empty_output><jupyter_text>Calculate Speeds Base speeds from road type and terrain<jupyter_code>terrain_current_class_filter = [prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Current_Road_Class'].str.contains('Collector'),
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Current_Road_Class'].str.contains('Collector'),
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Current_Road_Class'].str.contains('Collector')]
                        
terrain_upgrade_class_filter = [prop_rds['Terrain'].str.contains('Plains') & prop_rds['Upgrade_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Upgrade_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Upgrade_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Plains') & prop_rds['Upgrade_Road_Class'].str.contains('Collector'),
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Upgrade_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Upgrade_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Upgrade_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Hills') & prop_rds['Upgrade_Road_Class'].str.contains('Collector'),
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Upgrade_Road_Class'].str.contains('Provincial'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Upgrade_Road_Class'].str.contains('District'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Upgrade_Road_Class'].str.contains('Access'),    
    prop_rds['Terrain'].str.contains('Mountains') & prop_rds['Upgrade_Road_Class'].str.contains('Collector')]
                        
# Corresponding list of speeds, with line breaks by terrain type for editing and readability
speeds_lst = [80,50,40,20,\
              60,40,30,15,\
              40,30,20,10]
# initial dry speeds to calculate msn speeds
prop_rds['current_base_speed'] = np.select(terrain_current_class_filter,speeds_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road
prop_rds['upgrade_base_speed'] = np.select(terrain_upgrade_class_filter,speeds_lst,default=0.5)<jupyter_output><empty_output><jupyter_text>**Modify based on condition + season** Define a filter based on surface + condition that we'll use to assign speed modifiers in each season. Note that we add an `Asphaltic` surface type unique to proposed roads because this was specified by the client. We assume it's very good and has minimal speed reductions even under difficult conditions<jupyter_code>current_class_filter = [prop_rds['Current_Surface'].str.contains('Earthen') & prop_rds['Current_Road_Cond'].str.contains('Good'),    
                        prop_rds['Current_Surface'].str.contains('Earthen') & prop_rds['Current_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Current_Surface'].str.contains('Earthen') & prop_rds['Current_Road_Cond'].str.contains('Poor'),
                        prop_rds['Current_Surface'].str.contains('Gravel') & prop_rds['Current_Road_Cond'].str.contains('Good'),    
                        prop_rds['Current_Surface'].str.contains('Gravel') & prop_rds['Current_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Current_Surface'].str.contains('Gravel') & prop_rds['Current_Road_Cond'].str.contains('Poor'),
                        prop_rds['Current_Surface'].str.contains('Paved') & prop_rds['Current_Road_Cond'].str.contains('Good'),    
                        prop_rds['Current_Surface'].str.contains('Paved') & prop_rds['Current_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Current_Surface'].str.contains('Paved') & prop_rds['Current_Road_Cond'].str.contains('Poor'),
                        prop_rds['Current_Surface'].str.contains('Asphaltic') & prop_rds['Current_Road_Cond'].str.contains('Good'),    
                        prop_rds['Current_Surface'].str.contains('Asphaltic') & prop_rds['Current_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Current_Surface'].str.contains('Asphaltic') & prop_rds['Current_Road_Cond'].str.contains('Poor')]

upgrade_class_filter = [prop_rds['Upgrade_Surface'].str.contains('Earthen') & prop_rds['Upgrade_Road_Cond'].str.contains('Good'),    
                        prop_rds['Upgrade_Surface'].str.contains('Earthen') & prop_rds['Upgrade_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Upgrade_Surface'].str.contains('Earthen') & prop_rds['Upgrade_Road_Cond'].str.contains('Poor'),
                        prop_rds['Upgrade_Surface'].str.contains('Gravel') & prop_rds['Upgrade_Road_Cond'].str.contains('Good'),    
                        prop_rds['Upgrade_Surface'].str.contains('Gravel') & prop_rds['Upgrade_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Upgrade_Surface'].str.contains('Gravel') & prop_rds['Upgrade_Road_Cond'].str.contains('Poor'),
                        prop_rds['Upgrade_Surface'].str.contains('Paved') & prop_rds['Upgrade_Road_Cond'].str.contains('Good'),    
                        prop_rds['Upgrade_Surface'].str.contains('Paved') & prop_rds['Upgrade_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Upgrade_Surface'].str.contains('Paved') & prop_rds['Upgrade_Road_Cond'].str.contains('Poor'),
                        prop_rds['Upgrade_Surface'].str.contains('Asphaltic') & prop_rds['Upgrade_Road_Cond'].str.contains('Good'),    
                        prop_rds['Upgrade_Surface'].str.contains('Asphaltic') & prop_rds['Upgrade_Road_Cond'].str.contains('Fair'),    
                        prop_rds['Upgrade_Surface'].str.contains('Asphaltic') & prop_rds['Upgrade_Road_Cond'].str.contains('Poor')]<jupyter_output><empty_output><jupyter_text>Dry season<jupyter_code># Corresponding list of dry season speed modifiers, with line breaks by terrain type for editing and readability

dry_mods_lst = [ 0.6, 0.5, 0.4,\
                 0.75, 0.6, 0.5,\
                 1, 0.8, 0.6, 
                 1, 1, 0.8]

prop_rds['current_dry_mod'] = np.select(current_class_filter,dry_mods_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road
prop_rds['upgrade_dry_mod'] = np.select(upgrade_class_filter,dry_mods_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road

prop_rds['current_dry_speed'] = prop_rds.current_dry_mod * prop_rds.current_base_speed
prop_rds['upgrade_dry_speed'] = prop_rds.upgrade_dry_mod * prop_rds.upgrade_base_speed
# check that the last values are not null (indicating all speeds were successfully calculated)
prop_rds.sort_values('current_dry_speed').current_dry_speed.tail()
# check that the last values are not null (indicating all speeds were successfully calculated)
prop_rds.sort_values('upgrade_dry_speed').upgrade_dry_speed.tail()<jupyter_output><empty_output><jupyter_text>Monsoon<jupyter_code># Corresponding list of monsoon speed modifiers, with line breaks by terrain type for editing and readability

msn_mods_lst = [ 0.4, 0.3, 0.2,\
                 0.7, 0.5, 0.4,\
                 0.9, 0.75, 0.5,
                 1, 0.9, 0.75]

prop_rds['current_msn_mod'] = np.select(current_class_filter,msn_mods_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road
prop_rds['upgrade_msn_mod'] = np.select(upgrade_class_filter,msn_mods_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road

prop_rds['current_msn_speed'] = prop_rds.current_msn_mod * prop_rds.current_base_speed
prop_rds['upgrade_msn_speed'] = prop_rds.upgrade_msn_mod * prop_rds.upgrade_base_speed<jupyter_output><empty_output><jupyter_text>Winter<jupyter_code># Corresponding list of monsoon speed modifiers, with line breaks by terrain type for editing and readability

winter_mods_lst = [ 0.4, 0.3, 0.2,\
                    0.7, 0.5, 0.4,\
                    1, 0.75, 0.5,
                    1, 0.9, 0.75]

prop_rds['current_winter_mod'] = np.select(current_class_filter,winter_mods_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road
prop_rds['upgrade_winter_mod'] = np.select(upgrade_class_filter,winter_mods_lst,default=0.5) # very low default as path speeds will be re-calculated separately using terrain down the road
# Define which winter mods to apply to which road based on their terrain

winter_cutoff_lst = [prop_rds['Terrain'] == 'Plains', \
                     prop_rds['Terrain'] == 'Hills', \
                     prop_rds['Terrain'] == 'Mountains']
winter_current_mod_revised =  [ prop_rds['current_dry_speed'], \
                       (prop_rds['current_dry_speed'] * (prop_rds['current_winter_mod'] + (1 - prop_rds['current_winter_mod']) /2) ), \
                       (prop_rds['current_base_speed'] * prop_rds['current_winter_mod'])]

winter_upgrade_mod_revised =  [ prop_rds['upgrade_dry_speed'], \
                       (prop_rds['upgrade_dry_speed'] * (prop_rds['upgrade_winter_mod'] + (1 - prop_rds['upgrade_winter_mod']) /2) ), \
                       (prop_rds['upgrade_base_speed'] * prop_rds['upgrade_winter_mod'])]
prop_rds['current_winter_speed'] = np.select(winter_cutoff_lst, winter_current_mod_revised, default=prop_rds['current_dry_speed'])
prop_rds['upgrade_winter_speed'] = np.select(winter_cutoff_lst, winter_upgrade_mod_revised, default=prop_rds['upgrade_dry_speed'])
# check that the last values are not null (indicating all speeds were successfully calculated)
prop_rds.sort_values('current_winter_speed').current_winter_speed.tail()
# check that the last values are not null (indicating all speeds were successfully calculated)
prop_rds.sort_values('upgrade_dry_speed').upgrade_dry_speed.tail()<jupyter_output><empty_output><jupyter_text>Check<jupyter_code>test = prop_rds[['SN','Current_Road_Cond','Current_Surface','Current_Road_Class',\
              'Upgrade_Road_Cond','Upgrade_Surface','Upgrade_Road_Class',\
              'Terrain',
              'current_dry_mod','current_msn_mod','current_winter_mod',
              'current_base_speed','current_dry_speed','current_msn_speed','current_winter_speed',\
              'upgrade_base_speed','upgrade_dry_speed','upgrade_msn_speed','upgrade_winter_speed']]
# if all results are 0 or positive this means dry speeds are >= than winter speeds, as they should be
test['spd_dif'] = test['current_dry_speed'] - test['current_winter_speed']<jupyter_output>C:\Users\wb467985\AppData\Local\Temp/ipykernel_17920/1935561870.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  test['spd_dif'] = test['current_dry_speed'] - test['current_winter_speed']<jupyter_text>Calculate the distance to the Nth nearest destination of each type from the end point of the line This is necessary to clip the later friction and access surface outputs -- otherwise they will be enormous and difficult to work with<jupyter_code># adapted from: https://automating-gis-processes.github.io/site/notebooks/L3/nearest-neighbor-faster.html
# more info here https://towardsdatascience.com/using-scikit-learns-binary-trees-to-efficiently-find-latitude-and-longitude-neighbors-909979bd929b

from sklearn.neighbors import BallTree

def get_nth_nearest(src_points, candidates, k_neighbors=5):
    """Find nearest K neighbors for all source points from a set of candidate points, return Kth neighbor. If smaller than K, return nearest."""

    # Create tree from the candidate points
    tree = BallTree(candidates, leaf_size=15, metric='haversine')

    # Find closest points and distances
    
    if len(candidates) < k_neighbors:
        distances, indices = tree.query(src_points, k=1)
    else:
        distances, indices = tree.query(src_points, k=k_neighbors)

    # Transpose to get distances and indices into arrays
    distances = distances.transpose()
    indices = indices.transpose()

    # Get closest indices and distances (i.e. array at index 0)
    # note: for the second closest points, you would take index 1, etc.
    # Robert -- setting to the 5th nearest feature as an arbitrary cutoff
    
    if len(candidates) < k_neighbors:
        closest = indices[0]
        closest_dist = distances[0]
    else:
        closest = indices[(k_neighbors - 1)]
        closest_dist = distances[(k_neighbors - 1)]

    # Return indices and distances
    return (closest, closest_dist)

def nearest_neighbors(left_gdf, right_gdf, return_dist=False, k_neighbors=5):
    """
    For each point in left_gdf, find closest point in right GeoDataFrame and return them.

    NOTICE: Assumes that the input Points are in WGS84 projection (lat/lon).
    """

    left_geom_col = left_gdf.geometry.name
    right_geom_col = right_gdf.geometry.name

    # Ensure that index in right gdf is formed of sequential numbers
    right = right_gdf.copy().reset_index(drop=True)

    # Parse coordinates from points and insert them into a numpy array as RADIANS
    left_radians = np.array(left_gdf[left_geom_col].apply(lambda geom: (geom.y * np.pi / 180, geom.x * np.pi / 180)).to_list())
    right_radians = np.array(right[right_geom_col].apply(lambda geom: (geom.y * np.pi / 180, geom.x * np.pi / 180)).to_list())

    # Find the nearest points
    # -----------------------
    # closest ==> index in right_gdf that corresponds to the closest point
    # dist ==> distance between the nearest neighbors (in meters)

    closest, dist = get_nth_nearest(src_points=left_radians, candidates=right_radians, k_neighbors=k_neighbors)

    # Return points from right GeoDataFrame that are closest to points in left GeoDataFrame
    closest_points = right.loc[closest]

    # Ensure that the index corresponds the one in left_gdf
    closest_points = closest_points.reset_index(drop=True)

    # Add distance if requested
    if return_dist:
        # Convert to meters from radians
        earth_radius = 6371000  # meters
        closest_points['distance'] = dist * earth_radius

    return closest_points

# from https://gis.stackexchange.com/questions/378806/multi-part-geometries-do-not-provide-a-coordinate-sequence-error-when-extracti
# only works on LineStrings, not MultiLineStrings

def add_start_end_nodes_to_gdf(gdf):
    
# This function adds start and end nodes to geodataframe

    gdf['start_node'] = None
    gdf['end_node'] = None

    for index, row in gdf.iterrows():
        coords = [(coords) for coords in list(row['geometry'].coords)]
        start_node, end_node = [coords[i] for i in (0, -1)]
        gdf.at[index, 'start_node'] = start_node
        gdf.at[index, 'end_node'] = end_node<jupyter_output><empty_output><jupyter_text>Set up a dict of destinations we can loop over<jupyter_code>dest_fils = {
    re.findall(r'KP_(.*?).gpkg',fil)[0]: gpd.clip(gpd.read_file(os.path.join(data_dir,dest_dir,'reference',fil))\
                                                      .reset_index(drop=True)\
                                                      .set_crs(4326).to_crs(dest_crs),\
                                                  kp) \
    for fil in os.listdir(os.path.join(data_dir,dest_dir,'reference')) if fil.endswith(".gpkg")
}
dest_fils.keys()<jupyter_output><empty_output><jupyter_text>Calculate distance from each line segment to the Nth destination point (or nearest if < N features)<jupyter_code># sort notebook

prop_rds = prop_rds.sort_values('SN').reset_index(drop=True)

# create a temporary df we'll use to perform calculations

prop_rds_pts = prop_rds.copy()

# create single LineString objects so we can calculate their endpoints and measure from those
# NOTE this is hacky and problematic if the multilinestrings are really long
# TO DO: find the center point of the Linestring instead

prop_rds_pts.geometry = pd.Series([i[0]  if i.type == 'MultiLineString' else i for i in prop_rds_pts.geometry])
prop_rds_pts = prop_rds_pts.set_geometry('geometry')

# calculate their start / end points

add_start_end_nodes_to_gdf(prop_rds_pts)

# make the end point the geometry

prop_rds_pts['pt_geometry'] = pd.Series([Point(i) for i in prop_rds_pts.end_node])
prop_rds_pts = prop_rds_pts.set_geometry('pt_geometry')

# calculate the distance to the Nth feature on the prop_rds GDF using this pts GDF. 
# N =5 for most features, 2 for very dispersed features like markets and HQs. You may wish to vary these parameters based on your context

for key, gdf in dest_fils.items():
    
    if key in ['District_HQs','markets_All','markets_Central']:
    # pad the distance calculated by the road's length and an extra half kilometer to ensure no part of the road is ever excluded
    
        prop_rds[f'{key}_nth_nearest_distance'] = (nearest_neighbors(prop_rds_pts.to_crs(4326),\
                                                                     gdf.to_crs(4326),\
                                                                     return_dist=True,
                                                                      k_neighbors=2)['distance'] + prop_rds_pts['geometry'].length + 500).astype(int)
    else: 
    # pad the distance calculated by the road's length and an extra half kilometer to ensure no part of the road is ever excluded
    
        prop_rds[f'{key}_nth_nearest_distance'] = (nearest_neighbors(prop_rds_pts.to_crs(4326),\
                                                                     gdf.to_crs(4326),\
                                                                     return_dist=True,
                                                                      k_neighbors=5)['distance'] + prop_rds_pts['geometry'].length + 500).astype(int)

# add start and end nodes in case they're useful later

prop_rds.start_node = prop_rds_pts.start_node.astype(str)
prop_rds.end_node = prop_rds_pts.end_node.astype(str)

prop_rds.columns
prop_rds[prop_rds['SN'] == 60]['health_secondary_nth_nearest_distance']<jupyter_output><empty_output><jupyter_text>Export<jupyter_code>prop_rds.Terrain = prop_rds.Terrain.astype(str)
# prop_rds.to_file(os.path.join(data_dir,rd_dir,r'Proposed_final//Proposed_roads_processed_211021_second_dist.gpkg'),driver="GPKG")
prop_rds.to_file(os.path.join(data_dir,rd_dir,f'Proposed_final//Proposed_roads_processed_{today}_test.gpkg'),driver="GPKG")<jupyter_output><empty_output><jupyter_text>Manually merge fifth and second nearest distance columns Hacky solution to the problem that the 5th nearest distance is way too far for a few features -- Distritc HQs and all types of markets. This pointlessly extends processing time. TO DO make this flexibility part of the core code above.<jupyter_code>pr_second = gpd.read_file(os.path.join(data_dir,rd_dir,r'Proposed_final//Proposed_roads_processed_211021_second_dist.gpkg'),driver="GPKG")
pr_fifth = gpd.read_file(os.path.join(data_dir,rd_dir,r'Proposed_final//Proposed_roads_processed_211014.gpkg'),driver="GPKG")
dif_cols = ['District_HQs_fifth_nearest_distance','markets_All_fifth_nearest_distance','markets_Central_fifth_nearest_distance']
# pr_second.columns
pr_fifth.drop(dif_cols,axis=1,inplace=True)
pr_final = pd.merge(pr_fifth,pr_second[['SN','District_HQs_fifth_nearest_distance','markets_All_fifth_nearest_distance','markets_Central_fifth_nearest_distance']],how='left',on='SN')
pr_final.to_file(os.path.join(data_dir,rd_dir,r'Proposed_final//Proposed_roads_processed_211021.gpkg'),driver="GPKG")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12200557103064066, 'fraction_numerical': 0.010933147632311978, 'mean_word_length': 4.618348982785602, 'pattern_counts': {'":': 0, '<': 107, '<?xml version=': 0, '>': 107, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '318291', 'license': 'CC0-1.0', 'n_tokens_mistral': 9899, 'n_tokens_neox': 9278, 'n_words': 1777, 'path': 'notebooks/Roadbuilding_Scenario_Evaluation/Step 1 - Proposed Roads Preparation.ipynb', 'repo_name': 'worldbank/Khyber-Pakhtunkhwa-Accessibility-Analysis', 'size': '28720'}