In [29]:
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
import altair as alt
In [30]:
pd.options.display.max_columns = 999
In [31]:
plt.rcParams['figure.figsize'] = (10,6)
Week 9A: OpenStreetMap, Urban Networks, and Interactive Web Maps¶
Nov 1, 2021
Happy day after Halloween¶
Housekeeping¶
- HW #4 (optional) due today by the end of the day
- HW #5 (optional) posted last Wednesday — due on Wed Nov. 10
- Final project details posted: https://github.com/MUSA-550-Fall-2021/final-project
Important: Update your local environment¶
- Small update to course's Python environment relevant to today's lecture
- Update the environment on your laptop using these instructions on course website
Today: OpenStreetMap (OSM)¶
- Two tools that make working with OSM data very easy
- What kind of questions can we answer?
- Street orientations
- Mapping event points to streets: car crashes
- Mapping amenities
- Network-constrained distances: accessibility
OSM: what is it?¶
- Collaborative mapping
- A free editable map of the World
- Sort of like Wikipedia for maps
Great source of data: street networks and a wealth of amenity information
Working with OSM data¶
- Raw data is very messy
- Two relatively new, amazing Python packages greatly simply the process
- Related, but complementary features
- OSMnx: downloading and manipulating streets as networks
- Pandana: networks focused on accessibility of amenities
Related: interactive web maps in Python
Part 1: OSMnx¶
Several key features:
- Downloading political boundaries for cities, states, countries, etc
- Downloading street networks
- Analyzing networks: routing, visualization, statistics
In [32]:
import osmnx as ox
In [33]:
philly = ox.geocode_to_gdf('Philadelphia, PA, USA')
philly.head()
Out[33]:
We can plot it just like any other GeoDataFrame
In [34]:
# Project it to Web Mercator first and plot
ax = philly.to_crs(epsg=3857).plot(facecolor='none', edgecolor='black')
ax.set_axis_off()
1.2 Projecting with OSMnx¶
Key function: project_gdf() (docs)
Automatically projects to the Universal Transverse Mercator (UTM) CRS for the UTM zone that the centroid of the polygon lies in
A good, general projection that works for most latitudes except very northern locations.
In [35]:
ax = ox.project_gdf(philly).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
Some more examples:
In [36]:
# Some examples
place1 = ox.geocode_to_gdf('Manhattan, New York City, New York, USA')
place2 = ox.geocode_to_gdf('Miami-Dade County, Florida')
place3 = ox.geocode_to_gdf('Florida, USA')
place4 = ox.geocode_to_gdf('Spain')
In [37]:
# Manhattan
ax = ox.project_gdf(place1).plot(fc="lightblue", ec="gray");
ax.set_axis_off()
In [38]:
# Miami-Dade County
ax = ox.project_gdf(place2).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
In [39]:
# Florida
ax = ox.project_gdf(place3).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
In [40]:
# Spain
ax = ox.project_gdf(place4).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
1.3 Downloading OSM features¶
Key functions: geometries_from_*
geometries_from_place()(docs)- Download geometries within an OSM place boundary
geometries_from_address()(docs)- Download geometries within a certain distance of an address
geometries_from_bbox()(docs)- Download geometries within a N, S, E, W bounding box
geometries_from_point()(docs)- Download geometries within a certain distance of a specified point
geometries_from_polygon()(docs)- Download geometries within a polygon object
About OSM features¶
Important reference: https://wiki.openstreetmap.org/wiki/Map_Features
- OSM uses a tagging system to categorize different map features
- The main feature categories are available on the OSM Wikipedia
- Examples: 'amenity', 'building', 'landuse', 'highway'
- There are specific sub-categories for each feature type too:
- Amenity examples: 'bar', 'college', 'library'
In the language of OSM, the "key" is the main feature category (e.g., amenity) and the "value" is the sub-category type (e.g., "bar")
OSMnx mirrors the key/value syntax¶
Use a dict to specify the features you want:
In [41]:
# Get all amenities in Philadelphia
amenities = ox.geometries_from_place("Philadelphia, PA", tags={"amenity": True})
len(amenities)
Out[41]:
In [42]:
amenities.head()
Out[42]:
In [43]:
# Get all bars in philadelphia
bars = ox.geometries_from_place("Philadelphia, PA", tags={"amenity": "bar"})
len(bars)
Out[43]:
In [44]:
bars.head()
Out[44]:
In [45]:
# Get bar, pub, and restaurant features in Philadelphia
food_and_drink = ox.geometries_from_place("Philadelphia, PA", tags={"amenity": ["pub", "bar", "restaurant"]})
len(food_and_drink)
Out[45]:
In [46]:
food_and_drink.head()
Out[46]:
In [47]:
# Get higway bus stop features
bus_stops = ox.geometries_from_place("Philadelphia, PA", tags={"highway": "bus_stop"})
len(bus_stops)
Out[47]:
In [48]:
bus_stops.head()
Out[48]:
In [49]:
# Get commercial and retail landuse features
landuse = ox.geometries_from_place("Philadelphia, PA", tags={"landuse": ["commercial", "retail"]})
len(landuse)
Out[49]:
In [50]:
fig, ax = plt.subplots(figsize=(10, 6))
ax = landuse.plot(ax=ax)
ax.set_axis_off()
1.4 Downloading street networks¶
Key functions: graph_from_*
graph_from_place()(docs)- Download street network within an OSM place boundary
graph_from_address()(docs)- Download street network within a certain distance of an address
graph_from_bbox()(docs)- Download street network within a N, S, E, W bounding box
graph_from_point()(docs)- Download street network within a certain distance of a specified point
graph_from_polygon()(docs)- Download street network within a polygon object
Street network around an address¶
Get streets within 500 meters of the center of Northern Liberties
In [51]:
G = ox.graph_from_address("Northern Liberties, Philadelphia, PA", dist=500)
Project and plot it:
In [52]:
G_projected = ox.project_graph(G)
ox.plot_graph(G_projected);
Remove the nodes:
In [53]:
ox.plot_graph(G_projected, node_size=0);
Let's zoom out to 2,000 meters. This will take a little longer.
In [54]:
G = ox.graph_from_address("Northern Liberties, Philadelphia, PA", dist=2000)
G_projected = ox.project_graph(G)
In [55]:
ox.plot_graph(G_projected, node_size=0);
Getting different network types¶
drive- get drivable public streets (but not service roads)drive_service- get drivable streets, including service roadswalk- get all streets and paths that pedestrians can use (this network type ignores one-way directionality)bike- get all streets and paths that cyclists can useall- download all non-private OSM streets and pathsall_private- download all OSM streets and paths, including private-access ones (default)
In [56]:
# the "drive" network
G = ox.graph_from_address("Northern Liberties, Philadelphia, PA",
network_type='drive')
ox.plot_graph(G);
In [57]:
# the "walk" network
G = ox.graph_from_address("Northern Liberties, Philadelphia, PA",
network_type='walk')
ox.plot_graph(ox.project_graph(G));
Street network within a place boundary¶
Use graph_from_place() to get the streets within a specific OSM place.
Note: the place query has to be resolved by OSM.
In [58]:
berkeley = ox.graph_from_place("Berkeley, California, USA",)
In [59]:
ox.plot_graph(ox.project_graph(berkeley), node_size=0);
Streets within a specific polygon¶
Example: all streets within Northern Liberties and Fishtown
First, let's wrangle some Zillow neighborhood boundaries¶
In [60]:
url = "https://github.com/azavea/geo-data/raw/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson"
hoods = gpd.read_file(url).rename(columns={"mapname": "neighborhood"})
In [61]:
hoods.head()
Out[61]:
Trim to Fishtown and Northern Liberties
In [62]:
sel = hoods['neighborhood'].isin(['Fishtown - Lower Kensington', 'Northern Liberties'])
nolibs_fishtown = hoods.loc[sel]
In [63]:
nolibs_fishtown.head()
Out[63]:
In [64]:
ax = ox.project_gdf(nolibs_fishtown).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
Extract streets within these polygons¶
- Take the union of the polygons:
unary_union - Use
ox.graph_from_polygon()
In [65]:
nolibs_fishtown_outline = nolibs_fishtown.geometry.unary_union
nolibs_fishtown_outline
Out[65]:
In [66]:
type(nolibs_fishtown_outline)
Out[66]:
In [67]:
# get the graph
G_nolibs_fishtown = ox.graph_from_polygon(nolibs_fishtown_outline, network_type='drive')
In [68]:
# viola!
ox.plot_graph(ox.project_graph(G_nolibs_fishtown), node_size=0);
We could also use unary_union.convex_hull. This will be an encompassing polygon around any set of geometries.
In [69]:
nolibs_fishtown.geometry.unary_union
Out[69]:
In [70]:
nolibs_fishtown.geometry.unary_union.convex_hull
Out[70]:
In [71]:
# only get the edges
nolibs_edges = ox.graph_to_gdfs(G_nolibs_fishtown,
edges=True,
nodes=False)
In [72]:
# we have lots of data associated with each edge!
nolibs_edges.head()
Out[72]:
In [73]:
# plot it like any old GeoDataFrame
ax = nolibs_edges.to_crs(epsg=3857).plot(color='gray')
# add the neighborhood boundaries
boundary = gpd.GeoSeries([nolibs_fishtown_outline], crs='EPSG:4326')
boundary.to_crs(epsg=3857).plot(ax=ax, facecolor='none', edgecolor='red', linewidth=3, zorder=2)
ax.set_axis_off()
1.6 What can we do with the graph?¶
- Network-based statistics
- Routing
- Street orientations
- Visualizing crashes
And much more: see the OSMnx repository of Jupyter notebook examples
Network statistics¶
We can use the ox.basic_stats() to get some basic network statistics
In [74]:
ox.basic_stats(G_nolibs_fishtown)
Out[74]:
In [75]:
import networkx as nx
Let's calculate the shortest route between Frankford Hall (a bar in Fishtown) and the Spring Garden subway station:
In [76]:
# Get all food and drink places within our Fishtown/No Libs polygon
fishtown_food_drink = ox.geometries_from_polygon(nolibs_fishtown_outline,
tags={"amenity": ["restaurant", "bar", "pub"]})
In [77]:
fishtown_food_drink.head()
Out[77]:
In [78]:
# Select Frankford Hall based on the name column
frankford_hall = fishtown_food_drink.query("name == 'Frankford Hall'").squeeze()
# Get the x/y coordinates
frankford_hall_x = frankford_hall.geometry.x
frankford_hall_y = frankford_hall.geometry.y
In [79]:
# Get all subway stations within our Fishtown/No Libs polygon
fishtown_subway_stops = ox.geometries_from_polygon(nolibs_fishtown_outline,
tags={"railway": "station"})
In [80]:
# Select Spring Garden based on the name column
spring_garden_station = fishtown_subway_stops.query("name == 'Spring Garden'").squeeze()
# Get the x/y coordinates
spring_garden_x = spring_garden_station.geometry.x
spring_garden_y = spring_garden_station.geometry.y
Find the nearest nodes on our OSMnx graph:
In [81]:
# Get the origin node
orig_node = ox.distance.nearest_nodes(G_nolibs_fishtown, spring_garden_x, spring_garden_y)
# Get the destination node
dest_node = ox.distance.nearest_nodes(G_nolibs_fishtown, frankford_hall_x, frankford_hall_y)
Use networkx to find the shortest path between these graph nodes
In [82]:
# get the shortest path --> just a list of node IDs
route = nx.shortest_path(G_nolibs_fishtown,
orig_node,
dest_node,
weight='length')
route
Out[82]:
Use ox.plot_graph_route() to plot a graph and highlight a specific route
In [83]:
ox.plot_graph_route(G_nolibs_fishtown, route, node_size=0);
Part 2: Pandana¶
"Pandas Network Analysis - dataframes of network queries, quickly"
A complementary set of OSM-related features:
- Downloading OSM-based networks
- Extracting amenity data (so-called "Points of Interest")
- Calculating network-constrained distances
In [84]:
import pandana as pnda
from pandana.loaders import osm
Step 1: Get amenity data¶
Key function: osm.node_query()
- This will extract amenities within a given bounding box.
- Similar to the
ox.geometries_from_bbox()function in OSMnx, but we slightly different syntax.
In [85]:
osm.node_query?
Get the bounding box for Northern Liberties / Fishtown:
In [86]:
boundary = nolibs_fishtown_outline.bounds
boundary
Out[86]:
In [87]:
[lng_min, lat_min, lng_max, lat_max] = boundary
In [88]:
# query OSM
poi_df = osm.node_query(lat_min, lng_min, lat_max, lng_max)
# remove missing data
poi_df = poi_df.dropna(subset=['amenity'])
In [89]:
poi_df.head()
Out[89]:
In [90]:
len(poi_df)
Out[90]:
In [91]:
poi_df[['lat', 'lon', 'amenity']].head(10)
Out[91]:
Explore the amenities in this region¶
For the full list of amenities, see the OSM Wikipedia
In [92]:
chart = (
alt.Chart(poi_df)
.mark_bar()
.encode(y=alt.Y("amenity", sort="-x"), x="count()", tooltip=["amenity", "count()"])
)
chart
Out[92]:
Step 2: Create a Pandana network¶
- Key function:
pdna_network_from_bbox() - It takes a bounding box and returns the OSM network within that region.
- Multiple network types: 'walk' and 'drive'
In [93]:
net = osm.pdna_network_from_bbox(
lat_min, lng_min, lat_max, lng_max, network_type="walk"
)
Step 3: Tell the network the location of amenities¶
Key function: network.set_pois()
- Today, we'll explore these four amenities: "restaurant", "bar", "school", "car_sharing"
- IMPORTANT: if you want to explore other amenities, you'll need to run the code below for your amenities of interest
In [94]:
# sensible defaults
max_distance = 2000 # in meters
num_pois = 10 # only need the 10 nearest POI to each point in the network
AMENITIES = ["restaurant", "bar", "school", "car_sharing"]
for amenity in AMENITIES:
# get the subset of amenities for this type
pois_subset = poi_df[poi_df["amenity"] == amenity]
# set the POI, using the longitude and latitude of POI
net.set_pois(
amenity, max_distance, num_pois, pois_subset["lon"], pois_subset["lat"]
)
In [95]:
# keyword arguments to pass for the matplotlib figure
bbox_aspect_ratio = (lat_max - lat_min) / (lng_max - lng_min)
fig_kwargs = {"facecolor": "w", "figsize": (10, 10 * bbox_aspect_ratio)}
# keyword arguments to pass for scatter plots
plot_kwargs = {"s": 20, "alpha": 0.9, "cmap": "viridis_r", "edgecolor": "none"}
Step 4: Plot the walking distance to the nearest POI¶
For every point on the network, find the nth nearest POI, calculate the distance, and color that point according to the distance.
- Use
network.nearest_poi()to get distances from nodes to nearest POIs - Merge coordinates of network nodes with distances to nearest POIs
- Plot the node coordinates, colored by distance to nth nearest POI
1. Use network.nearest_poi() to get distances from nodes to nearest POIs¶
In [96]:
num_pois
Out[96]:
In [97]:
amenity = 'bar'
access = net.nearest_pois(distance=1000,
category=amenity,
num_pois=num_pois)
In [98]:
access.head(n=20)
Out[98]:
2. Merge coordinates of network nodes with distances to nearest POIs¶
In [99]:
net.nodes_df.head()
Out[99]:
In [100]:
access.head()
Out[100]:
In [101]:
# Merge the nodes and the distance to POIs
nodes = pd.merge(net.nodes_df, access, left_index=True, right_index=True)
# Make into a geodataframe
nodes = gpd.GeoDataFrame(
nodes, geometry=gpd.points_from_xy(nodes["x"], nodes["y"]), crs="EPSG:4326"
)
In [102]:
nodes.head()
Out[102]:
And now plot it!¶
Let's define a function to do this for us, since we'll repeat this plot multiple times
In [103]:
def plot_walking_distance(net, amenity, distance=1000, n=1):
"""
Plot the walking distance to the specified amenity
"""
from mpl_toolkits.axes_grid1 import make_axes_locatable
# subset of POI
poi_subset = poi_df[poi_df["amenity"] == amenity]
# get the distances to nearest num_pois POI
access = net.nearest_pois(distance=1000, category=amenity, num_pois=num_pois)
# merge node positions and distances to nearest PO
nodes = pd.merge(net.nodes_df, access, left_index=True, right_index=True)
nodes = gpd.GeoDataFrame(
nodes, geometry=gpd.points_from_xy(nodes["x"], nodes["y"]), crs="EPSG:4326"
)
# Create the figure
fig, ax = plt.subplots(figsize=(10, 10))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
# plot the distance to the nth nearest amenity
ax = nodes.plot(ax=ax, cax=cax, column=nodes[n], legend=True, **plot_kwargs)
# add the amenities as stars
for i, row in poi_subset.iterrows():
ax.scatter(row["lon"], row["lat"], color="red", s=100, marker="*")
# format
ax.set_facecolor("black")
ax.figure.set_size_inches(fig_kwargs["figsize"])
# set extent
[xmin, ymin, xmax, ymax] = nodes.geometry.total_bounds
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
return ax
Evaluating amenity choice¶
The difference between maps to the nearest amenity and for example, the 5th nearest amenity tells us about the options consumers have
Example: bars¶
In [104]:
ax = plot_walking_distance(net, "bar", n=1)
ax.set_title("Walking distance to the nearest bar", fontsize=18);
In [105]:
ax = plot_walking_distance(net, "bar", n=3)
ax.set_title("Walking distance to the 3rd nearest bar", fontsize=18);
Example: schools¶
In [106]:
ax = plot_walking_distance(net, "school", n=1)
ax.set_title("Walking distance to the nearest school", fontsize=18);
In [107]:
ax = plot_walking_distance(net, "school", n=3)
ax.set_title("Walking distance to the 3rd nearest school", fontsize=18);
Example: restaurants¶
In [108]:
ax = plot_walking_distance(net, "restaurant", n=1)
ax.set_title("Walking distance to the nearest restaurant", fontsize=18);
In [109]:
ax = plot_walking_distance(net, "restaurant", n=5)
ax.set_title("Walking distance to the 5th nearest restaurant", fontsize=18);
Example: car sharing¶
In [110]:
ax = plot_walking_distance(net, "car_sharing", n=1)
ax.set_title("Walking distance to the nearest car sharing", fontsize=18);
In [111]:
ax = plot_walking_distance(net, "car_sharing", n=5)
ax.set_title("Walking distance to the 5th nearest car sharing", fontsize=18);
At-home exercise: Explore amenities in the neighborhood of your choice¶
Many, many more amenities are logged throughout the city. Pick your favorite neighborhood and explore.
See this page for the full list of amenities.
Final project idea: With this kind of analysis, you can look at amenity-based influence in housing, neighborhood selection, etc. or something similar to the Walk Score.
In [ ]: