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
#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
#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)
#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 ')
#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]
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
#df1 = df.head(10)
df['Direct_Localities'] = df['geometry'].apply(lambda x : direct_localities(x))
#df[df.index == 18][['name', 'osm_id', 'Direct_Localities']]
#Final Output
'''
index_req = 18
values = df[df.index == index_req][['name', 'osm_id', 'Direct_Localities']].values
'''
#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)