This page was generated from /home/runner/work/esda/esda/docs/notebooks/adbscan_berlin_example.ipynb. Interactive online version: Binder badge

Cluster points and explore boundary blurriness with A-DBSCAN

In this example, we will illustrate how to use A-DBSCAN (Arribas-Bel et al., 2019) with a sample of AirBnb properties in Berlin. A-DBSCAN will allow us do two things:

  • Identify clusters of high density of AirBnb properties and delineate their boundaries

  • Explore the stability of such boundaries

%matplotlib inline

import pandas
import geopandas
import numpy as np
import contextily as cx
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from import alpha_shape_auto

import sys
    from esda.adbscan import ADBSCAN, get_cluster_boundary, remap_lbls
# This below can be removed once A-DBSCAN is merged into `esda`
    print("Import from local folder...")
    import sys
    from adbscan import ADBSCAN, get_cluster_boundary, remap_lbls


We will be using the Berlin extract from Inside Airbnb. This is the same dataset used in the Scipy 2018 tutorial on Geospatial data analysis with Python.

tab = pandas.read_csv("data/berlin-listings.csv")
Unnamed: 0 id listing_url scrape_id last_scraped name summary space description experiences_offered ... review_scores_value requires_license license jurisdiction_names instant_bookable cancellation_policy require_guest_profile_picture require_guest_phone_verification calculated_host_listings_count reviews_per_month
0 0 17260587 20170507222235 2017-05-08 Kunterbuntes Zimmer mit eigenem Bad für jedermann Meine Unterkunft ist gut für paare, alleinreis... NaN Meine Unterkunft ist gut für paare, alleinreis... none ... 10.0 f NaN NaN t flexible f f 3 2.00
1 1 17227881 20170507222235 2017-05-08 Modernes Zimmer in Berlin Pankow Es ist ein schönes gepflegtes und modernes Zim... Das Haus befindet sich direkt vor eine Tram Ha... Es ist ein schönes gepflegtes und modernes Zim... none ... 10.0 f NaN NaN t flexible f f 1 1.29

2 rows × 96 columns

The original dataset includes more than 20,000 observations:

(20053, 96)

To make the illustration run a bit speedier on any hardware, let us pick a random sample of 10% of the set; that is we’ll pick 2,000 properties at random:

tab = tab.sample(n=2000, random_state=1234)

For convenience, we convert into a GeoDataFrame where the geometries are built from the lon/lat columns in the original table:

db_ll = geopandas.GeoDataFrame(tab,
                               geometry=geopandas.points_from_xy(tab.longitude, tab.latitude),
                               crs={'init': 'epsg:4326'}
/home/serge/anaconda3/envs/analytical/lib/python3.7/site-packages/pyproj/crs/ FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes:
  return _prepare_from_string(" ".join(pjargs))

Because we are going to run an algorithm that relies on distances, we need to be able to calculate such distances in a projected plane. Instead of the original lon/lat coordinates, we use the `ETRS89 <>`__ projection, designed for Germany and expressed in metres:

db = db_ll.to_crs(epsg=5243)

Voila! We can now visualise our dataset:

ax = db.plot(markersize=0.1, color='orange')

And before we get going, we will also pull out the projected coordinates into separate columns:

db["X"] = db.geometry.x
db["Y"] = db.geometry.y

Identify clusters of AirBnb

A-DBSCAN, similar to the original DBSCAN algorithm, requires two parameters before it can be run:

  1. eps: maximum distance to look for neighbors from each location

  2. min_samples: minimum number of neighboring points required for a point to be considered part of a cluster

For this example, we will pick a 1% of the overall sample size for the min_samples parameter:

db.shape[0] * 0.01

We will use a maximum radious of 500m for the eps parameter. This implicitly implies we are looking at a density of at least about 25 properties per Sq. \(Km^2\) (\(dens = \frac{20}{\pi \, \times \, 0.5^2}\)).

Once we know the parameters we’ll use, we can go ahead and run A-DBSCAN:

# Get clusters
adbs = ADBSCAN(500, 20, pct_exact=0.5, reps=10, keep_solus=True)
CPU times: user 755 ms, sys: 3.36 ms, total: 758 ms
Wall time: 752 ms
ADBSCAN(eps=500, keep_solus=True, min_samples=20, pct_exact=0.5, reps=10)

Once fit(), we can extract the labels in a similar way that we would from scikit-learn, and then we can plot on a map:

ax = db.assign(lbls=adbs.votes["lbls"])\
             figsize=(12, 12)

This only displays each property colored based on its label assigned. Beyond this, we can also create polygons that represent a tight boundary around all the points in a given cluster. To do this, we use the \(\alpha\)-shapes algorithm through the helper function get_cluster_boundary:

%time polys = get_cluster_boundary(adbs.votes["lbls"], db,
CPU times: user 1.76 s, sys: 15.2 ms, total: 1.78 s
Wall time: 1.77 s

These polygons can also be plotted as any other GeoSeries object:

ax = polys.plot(alpha=0.5, color="red")

And, just for fun, we can create more zoomed-in views, dimming out areas outside each cluster:

f, axs = plt.subplots(1, 3, figsize=(18, 6))
for i, ax in enumerate(axs):
    # Plot the boundary of the cluster found
    ax = polys[[i]].plot(ax=ax,
    # Add basemap
    # Extra to dim non-cluster areas
    (minX, maxX), (minY, maxY) = ax.get_xlim(), ax.get_ylim()
    bb = Polygon([(minX, minY),
                  (maxX, minY),
                  (maxX, maxY),
                  (minX, maxY),
                  (minX, minY)
    ax.set_title(f"Cluster {polys[[i]].index[0]}")
/home/serge/anaconda3/envs/analytical/lib/python3.7/site-packages/ FutureWarning: The "url" option is deprecated. Please use the "source" argument instead.
  if sys.path[0] == '':
/home/serge/anaconda3/envs/analytical/lib/python3.7/site-packages/ FutureWarning: The "url" option is deprecated. Please use the "source" argument instead.
  if sys.path[0] == '':
/home/serge/anaconda3/envs/analytical/lib/python3.7/site-packages/ FutureWarning: The "url" option is deprecated. Please use the "source" argument instead.
  if sys.path[0] == '':

Explore boundary blurriness

One of the key benefits of A-DBSCAN is that, because it relies on an ensemble approach that generates several candidate solutions, it allows us to explore to what extent boundaries are stable and clearcut or blurrier. To do this, we need to extract from the ADBSCAN object each candidate solution and turn them into boundary lines (this may take a little bit to run):

solus_rl = remap_lbls(adbs.solus, db, n_jobs=-1)
lines = []
for rep in solus_rl:
    line = get_cluster_boundary(solus_rl[rep],
    line = line.boundary
    line = line.reset_index()\
               .rename(columns={0: "geometry",
                                "index": "cluster_id"}
lines = pandas.concat(lines)
lines = geopandas.GeoDataFrame(lines,
CPU times: user 665 ms, sys: 1.45 s, total: 2.12 s
Wall time: 4.83 s

Here is a first look at all the solutions drawn up from the simulation:

ax = lines.plot(color="#FFDB58",
/home/serge/anaconda3/envs/analytical/lib/python3.7/site-packages/ FutureWarning: The "url" option is deprecated. Please use the "source" argument instead.
  import sys

Having every single candidate boundary stored and labelled into a table, it allows us to make several queries. For example, here are all the solutions generated in the first replication:

lines.query("rep == 'rep-00'").plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fc0003522d0>

And here we have all the candidate solutions for cluster label No. 2:

lines.query("cluster_id == '2'").plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fc000452e90>

Finally, we can take this idea into an interactive context by building widgets that allow us to brush through replications:

from ipywidgets import interact, IntSlider
def plot_rep(rep):
    f, ax = plt.subplots(1, figsize=(9, 9))
    # Background points
    db[["X", "Y"]].plot.scatter("X", "Y", ax=ax, color="0.25", s=0.5)
    # Boundaries
    cs = lines.query(f"rep == 'rep-{str(rep).zfill(2)}'")
    cs.plot(ax=ax, color="red")
    # Cluster IDs
    for s, row in cs.iterrows():
    return None
reps = range(len(lines["rep"].unique()))
slider = IntSlider(min=min(reps), max=max(reps), step=1)
interact(plot_rep, rep=slider);

Or to pull out all solutions for a given cluster IDs:

def plot_cluster(cluster_id):
    f, ax = plt.subplots(1, figsize=(9, 9))
    # Background points
    db[["X", "Y"]].plot.scatter("X",
    # Boundaries
    lines.query(f"cluster_id == '{cluster_id}'").plot(ax=ax,
    return None
interact(plot_cluster, cluster_id=lines["cluster_id"].unique());