Hello, World! Welcome to my first Python project.
In this project, we are going to look at OSMnx, a Python library that allows to download spatial geometries and construct, project, analyze, and visualize street networks using OpenStreetMap data. In addition, there is also a scientific articles available describing the package:
Boeing, G. 2017. “OSMnx: New Methods for Acquiring, Constructing, Analyzing, and Visualizing Complex Street Networks.” Computers, Environment and Urban Systems 65, 126-139. doi:10.1016/j.compenvurbsys.2017.05.004
Finally, I'd like to thank Mr Ujaval Gandhi from Spatial Thoughts, that encouraged this activity after attending his great course Python foundation for spatial analysis
# import necessary modules
import osmnx as ox
import networkx as nx
import geopandas as gpd
import matplotlib.pyplot as plt
# command to display matplotlib plots inline within the ipython notebook
%matplotlib inline
Let's start by specifying the place from where to download the data. For this project, we download data from Crescenzago/Padova district in Milan, Italy.
Note: the place name should be geocodable. To find the proper name, you can use Nominatim which is a debugging interface for the search engine used on the OpenStreetMap website.
Tag: "Padova, Municipio 2, Milano, Lombardia, Italia"
# Specify the name that is used to seach for the data
place_name = "Padova, Municipio 2, Milano, Lombardia, Italia"
OSMnx lets you download street network data and build topologically-corrected street networks, project and plot the networks, and save the street network as geopackage or shapefiles for later use.
You can specify several different network types:
Or simplify it, which will clean up nodes that are not intersections or dead-ends while retaining the complete edge geometry.
For a more in-depth demonstration of creating street networks, see OSMnx examples.
# Fetch OSM street network from the location
graph = ox.graph_from_place(place_name, simplify=True, network_type= 'all_private')
# Check data type of the graph
type(graph)
#DiGraph is a data type that stores nodes and edges with optional data, or attributes.
networkx.classes.multidigraph.MultiDiGraph
#Project graph from its current CRS to UTM.
graph = ox.project_graph(graph)
#Plot the graph nodes and edges
fig, ax = ox.plot_graph(graph, edge_linewidth=0.5, node_size=10, node_color="#66cc66")
OSMnx lets you download spatial "place boundary" geometries from OpenStreetMap, save them to shapefiles, project them, and plot them. For a more in-depth demonstration of creating these shapefiles, see OSMnx examples.
Let's also plot the Polygon that represents our area of interest. We can retrieve the polygon geometry using the geocode_to_gdf() function
# returns a GeoDataFrame based on the specified place name query
area = ox.geocode_to_gdf(place_name)
# Check the data values
area.head()
geometry | bbox_north | bbox_south | bbox_east | bbox_west | place_id | osm_type | osm_id | lat | lon | display_name | class | type | importance | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | POLYGON ((9.21711 45.50097, 9.21765 45.50050, ... | 45.511538 | 45.492533 | 9.256124 | 9.217106 | 198748861 | way | 404140221 | 45.501935 | 9.230626 | Padova, Municipio 2, Milan, Lombardy, Italy | place | quarter | 0.56 |
# get the polygon for Padova district, project it to UTM, and plot it
area_proj = ox.project_gdf(area)
ax = area_proj.plot (fc="grey", ec="none")
_ = ax.axis("off")
It is also possible to download other types of OSM data features with OSMnx such as buildings or points of interest (POIs). Let’s download the buildings with OSMnx geometries_from_place() function and plot them.
buildings = ox.geometries_from_place(place_name, tags={'building':True})
# Check how many building footprints we received
len(buildings)
1442
buildings.head()
addr:city | addr:country | addr:housenumber | addr:street | geometry | addr:postcode | source | name | access | amenity | ... | height | start_date | contact:facebook | contact:instagram | contact:linkedin | tower:type | building:use | grades | ways | type | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
element_type | osmid | |||||||||||||||||||||
node | 9297290741 | Milano | NaN | 14 | Via Rovigno | POINT (9.21947 45.50296) | 20125 | NaN | Bocciofila Martesana | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
way | 25666269 | NaN | NaN | NaN | NaN | POLYGON ((9.23091 45.50313, 9.23094 45.50312, ... | NaN | UDB60 | Anfiteatro Martesana | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
25678877 | NaN | NaN | NaN | NaN | POLYGON ((9.22938 45.49910, 9.22941 45.49897, ... | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
25678878 | NaN | NaN | NaN | NaN | POLYGON ((9.23009 45.49913, 9.23005 45.49928, ... | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
25678879 | NaN | NaN | NaN | NaN | POLYGON ((9.23010 45.49931, 9.23006 45.49947, ... | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
5 rows × 71 columns
buildings.columns
Index(['addr:city', 'addr:country', 'addr:housenumber', 'addr:street', 'geometry', 'addr:postcode', 'source', 'name', 'access', 'amenity', 'operator', 'brand', 'brand:wikidata', 'brand:wikipedia', 'tourism', 'website', 'building', 'leisure', 'sport', 'nodes', 'building:levels', 'created_by', 'isced:level', 'office', 'roof:levels', 'contact:phone', 'ref:vatin', 'shop', 'roof:shape', 'email', 'fax', 'phone', 'rooms', 'stars', 'addr:housename', 'denomination', 'religion', 'wikidata', 'man_made', 'name:zh', 'air_conditioning', 'opening_hours', 'payment:cash', 'payment:credit_cards', 'payment:debit_cards', 'wheelchair', 'abandoned', 'roof:colour', 'building:colour', 'building:material', 'wikipedia', 'internet_access', 'operator:type', 'ref:isil', 'old_name', 'cuisine', 'layer', 'ref:mise', 'outdoor_seating', 'smoking', 'architect', 'height', 'start_date', 'contact:facebook', 'contact:instagram', 'contact:linkedin', 'tower:type', 'building:use', 'grades', 'ways', 'type'], dtype='object')
# project to UTM
buildings_proj= ox.project_gdf(buildings)
# Plot
fig, ax = ox.plot_footprints(buildings_proj, figsize=(5, 5), color='black', bgcolor='w')
_ = ax.axis("off")
# Retrieve metro stop
stops = ox.geometries_from_place(place_name, tags={"station": "subway"})
# How many stops do we have?
len(stops)
1
# available colums
stops.columns
Index(['description', 'name', 'network', 'operator', 'public_transport', 'railway', 'start_date', 'station', 'wheelchair', 'wikidata', 'wikipedia', 'geometry'], dtype='object')
# Select some useful columns and print
stop_name = ['name', 'description', 'station']
# Print
stops[stop_name]
name | description | station | ||
---|---|---|---|---|
element_type | osmid | |||
node | 9897478677 | Turro | Area STIBM: Mi1 | subway |
We downloaded Turro subway station from M1 Line.
We can now plot all these different OSM layers by using geopandas. However, the street network data is not a GeoDataFrame, but a graph object. OSMnx provides the function graph_to_gdfs() that can convert the graph into two separate GeoDataFrames where the first one contains the information about the nodes and the second one about the edge.
Let’s extract the nodes and edges from the graph as GeoDataFrames:
# Retrieve nodes and edges
nodes, edges = ox.graph_to_gdfs(graph)
nodes.head()
y | x | highway | street_count | lon | lat | geometry | |
---|---|---|---|---|---|---|---|
osmid | |||||||
27551137 | 5.038073e+06 | 518401.541084 | crossing | 4 | 9.235517 | 45.495944 | POINT (518401.541 5038072.902) |
2533539843 | 5.038070e+06 | 518398.892171 | NaN | 3 | 9.235483 | 45.495921 | POINT (518398.892 5038070.305) |
2533541845 | 5.038103e+06 | 518415.914145 | NaN | 3 | 9.235703 | 45.496218 | POINT (518415.914 5038103.362) |
2533542537 | 5.038074e+06 | 518413.938024 | crossing | 4 | 9.235676 | 45.495953 | POINT (518413.938 5038073.849) |
27551139 | 5.038721e+06 | 518985.296557 | NaN | 3 | 9.243014 | 45.501765 | POINT (518985.297 5038721.309) |
edges.head()
osmid | highway | oneway | reversed | length | geometry | name | maxspeed | lanes | width | junction | access | service | bridge | tunnel | |||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
u | v | key | |||||||||||||||
27551137 | 2533539843 | 0 | 369654598 | footway | False | False | 3.706 | LINESTRING (518401.541 5038072.902, 518398.892... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
2533541845 | 0 | 8421642 | unclassified | True | False | 35.300 | LINESTRING (518401.541 5038072.902, 518402.028... | Via Palmanova | 30 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
2533542537 | 0 | [369654597, 369654598] | footway | False | True | 22.011 | LINESTRING (518401.541 5038072.902, 518410.804... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
2533539843 | 2533539844 | 0 | 369654598 | footway | False | False | 8.485 | LINESTRING (518398.892 5038070.305, 518392.705... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
27551137 | 0 | 369654598 | footway | False | True | 3.706 | LINESTRING (518398.892 5038070.305, 518401.541... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
Let’s create a map of Crescenzago/Padova district using the area, streets and buildings. We can plot in different ways
fig, ax = plt.subplots(1, 3, figsize=(18,4))
ax0, ax1, ax2 = ax
# Plot the area
area.plot(ax=ax0, facecolor='black')
fig
<Figure size 432x288 with 0 Axes>
# Plot the network edges
edges.plot(ax=ax1 ,linewidth=1, edgecolor='#BC8F8F')
fig
<Figure size 432x288 with 0 Axes>
# Plot the buildings
buildings.plot(ax=ax2, facecolor='khaki')
fig
<Figure size 432x288 with 0 Axes>
# Remove axis and add titles
ax0.axis('off')
ax0.set_title('Padova District area', y=-0.1)
ax1.axis('off')
ax1.set_title('Padova District network', y=-0.1)
ax2.axis('off')
ax2.set_title('Padova District buildings', y=-0.1)
fig
from pyproj import CRS
# Set projection
projection = CRS.from_epsg(32632)
# Re-project layers
area = area.to_crs(projection)
edges = edges.to_crs(projection)
buildings = buildings.to_crs(projection)
stops = stops.to_crs(projection)
# Create a new plot with the re-projected layers overlapped:
fig, ax = plt.subplots(figsize=(18,4))
ax.axis('off')
# Plot the area
area.plot(ax=ax, facecolor='#b1b1b1')
# plot the network
edges.plot(ax=ax, linewidth=0.8, edgecolor='white')
# plot the buildings
buildings.plot(ax=ax, facecolor='#6e6e6e')
#plot the statio
stops.plot(ax=ax, color='red')
plt.tight_layout()
Finally, we can export our data in different formats, such as Geopackage, to use in GIS.
# save street network as geopackage to work with it in GIS
ox.save_graph_geopackage(graph, filepath="./data/ graph.gpkg")
# save the buildings, area, and stops as a geopackage to work in GIS
buildings.geometry.to_file("./data/padova.gpkg", driver='GPKG', layer='buildings')
area.geometry.to_file("./data/padova.gpkg", driver='GPKG', layer='area')
stops.geometry.to_file("./data/padova.gpkg", driver='GPKG', layer='stops')
In this project, we explored a very cool library, OSMnx, which makes working with the OpenStreetMap source very easy. Additionally, I suggest you read the official documentation to see the library's full functionality.
Thanks!