In [ ]:
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import ogr, gdal
import os
import sqlite3
import json
import requests
import shapely
import math
from shapely.geometry import shape 
In [ ]:
#config
data_folder = 'data/'
layertypes = ["multipolygons", "points"] #Type of features
sqlite_path = data_folder + 'sqlite.db' #location of SQLite file
csv_path = data_folder + 'csv/' #location of csv's on folder
shp_path = data_folder + 'shp/' #location of csv's on folder
geojson_path = data_folder + 'geojson/' #location of csv's on folder
pbf_path = data_folder + 'map.osm' #location of PBF file
osrm_server = 'http://192.168.0.5:5000' #osrm
filter_type = ['village', 'town'] #type of locality
degree_to_meter = 100977 #degree to Meter for romania
buffer_point = 500 #Point to polygon
buffer_point = buffer_point/degree_to_meter #Meter to degree for romania
In [ ]:
#Define Input datasource
driver = ogr.GetDriverByName('OSM') # Format
data_source = ogr.Open(pbf_path, 0) #open Data
LayerCount = data_source.GetLayerCount()

layers = []
# get Required Layer
for i in range(LayerCount):
    layers.append(data_source.GetLayer(i).GetName())
print ("Layers in OSM File : " , layers)
In [ ]:
#GeoJSON Example
if os.path.exists(geojson_path):
    print ( 'Using Already Converted Files in '+ geojson_path + ' Folder' )
else :
    print( 'Converting File ' + pbf_path + " to " + geojson_path )
    os.makedirs(geojson_path)
    for layer in layers:
        dir(os.system('ogr2ogr -f "GeoJSON" ' + geojson_path + layer + '.geojson ' + pbf_path + ' ' + layer))
        print ('Done Converting File ' + pbf_path + " to " +  geojson_path + layer + '.geojson ')
In [ ]:
#Intialise master
df = None

#Get the layers in master
for layertype in layertypes:
    #spatialReference
    WKT_SR = data_source.GetLayer(layertype).GetSpatialRef().ExportToWkt()
    layer_file = geojson_path + layertype + '.geojson'

    #get GeoDataFrame from geojson
    gdf = gpd.read_file(layer_file)
    gdf.crs = WKT_SR
    
    #Filter Layer
    gdf = gdf[gdf['place'].isin(filter_type)]
    
    #If layertype is point buffer it
    if (layertype == "points"):
        #Buffer Layer
        gdf['geometry'] = gdf.geometry.buffer(buffer_point)
        #gdf.to_file(driver='ESRI Shapefile', filename='data.shp')
    elif (layertype == "multipolygons"):
        gdf['osm_id'] = gdf['osm_way_id']
    #Join to df
    if df is not None:
        df = pd.concat([df, gdf], axis=0)
    else:
        df = gdf

#Drop Duplicates
df = df.drop_duplicates(subset=['name'], keep=False)

#Drop Villages with no Name
df = df[df['name'] != None]
In [ ]:
route = {
    "type": "FeatureCollection",
    "name": "route",
    "crs": gdf.crs,
    "features": [
        { "type": "Feature", "properties": { }, "geometry": "" }
    ]
}

#limits for filtering reults
limit_dist_high = 7000 #in meters
limit_dist_low = 100 #in meters
def check_direct_localities(source, geometry):
    #get centroid
    centroid = str(geometry.centroid.x) + "," + str(geometry.centroid.y) 
    
    #If Far away dont call the Api
    src = source.split(",")
    distance_src_des = math.sqrt((float(src[0].strip()) - geometry.centroid.x)**2 + (float(src[1].strip()) - geometry.centroid.y)**2)
    distance_src_des = distance_src_des*degree_to_meter
    if ( (distance_src_des > limit_dist_high) or (distance_src_des < limit_dist_low)):
        return 0
    
    #Call Api
    url = '/route/v1/car/'+ source + ';' + centroid + '?alternatives=0&overview=full&steps=false&annotations=false&geometries=geojson'
    response = requests.get(osrm_server + url) #Get Route
    
    #Parse Response
    try :
        response = response.json() # Route to Json
        distance = response['routes'][0]['distance'] #Route Distance  
        if ( (distance > limit_dist_high) or (distance < limit_dist_low)):
            return 0
        #print ('distance:', distance)
        route["features"][0]["geometry"] = response['routes'][0]['geometry'] #convert to json
        route_df = gpd.GeoDataFrame.from_features(route['features']) #convert to GeoDataFrame
        route_df.crs = route['crs']
        if (gpd.sjoin(df, route_df, how="inner", op='intersects').shape[0] > 2): #check if more localities are on route
            return 0
        else :
            return 1
    except:
        return 0

def direct_localities(polygon):
    source = str(polygon.centroid.x) + "," + str(polygon.centroid.y)
    df['direct_route'] = df['geometry'].apply(lambda row : check_direct_localities(source, row)) #localities are directly connected
    return df[df['direct_route'] > 0][['name', 'osm_id']].values
In [ ]:
#df1 = df.head(10)
df['Direct_Localities'] = df['geometry'].apply(lambda x : direct_localities(x))
#df[df.index == 18][['name', 'osm_id', 'Direct_Localities']]
In [ ]:
#Final Output
'''
index_req = 18
values = df[df.index == index_req][['name', 'osm_id', 'Direct_Localities']].values
'''
In [ ]:
#Final Output
for index, row in df.iterrows():
    print_text = row['name'] + "," + str(row['geometry'].centroid.x) + "," + str(row['geometry'].centroid.y) + "," + str(len(row['Direct_Localities']))
    print(print_text)