In [1]:
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
np.random.seed(42)
In [2]:
plt.rcParams['figure.figsize'] = (10,6)
Lecture 10B: Clustering Analysis in Python¶
Nov 10, 2021
Housekeeping¶
- Assignment #5 (optional) due by the end of the day today
- Assignment #6 (optional) assigned today and due in two weeks (11/24)
- Covers urban street networks and osmnx
- Plotting a folium heatmap of a dataset of your choosing — similar to building permit example from last lecture
- Remember, you have to do one of homeworks #4, #5, or #6
Last lecture (10A)¶
- Making a Folium Choropleth the "hard way" — using folium.GeoJson()
- Folium heatmap of new construction permits in Philadelphia
- Introduction to clustering with scikit-learn
- K-means algorithm for non-spatial clustering
Today¶
- Exercise on non-spatial clustering with k-means: Airbnb data in Philly neighborhoods
- Spatial clustering with the DBSCAN algorithm
- Exercise on spatial clustering of NYC taxi trips
In [3]:
import altair as alt
from vega_datasets import data as vega_data
In [4]:
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
Exercise: Clustering neighborhoods by Airbnb stats¶
I've extracted neighborhood Airbnb statistics for Philadelphia neighborhoods from Tom Slee's website.
The data includes average price per person, overall satisfaction, and number of listings.
Two good references for Airbnb data¶
Original research study: How Airbnb's Data Hid the Facts in New York City
Step 1: Load the data with pandas¶
The data is available in CSV format ("philly_airbnb_by_neighborhoods.csv") in the "data/" folder of the repository.
In [5]:
airbnb = pd.read_csv("data/philly_airbnb_by_neighborhoods.csv")
airbnb.head()
Out[5]:
Step 2: Perform the K-Means fit¶
- Use our three features:
price_per_person,overall_satisfaction,N - I used 5 clusters, but you are welcome to experiment with different values!
- Scaling the features is recommended, but if the scales aren't too different, so probably isn't necessary in this case
In [6]:
# Initialize the Kmeans object
kmeans = KMeans(n_clusters=5, random_state=42)
# Scale the data features we want
scaler = StandardScaler()
scaled_airbnb_data = scaler.fit_transform(airbnb[['price_per_person', 'overall_satisfaction', 'N']])
In [7]:
# Run the fit!
kmeans.fit(scaled_airbnb_data)
# Save the cluster labels
airbnb['label'] = kmeans.labels_
In [8]:
# New column "label"!
airbnb.head()
Out[8]:
Step 3: Calculate average features per cluster¶
To gain some insight into our clusters, after calculating the K-Means labels:
- group by the
labelcolumn - calculate the
mean()of each of our features - calculate the number of neighborhoods per cluster
In [9]:
airbnb.groupby('label', as_index=False).size()
Out[9]:
In [10]:
airbnb.groupby('label', as_index=False).mean().sort_values(by='price_per_person')
Out[10]:
In [11]:
airbnb.loc[airbnb['label'] == 2]
Out[11]:
In [12]:
airbnb.loc[airbnb['label'] == 3]
Out[12]:
In [13]:
airbnb.loc[airbnb['label'] == 4]
Out[13]:
Step 4: Plot a choropleth, coloring neighborhoods by their cluster label¶
- Part 1: Load the Philadelphia neighborhoods available in the data directory:
./data/philly_neighborhoods.geojson
- Part 2: Merge the Airbnb data (with labels) and the neighborhood polygons
- Part 3: Use geopandas to plot the neighborhoods
- The
categorical=Trueandlegend=Truekeywords will be useful here
- The
In [14]:
hoods = gpd.read_file("./data/philly_neighborhoods.geojson")
hoods.head()
Out[14]:
In [15]:
airbnb.head()
Out[15]:
In [16]:
# do the merge
airbnb2 = hoods.merge(airbnb, left_on='Name', right_on='neighborhood', how='left')
# assign -1 to the neighborhoods without any listings
airbnb2['label'] = airbnb2['label'].fillna(-1)
In [17]:
# plot the data
airbnb2 = airbnb2.to_crs(epsg=3857)
# setup the figure
f, ax = plt.subplots(figsize=(10, 8))
# plot, coloring by label column
# specify categorical data and add legend
airbnb2.plot(
column="label",
cmap="Dark2",
categorical=True,
legend=True,
edgecolor="k",
lw=0.5,
ax=ax,
)
ax.set_axis_off()
plt.axis("equal");
Step 5: Plot an interactive map¶
Use altair to plot the clustering results with a tooltip for neighborhood name and tooltip.
Hint: See week 3B's lecture on interactive choropleth's with altair
In [18]:
# plot map, where variables ares nested within `properties`,
alt.Chart(airbnb2.to_crs(epsg=4326)).mark_geoshape().properties(
width=500, height=400,
).encode(
tooltip=["Name:N", "label:N"],
color=alt.Color("label:N", scale=alt.Scale(scheme="Dark2")),
)
Out[18]:
Based on these results, where would you want to stay?¶
Group 2!
Why 5 groups?¶
Use the "Elbow" method)...
In [19]:
# Number of clusters to try out
n_clusters = list(range(2, 10))
# Run kmeans for each value of k
inertias = []
for k in n_clusters:
# Initialize and run
kmeans = KMeans(n_clusters=k)
kmeans.fit(scaled_airbnb_data)
# Save the "inertia"
inertias.append(kmeans.inertia_)
# Plot it!
plt.plot(n_clusters, inertias, marker='o', ms=10, mfc='white', lw=4, mew=3)
Out[19]:
The kneed package¶
There is also a nice package called kneed that can determine the “knee” point quantitatively, using the kneedle algorithm.
In [20]:
from kneed import KneeLocator
# Initialize the knee algorithm
kn = KneeLocator(n_clusters, inertias, curve='convex', direction='decreasing')
# Print out the knee
print(kn.knee)
Part 2: Spatial clustering¶
Now on to the more traditional view of "clustering"...
DBSCAN¶
"Density-Based Spatial Clustering of Applications with Noise"
- Clusters are areas of high density separated by areas of low density.
- Can identify clusters of any shape
- Good at separating core samples in high-density regions from low-density noise samples
- Best for spatial data
Two key parameters¶
- eps: The maximum distance between two samples for them to be considered as in the same neighborhood.
- min_samples: The number of samples in a neighborhood for a point to be considered as a core point (including the point itself).
Example Scenario¶
Example Scenario¶
min_samples= 4.- Point A and the other red points are core points
- There are at least
min_samples(4) points (including the point itself) within a distance ofepsfrom each of these points. - These points are all reachable from one another, so they for a cluster
- There are at least
- Points B and C are edge points of the cluster
- They are reachable (within a distance of
eps) from the core points so they are part of the cluster - They do not have at least
min_ptswithin a distance ofepsso they are not core points
- They are reachable (within a distance of
- Point N is a noise point — it not within a distance of
epsfrom any of the cluster points
Importance of parameter choices¶
Higher min_samples or a lower eps requires a higher density necessary to form a cluster.
Example: OpenStreetMap GPS traces in Philadelphia¶
- Data extracted from the set of 1 billion GPS traces from OSM.
- CRS is EPSG=3857 —
xandyin units of meters
In [22]:
coords = gpd.read_file('./data/osm_gps_philadelphia.geojson')
coords.head()
Out[22]:
In [23]:
num_points = len(coords)
print(f"Total number of points = {num_points}")
DBSCAN basics¶
In [24]:
from sklearn.cluster import dbscan
In [25]:
dbscan?
In [32]:
# some parameters to start with
eps = 50 # in meters
min_samples = 50
cores, labels = dbscan(coords[["x", "y"]], eps=eps, min_samples=min_samples)
The function returns two objects, which we call cores and labels. cores contains the indices of each point which is classified as a core.
In [33]:
# The first 5 elements
cores[:5]
Out[33]:
The length of cores tells you how many core samples we have:
In [34]:
num_cores = len(cores)
print(f"Number of core samples = {num_cores}")
The labels tells you the cluster number each point belongs to. Those points classified as noise receive a cluster number of -1:
In [35]:
# The first 5 elements
labels[:5]
Out[35]:
The labels array is the same length as our input data, so we can add it as a column in our original data frame
In [36]:
# Add our labels to the original data
coords['label'] = labels
In [37]:
coords.head()
Out[37]:
The number of clusters is the number of unique labels minus one (because noise has a label of -1)
In [38]:
num_clusters = coords['label'].nunique() - 1
print(f"number of clusters = {num_clusters}")
We can group by the label column to get the size of each cluster:
In [39]:
cluster_sizes = coords.groupby('label', as_index=False).size()
cluster_sizes
Out[39]:
In [40]:
# All points get assigned a cluster label (-1 reserved for noise)
cluster_sizes['size'].sum() == num_points
Out[40]:
The number of noise points is the size of the cluster with label "-1":
In [41]:
num_noise = cluster_sizes.iloc[0]['size']
print(f"number of noise points = {num_noise}")
If points aren't noise or core samples, they must be edges:
In [42]:
num_edges = num_points - num_cores - num_noise
print(f"Number of edge points = {num_edges}")
Now let's plot the noise and clusters¶
- Extract each cluster: select points with the same label number
- Plot the cluster centers: the mean
xand meanyvalue for each cluster
In [43]:
# Setup figure and axis
f, ax = plt.subplots(figsize=(10, 10), facecolor="black")
# Plot the noise samples in grey
noise = coords.loc[coords["label"] == -1]
ax.scatter(noise["x"], noise["y"], c="grey", s=5, linewidth=0)
# Loop over each cluster number
for label_num in range(0, num_clusters):
# Extract the samples with this label number
this_cluster = coords.loc[coords["label"] == label_num]
# Calculate the mean (x,y) point for this cluster in red
x_mean = this_cluster["x"].mean()
y_mean = this_cluster["y"].mean()
# Plot this centroid point in red
ax.scatter(x_mean, y_mean, linewidth=0, color="red")
# Format
ax.set_axis_off()
ax.set_aspect("equal")
Extending DBSCAN beyond just spatial coordinates¶
DBSCAN can perform high-density clusters from more than just spatial coordinates, as long as they are properly normalized
Exercise: Extracting patterns from NYC taxi rides¶
I've extracted data for taxi pickups or drop offs occurring in the Williamsburg neighborhood of NYC from the NYC taxi open data.
Includes data for:
- Pickup/dropoff location
- Fare amount
- Trip distance
- Pickup/dropoff hour
Goal: identify clusters of similar taxi rides that are not only clustered spatially, but also clustered for features like hour of day and trip distance
Inspired by this CARTO blog post
Step 1: Load the data¶
In [44]:
taxi = pd.read_csv("./data/williamsburg_taxi_trips.csv")
taxi.head()
Out[44]:
Step 2: Extract and normalize several features¶
We will focus on the following columns:
pickup_xandpickup_ydropoff_xanddropoff_ytrip_distancepickup_hour
Use the StandardScaler to normalize these features.
In [45]:
feature_columns = [
"pickup_x",
"pickup_y",
"dropoff_x",
"dropoff_y",
"trip_distance",
"pickup_hour",
]
features = taxi[feature_columns].copy()
In [46]:
features.head()
Out[46]:
In [47]:
# Scale these features
scaler = StandardScaler()
scaled_features = scaler.fit_transform(features)
scaled_features
Out[47]:
In [48]:
print(scaled_features.shape)
print(features.shape)
Step 3: Run DBSCAN to extract high-density clusters¶
- We want the highest density clusters, ideally no more than about 30-50 clusters.
Run the DBSCAN and experiment with different values of
epsandmin_samples- I started with
epsof 0.25 andmin_samplesof 50
- I started with
Add the labels to the original data frame and calculate the number of clusters. It should be less than 50 or so.
Hint: If the algorithm is taking a long time to run (more than a few minutes), the eps is probably too big!
In [49]:
# Run DBSCAN
cores, labels = dbscan(scaled_features, eps=0.25, min_samples=50)
# Add the labels back to the original (unscaled) dataset
features['label'] = labels
In [50]:
# Extract the number of clusters
num_clusters = features['label'].nunique() - 1
print(num_clusters)
In [51]:
features
Out[51]:
Step 4: Identify the 5 largest clusters¶
Group by the label, calculate and sort the sizes to find the label numbers of the top 5 largest clusters
In [52]:
N = features.groupby('label').size()
print(N)
In [53]:
# sort from largest to smallest
N = N.sort_values(ascending=False)
N
Out[53]:
In [59]:
# extract labels (ignoring label -1 for noise)
top5 = N.iloc[1:6]
top5
Out[59]:
In [60]:
top5_labels = top5.index.tolist()
top5_labels
Out[60]:
Step 5: Get mean statistics for the top 5 largest clusters¶
To better identify trends in the top 5 clusters, calculate the mean trip distance and pickup_hour for each of the clusters.
In [61]:
# get the features for the top 5 labels
selection = features['label'].isin(top5_labels)
# select top 5 and groupby by the label
grps = features.loc[selection].groupby('label')
# calculate average pickup hour and trip distance per cluster
avg_values = grps[['pickup_hour', 'trip_distance']].mean()
avg_values
Out[61]:
Step 6a: Visualize the top 5 largest clusters¶
Now visualize the top 5 largest clusters:
- plot the dropoffs and pickups (same color) for the 5 largest clusters
- include the "noise" samples, shown in gray
Hints:
- For a given cluster, plot the dropoffs and pickups with the same color so we can visualize patterns in the taxi trips
- A good color scheme for a black background is given below
In [62]:
# a good color scheme for a black background
colors = ['aqua', 'lime', 'red', 'fuchsia', 'yellow']
In [63]:
# EXAMPLE: enumerating a list
for i, label_num in enumerate([0, 1, 2, 3, 4]):
print(f"i = {i}")
print(f"label_num = {label_num}")
In [64]:
# Setup figure and axis
f, ax = plt.subplots(figsize=(10, 10), facecolor="black")
# Plot noise in grey
noise = features.loc[features["label"] == -1]
ax.scatter(noise["pickup_x"], noise["pickup_y"], c="grey", s=5, linewidth=0)
# specify colors for each of the top 5 clusters
colors = ["aqua", "lime", "red", "fuchsia", "yellow"]
# loop over top 5 largest clusters
for i, label_num in enumerate(top5_labels):
print(f"Plotting cluster #{label_num}...")
# select all the samples with label equals "label_num"
this_cluster = features.loc[features["label"] == label_num]
# plot pickups
ax.scatter(
this_cluster["pickup_x"],
this_cluster["pickup_y"],
linewidth=0,
color=colors[i],
s=5,
alpha=0.3,
)
# plot dropoffs
ax.scatter(
this_cluster["dropoff_x"],
this_cluster["dropoff_y"],
linewidth=0,
color=colors[i],
s=5,
alpha=0.3,
)
# Display the figure
ax.set_axis_off()
If you're feeling ambitious, and time-permitting...¶
Step 6b: Visualizing one cluster at a time¶
Another good way to visualize the results is to explore the other clusters one at a time, plotting both the pickups and dropoffs to identify the trends.
Use different colors for pickups/dropoffs to easily identify them.
Make it a function so we can repeat it easily:
In [ ]:
def plot_taxi_cluster(label_num):
"""
Plot the pickups and dropoffs for the input cluster label
"""
# Setup figure and axis
f, ax = plt.subplots(figsize=(10, 10), facecolor='black')
# Plot noise in grey
noise = features.loc[features['label']==-1]
ax.scatter(noise['pickup_x'], noise['pickup_y'], c='grey', s=5, linewidth=0)
# get this cluster
this_cluster = features.loc[features['label']==label_num]
# Plot pickups in fuchsia
ax.scatter(this_cluster['pickup_x'], this_cluster['pickup_y'],
linewidth=0, color='fuchsia', s=5, alpha=0.3)
# Plot dropoffs in aqua
ax.scatter(this_cluster['dropoff_x'], this_cluster['dropoff_y'],
linewidth=0, color='aqua', s=5, alpha=0.3)
# Display the figure
ax.set_axis_off()
plt.show()
In [ ]:
# plot specific label
plot_taxi_cluster(label_num=6)
Step 7: an interactive map of clusters with hvplot + datashader¶
- We'll plot the pickup/dropoff locations for the top 5 clusters
- Use the
.hvplot.scatter()function to plot the x/y points - Specify the
c=keyword as the column holding the cluster label - Use the below
color_keyas thecmap=keyword - Specify the aggregator as
ds.count_cat("label")— this will color the data points using our categorical color map - Use the
datashade=Truekeyword to tell hvplot to use datashader - Add a background tile using Geoviews
- Combine the pickups, dropoffs, and background tile into a single interactive map
In [65]:
# map colors to the top 5 cluster labels
color_key = dict(zip(top5_labels, ['aqua', 'lime', 'red', 'fuchsia', 'yellow']))
print(color_key)
In [66]:
# extract the features for the top 5 clusters
top5_features = features.loc[features['label'].isin(top5_labels)]
In [67]:
import hvplot.pandas
import datashader as ds
import geoviews as gv
In [68]:
# Pickups using a categorical color map
pickups = top5_features.hvplot.scatter(
x="pickup_x",
y="pickup_y",
c="label",
aggregator=ds.count_cat("label"),
datashade=True,
cmap=color_key,
width=800,
height=600,
)
# Dropoffs using a categorical color map
dropoffs = top5_features.hvplot.scatter(
x="dropoff_x",
y="dropoff_y",
c="label",
aggregator=ds.count_cat("label"),
datashade=True,
cmap=color_key,
width=800,
height=600,
)
pickups * dropoffs * gv.tile_sources.CartoDark
Out[68]:
In [ ]: