In [1]:
%matplotlib inline
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
from shapely.geometry import shape
In [2]:
#config
data_folder = 'data/'
layertype = "multipolygons" #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 + 'romania-latest.osm.pbf' #location of PBF file
osrm_server = 'http://192.168.0.5:5000' #osrm
In [3]:
#Define Input datasource
driver = ogr.GetDriverByName('PBF') # 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 PBF File : " , layers)

Layers in PBF File :  ['points', 'lines', 'multilinestrings', 'multipolygons', 'other_relations']
In [4]:
#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 ')

Using Already Converted Files in data/geojson/ Folder
In [5]:
#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

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-5-bf0851e35875> in <module>()
      4 
      5 #get GeoDataFrame from geojson
----> 6 gdf = gpd.read_file(layer_file)
      7 gdf.crs = WKT_SR

C:\ProgramData\Anaconda3\lib\site-packages\geopandas\io\file.py in read_file(filename, **kwargs)
     18     bbox = kwargs.pop('bbox', None)
     19     with fiona.open(filename, **kwargs) as f:
---> 20         crs = f.crs
     21         if bbox is not None:
     22             assert len(bbox)==4

C:\ProgramData\Anaconda3\lib\site-packages\fiona\collection.py in crs(self)
    201         """Returns a Proj4 string."""
    202         if self._crs is None and self.session:
--> 203             self._crs = self.session.get_crs()
    204         return self._crs
    205 

KeyboardInterrupt:
In [177]:
gdf1 = gdf[gdf['admin_level'] == '6'].head(10)
gdf1['geometry'].apply(lambda x : direct_localities(x))
Out[177]:
index_left
103     [1207838, 2021225, 2835502, 3277038, 3456473, ...
2444    [1207838, 2021225, 2835502, 3277038, 3456473, ...
2960                                                   []
3612    [1207838, 2021225, 2835502, 3277038, 3456473, ...
3702    [1207838, 2021225, 2835502, 3277038, 3456473, ...
3704    [1207838, 2021225, 2835502, 3277038, 3456473, ...
3706                                                   []
4273    [1207838, 2021225, 2835502, 3277038, 3456473, ...
4536    [1207838, 2021225, 2835502, 3277038, 3456473, ...
5101    [1207838, 2021225, 2835502, 3277038, 3456473, ...
Name: geometry, dtype: object
In [176]:
route = {
    "type": "FeatureCollection",
    "name": "route",
    "crs": gdf.crs,
    "features": [
        { "type": "Feature", "properties": { }, "geometry": "" }
    ]
}

def check_direct_localities(source, geometry):
    centroid = str(geometry.centroid.x) + "," + str(geometry.centroid.y) #get centroid 
    url = '/route/v1/car/'+ source + ';' + centroid + '?alternatives=0&overview=full&steps=false&annotations=false&geometries=geojson'
    response = requests.get(osrm_server + url) #Get Route
    try :
        response = response.json() # Route to Json
        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(gdf1, 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)
    gdf1['direct_route'] = gdf1['geometry'].apply(lambda row : check_direct_localities(source, geometry)) #localities are directly connected
    return gdf1[gdf1['direct_route'] == 1]['osm_id'].values
In [ ]: