Solution 05

Setup

import pandas as pd
import requests

import warnings
warnings.filterwarnings('ignore')
from pandas.core.common import SettingWithCopyWarning
warnings.simplefilter(action="ignore", category=SettingWithCopyWarning)

try:
  import pygeos
except ModuleNotFoundError as e:
  !pip install pygeos==0.13
  import pygeos

try:
  import mapclassify
except ModuleNotFoundError as e:
  !pip install mapclassify
  import mapclassify

if mapclassify.__version__ != "2.4.3":
  !pip install -U mapclassify==2.4.3

try:
  import geopandas as gpd
except ModuleNotFoundError as e:
  !pip install geopandas==0.12.1
  import geopandas as gpd

if gpd.__version__ != "0.12.1":
  !pip install -U geopandas==0.12.1
  import geopandas as gpd

try:
  import pyrosm
except ModuleNotFoundError as e:
  !pip install pyrosm==0.6.1
  import pyrosm

try:
  import osmnx  as ox
except ModuleNotFoundError as e:
  !pip install osmnx==1.2.2
  import osmnx  as ox
if ox.__version__ != "1.2.2":
  !pip install -U osmnx==1.2.2
  import osmnx  as ox

import pandas as pd
from matplotlib import pyplot as plt

Exercise

identify the shortest path by walking to reach the Castle of Trento from the main train station of Trento
identify the streets network orientation of the cities: Trento-Italy, Verona-Italy, Munich-Germany, Athens-Greece
locate the student residences of Trento in OpenStreetMap and identify services in an area of 15 minutes based on the concept of Carlos Moreno of the 15 minute city
“…a way to ensure that urban residents can fulfill six essential functions within a 15-minute walk [or bike] from their dwellings: living, working, commerce, healthcare, education and entertainment…”

You are free to use OSMnx, Pyrosm with networkx or pandas or igraph

identify the shortest path by walk to reach the Castle of Trento from the main train station

where is the castle in Trento
where is the main train station

where is the castle in Trento

THe name of the castle of Trento is Castello del Buonconsiglio

from geopy.geocoders import Nominatim

geolocator = Nominatim(user_agent="geospatial course")
location = geolocator.geocode("Museo del Castello del Buonconsiglio")

point_castle =  (location.latitude, location.longitude)

where is the main train station in Trento?

https://en.wikipedia.org/wiki/Trento_railway_station

location = geolocator.geocode("Piazza Dante 9, Trento")

point_train_station =  (location.latitude, location.longitude)

download the data for PyrOSM

url_download_trento_pbf = 'https://osmit-estratti.wmcloud.org/dati/poly/comuni/pbf/022205_Trento_poly.osm.pbf'
import urllib.request
urllib.request.urlretrieve(url_download_trento_pbf ,"trento_osm.pbf")    
osm = pyrosm.OSM("trento_osm.pbf")

create the streets walking graph

nodes, edges = osm.get_network(network_type="walking", nodes=True)

with networkx

network_onx = osm.to_graph(nodes, edges, graph_type="networkx")
%time

CPU times: user 4 µs, sys: 1e+03 ns, total: 5 µs
Wall time: 7.63 µs

point_nearest_train_station = ox.distance.nearest_nodes(network_onx,Y=point_train_station[0],X=point_train_station[1])
point_nearest_castle = ox.distance.nearest_nodes(network_onx,Y=point_castle[0],X=point_castle[1])

calculate the shortest path

shortest_path = ox.shortest_path(network_onx, point_nearest_train_station, point_nearest_castle, weight="length")

show it on the map

from shapely.geometry import Point, LineString

# function to create the route from the ids of the nodes
def route_nodes_to_line_networkx(nodelist, network,solution=1):
  points = []
  for idnode in nodelist:
    lon = network.nodes[idnode]['x']
    lat = network.nodes[idnode]['y']
    point = Point(lon,lat)
    points.append(point)
  path = LineString(points)
  route = gpd.GeoDataFrame(
    {"src_node": [nodelist[0]], "tgt_node": [nodelist[-1]], "solution": [solution]},
      geometry=[path],
      crs="epsg:4326"
  )
  return route

route_with_onx = route_nodes_to_line_networkx(shortest_path, network_onx)

route_with_onx

	src_node	tgt_node	solution	geometry
0	2591984564	6036162196	1	LINESTRING (11.11953 46.07204, 11.11957 46.072...

route_with_onx.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

import fiona
fiona.drvsupport.supported_drivers['KML'] = 'rw'

route_with_onx.to_file("path_trainstation2castle_trento.kml", driver="KML")

.. and now you can visualize it with Google Earth

paths = ox.k_shortest_paths(network_onx, point_nearest_train_station, point_nearest_castle, k=10, weight="distance")

i = 0
gdf_paths = gpd.GeoDataFrame()
for path in paths:
    new_route = route_nodes_to_line_networkx(path, network_onx,solution=str(i))
    if i == 0:
        gdf_paths = new_route
    else:
        gdf_paths = gdf_paths.append(new_route)
    i += 1

gdf_paths.explore(column='solution',cmap="tab10",tiles="Stamen Toner")

Make this Notebook Trusted to load map: File -> Trust Notebook

gdf_paths[gdf_paths.solution == "0"].explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

gdf_paths[gdf_paths.solution == "1"].explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

gdf_paths[gdf_paths.solution == "2"].explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

with igraph

nodes, edges = osm.get_network(nodes=True, network_type="walking")

network_igraph = osm.to_graph(nodes, edges,graph_type='igraph')

from shapely.ops import nearest_points

# look code solutions of lesson 2
def get_nearest_id(id, points):
    # Create a union of points (multipoint geometry
    multipoints = points[points.id != id]["geometry"].unary_union
    # identify the starting point
    point = points[points.id == id]
    # find the nearest points
    nearest_geoms = nearest_points(point['geometry'].values[0], multipoints)
    # get corresponding values of the nearest point
    # note: in the position 0 there is the starting point
    nearest_data = points[points["geometry"] == nearest_geoms[1]]
    # extract the id of the nearest point
    nearest_id = nearest_data['id'].values[0]
    return (nearest_id)

# create geodaframe
node_ids = []
ids = []
points = []
for v in network_igraph.vs:
    node_id = v['node_id']
    id = v['id']
    p = Point(v['lon'],v['lat'])
    ids.append(id)
    points.append(p)
    node_ids.append(node_id)
data = {}
data['id'] = ids
data['node_id'] = node_ids
data['geometry'] = points
gdf_vertex = gpd.GeoDataFrame(data,crs="epsg:4326")    

multipoints = gdf_vertex.unary_union

point_trainstation = Point(point_train_station[1],point_train_station[0])
point_castle = Point(point_castle[1],point_castle[0])

nearest_geoms_train = nearest_points(point_trainstation, multipoints)
nearest_geoms_castle = nearest_points(point_castle, multipoints)

pstation = gdf_vertex[gdf_vertex.geometry == nearest_geoms_train[1]]
pcastle = gdf_vertex[gdf_vertex.geometry == nearest_geoms_castle[1]]

pstation

	id	node_id	geometry
39180	2591984564	39180	POINT (11.11953 46.07204)

idnodes_path = network_igraph.get_shortest_paths(pstation.node_id.values[0],to=pcastle.node_id.values[0])[0]

gdf_nodes_path = gdf_vertex[gdf_vertex.node_id.isin(idnodes_path)]

gdf_nodes_path.explore(marker_kwds={"color": "green", "radius": "10"})

Make this Notebook Trusted to load map: File -> Trust Notebook

identify the streets network orientation of this list of cities

Trento-Italy, Verona-Italy, Munich-Germany, Athens-Greece

The suggestion is to read the blog Street Network Orientation post of Geoff Boeing (author of OSMnx) and adapt also che script present in the examples section on the OSMnx github repository.

From the script you have simple to change the names of the cities … and .. understand the code.

# define the study sites as label : query
places = {
    'Trento' : "Trento, Territorio Val d'Adige, Provincia di Trento, Trentino-Alto Adige/Südtirol, Italy",
    "Verona": "Verona, Veneto, Italy",
    "Munich": "Munich, Bavaria, Germany",
    "Athens": "Περιφερειακή Ενότητα Κεντρικού Τομέα Αθηνών"
}

# verify OSMnx geocodes each query to what you expect (i.e., a [multi]polygon geometry)
gdf = ox.geocode_to_gdf(list(places.values()))
gdf

	geometry	bbox_north	bbox_south	bbox_east	bbox_west	place_id	osm_type	osm_id	lat	lon	display_name	class	type	importance
0	POLYGON ((11.02247 46.05027, 11.02296 46.05019...	46.153011	45.977531	11.194823	11.022474	297532930	relation	46663	46.066423	11.125760	Trento, Territorio Val d'Adige, Provincia di T...	boundary	administrative	1.932529
1	POLYGON ((10.87685 45.46084, 10.87712 45.46017...	45.541837	45.349440	11.123900	10.876851	297101238	relation	44830	45.438496	10.992412	Verona, Veneto, Italy	boundary	administrative	0.976157
2	MULTIPOLYGON (((11.36078 48.15807, 11.36085 48...	48.248116	48.061624	11.722910	11.360777	297759456	relation	62428	48.137108	11.575382	Munich, Bavaria, Germany	boundary	administrative	1.076211
3	POLYGON ((23.68700 37.97844, 23.68749 37.97825...	38.059134	37.915813	23.823200	23.687005	298197509	relation	2604796	37.987486	23.740311	Regional Unit of Central Athens, Attica, Greece	boundary	administrative	0.419111

gdf.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

https://www.oraask.com/wiki/draw-a-polar-histogram-in-matplotlib

# import numpy as np
# n = len(places)
# ncols = int(np.ceil(np.sqrt(n)))
# nrows = int(np.ceil(n / ncols))
# figsize = (ncols * 5, nrows * 5)
# fig, axes = plt.subplots(nrows, ncols, figsize=figsize, subplot_kw={"projection": "polar"})
# # plot each city's polar histogram

# for ax, place in zip(axes.flat, sorted(places.keys())):
#     print(ox.utils.ts(), place)

#     # get undirected graphs with edge bearing attributes
#     G = ox.graph_from_place(place, network_type="drive")
#     Gu = ox.add_edge_bearings(ox.get_undirected(G))
#     fig, ax = ox.bearing.plot_orientation(Gu, ax=ax, title=place, area=True)

# # add figure title and save image
# suptitle_font = {
#     "family": "Arial",
#     "fontsize": 60,
#     "fontweight": "normal",
#     "y": 1,
# }
# fig.suptitle("City Street Network Orientation") #, **suptitle_font)
# fig.tight_layout()
# fig.subplots_adjust(hspace=0.35)
# fig.savefig("street-orientations.png", facecolor="w", dpi=100, bbox_inches="tight")
# plt.close()

Here the result

.. and here the result for all the main cities of Italy

Vladimir Agafonkin created the web verson

https://mourner.github.io/road-orientation-map/#13.49/45.43646/12.33136

Here a simple map with the information of some cities in Trentino long the river Adige.

http://umap.openstreetmap.fr/en/map/orientamento-comuni-valle-delladige_234308#11/45.9588/11.1295

locate the student residences of Trento in OpenStreetMap and identify services in an area of 15 minutes based on the concept of Carlos Moreno of the 15 minute city

“…a way to ensure that urban residents can fulfill six essential functions within a 15-minute walk [or bike] from their dwellings: living, working, commerce, healthcare, education and entertainment…”

Data from OpenStreetMap

identify the tags for:

commerce:
- building = commercial, industrial, kiosk, retail, supermarket, warehouse
- shop = *
- amenity = bar, pub, cafe, restaurant, fast_food, food_court, ice_cream
healthcare
- amenity = baby_hatch,clinic,dentist, doctors, hospital, nursing_home, pharmacy, social_facility, veterinary
education
- amenity = college, driving_school, kindergarten, language_school, library, toy_library, training, music_school, university
entertainment
- amenity = arts_centre, cinema, community_centre, events_venue, fountain, social_facility, gambling, planetarium, public_bookcase, social_centre, studio, theatre
living
- amenity = courthouse, fire_station, police, post_box, post_office, ranger_station, townhall

buildings = ox.geometries.geometries_from_place("Trento, Territorio Val d'Adige, Provincia di Trento, Trentino-Alto Adige/Südtirol, Italy", tags = {'building': True})
amenities = ox.geometries.geometries_from_place("Trento, Territorio Val d'Adige, Provincia di Trento, Trentino-Alto Adige/Südtirol, Italy", tags = {'amenity': True})
shops = ox.geometries.geometries_from_place("Trento, Territorio Val d'Adige, Provincia di Trento, Trentino-Alto Adige/Südtirol, Italy", tags = {'shop': True})

amenity = student_accommodation

student_accomodations = amenities[amenities.amenity.isin(["student_accommodation"])]

transform in POI (with some bias)

student_accomodations['geometry'] = student_accomodations.geometry.representative_point()

create isocrhones for each student accomodations

import networkx as nx
G_proj = ox.project_graph(network_onx)
travel_speed = 4  # walking speed in km/hour

# add an edge attribute for time in minutes required to traverse each edge
meters_per_minute = travel_speed * 1000 / 60  # km per hour to m per minute
for orig,dest, p, data in G_proj.edges(data=True, keys=True):
    data["time"] = data["length"] / meters_per_minute

crs_proj = ox.graph_to_gdfs(G_proj)[0].crs

# create the isochrone polygons
isochrone_polys = []
names = []
for i in range(len(student_accomodations.geometry)):
    center_node = ox.distance.nearest_nodes(network_onx,Y=student_accomodations.geometry[i].y,X=student_accomodations.geometry[i].x)
    name = student_accomodations[student_accomodations.geometry==student_accomodations.geometry[i]].name.values[0]
    subgraph = nx.ego_graph(G_proj, center_node, radius=15, distance="time")
    node_points = [Point((data["x"], data["y"])) for node, data in subgraph.nodes(data=True)]
    bounding_poly = gpd.GeoSeries(node_points).unary_union.convex_hull
    isochrone_polys.append(bounding_poly)
    names.append(name)

data = {}
data['name'] = names
data['geometry'] = isochrone_polys
isochrones_student_accomodations = gpd.GeoDataFrame(data,crs=crs_proj)    

commerce

commerce = ["commercial","industrial", "kiosk", "retail", "supermarket", "warehouse"]
commerce_amenity = ["bar","pub","cafe","restaurant","fast_food","food_court","ice_cream"]

pois_commerce = buildings[buildings.building.isin(commerce)]
pois_shops = shops 
pois_commerce_amenities = amenities[amenities.amenity.isin(commerce_amenity)]

pois_shops.plot()
plt.show()

png

pois_commerce.plot()
plt.show()

png

pois_commerce.explore(tooltip=['building'])

Make this Notebook Trusted to load map: File -> Trust Notebook

pois_commerce_amenities.plot()
plt.show()

png

pois_commerce_amenities.explore(tooltip=['amenity'])

Make this Notebook Trusted to load map: File -> Trust Notebook

healthcare

healthcare_amenity_values = ["baby_hatch","clinic,dentist","doctors","hospital","nursing_home","pharmacy","social_facility","veterinary"]
pois_healthcare = amenities[amenities.amenity.isin(healthcare_amenity_values)]

pois_healthcare.plot()
plt.show()

png

education

education_amenity_values = ["college","driving_school","kindergarten","language_school","library","toy_library","training","music_school","university"]
pois_education = amenities[amenities.amenity.isin(education_amenity_values)]

pois_education.plot()
plt.show()

png

entertainment

entertainment_amenity_values = ["arts_centre","cinema","community_centre","events_venue","fountain","social_facility","gambling","planetarium","public_bookcase","social_centre","studio","theatre"]
pois_entertainment = amenities[amenities.amenity.isin(entertainment_amenity_values)]

pois_entertainment.plot()
plt.show()

png

living

living_amenity_values = ["courthouse","fire_station","police","post_box","post_office","ranger_station","townhall"]
pois_living = amenities[amenities.amenity.isin(living_amenity_values)]

pois_living.plot()
plt.show()

png

cross the area

isochrones_student_accomodations = isochrones_student_accomodations.to_crs(epsg=4326)

… you can use different approach
this is made only to repeat some concepts presented in the previous lessons

def checkInsideArea(p,i):
    touch = False
    try:
        touch = p.within(i)
    except Exception as e:
        print(e)
        pass
    if touch == False:
        return False
    else:
        return True

names = {}
for idx, row in isochrones_student_accomodations.iterrows():
    name = row['name']
    label = "student_accomodation_" + str(idx)
    names[label] = name
    polygon = row.geometry
    pois_living[label] = pois_living.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_entertainment[label] = pois_entertainment.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_education[label] = pois_education.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_healthcare[label] = pois_healthcare.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_commerce[label] = pois_commerce.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_shops[label] = pois_healthcare.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_commerce[label] = pois_commerce.geometry.apply(lambda g: checkInsideArea(g,row.geometry))
    pois_commerce_amenities[label] = pois_commerce_amenities.geometry.apply(lambda g: checkInsideArea(g,row.geometry))

for name in names:
    living_places = pois_living[pois_living[name] == True].shape[0]
    entertainment_places = pois_entertainment[pois_entertainment[name] == True].shape[0]
    education_places = pois_education[pois_education[name] == True].shape[0]
    healthcare_places = pois_healthcare[pois_healthcare[name] == True].shape[0]
    commerce_places = pois_commerce[pois_commerce[name] == True].shape[0]
    commerce_places += pois_shops[pois_shops[name] == True].shape[0]
    commerce_places += pois_commerce_amenities[pois_commerce_amenities[name] == True].shape[0]
    total = commerce_places + living_places + education_places + healthcare_places + entertainment_places
    print(names[name])
    print("living places %s" % living_places)
    print("entertainment places %s" % entertainment_places)
    print("education places %s" % education_places)
    print("healthcare places %s" % healthcare_places)
    print("commerce places %s" % commerce_places)
    print("TOTAL %s" % total)
    print(" ")

Studentato San Bartolameo
living places 0
entertainment places 2
education places 0
healthcare places 0
commerce places 1
TOTAL 3
 
Collegio Bernardo Clesio
living places 2
entertainment places 1
education places 5
healthcare places 0
commerce places 58
TOTAL 66
 
Studentato Mayer
living places 3
entertainment places 1
education places 3
healthcare places 1
commerce places 13
TOTAL 21
 
NEST
living places 1
entertainment places 1
education places 1
healthcare places 0
commerce places 5
TOTAL 8

Share on

Twitter Facebook LinkedIn

Geospatial Lab

Setup

Exercise

identify the shortest path by walk to reach the Castle of Trento from the main train station

where is the castle in Trento

where is the main train station in Trento?

download the data for PyrOSM

create the streets walking graph

with networkx

calculate the shortest path

show it on the map

with igraph

identify the streets network orientation of this list of cities

Trento-Italy, Verona-Italy, Munich-Germany, Athens-Greece

locate the student residences of Trento in OpenStreetMap and identify services in an area of 15 minutes based on the concept of Carlos Moreno of the 15 minute city

Data from OpenStreetMap

Share on