Detection of hierarchical crowd activity structures in geographic point data

Natural cities

The Natural Cities algorithm proposed in Jiang & Miao (2015) looks to objectively define and delineate human settlements or activities at different scales from large amounts of geographic information. The algorithm exploits the heavy tailed distribution of sizes present in geographic data and uses Head-Tails breaks (Jiang, 2013) to iteratively extract the structures at successive scales. The algorithm can be summarized as follows:

• Use all the points and generate a Triangulated Irregular Network (TIN) using the Delaunay triangulation algorithm.

• Extract the edges of the triangulation and obtain their lengths.

• Calculate the average length.

• Use the average as a threshold to split the edges in two categories, the Head with those segments larger than the mean, and the Tail, with the segments smaller than the mean.

• Remove the Head.

• Generate the least amount of continuous polygons from the union of Tail edges.

• The points that are inside the obtained polygons are kept and the ones outside are consider noise.

• Repeat the procedure iteratively for each resulting polygon until the size distribution does not resemble a heavy tailed distribution.

At the end of the procedure, a hierarchical structure is obtained that contains structures within structures. This algorithms has been used to examine the spatial structure of road networks and their relation to social media messages (Jiang & Ren, 2019) and the evolution of Natural Cities from geo-tagged social media (Jiang & Miao, 2015).

DBSCAN

DBSCAN is an algorithm designed to detect arbitrarily shaped density clusters in databases with noise (Ester et al., 1996). In the following we provide a brief description of the algorithm along with its core assumptions. The definitions established here will serve as basis to the description of OPTICS (‘OPTICS’), HDBSCAN (‘HDBSCAN’) and Adaptative DBSCAN (‘Adaptative DBSCAN’).

DBSCAN uses a threshold distance (epsilon distance ɛ) and a minimum number of points minPts as initial parameters. From this, DBSCAN makes the following basic definitions:

Definition (Neighborhood). Let p ∈ P be a point in the data set P and ɛ > 0 a threshold distance, then the epsilon-neighborhood of p is ${\displaystyle {\mathcal{N}}_{\varepsilon }\left(p\right)=\left\{x|d\left(x,p\right)<\varepsilon \right\}}$ for x ∈ P

Definition (Reachable). A point q ∈ P is a reachable point from a point p ∈ P if there is a path of points {p₁, …p_n − 1, q} ⊂ P such that the distance between consecutive points is less than ɛ (d(p_i, p_i+1) < ɛ).

Definition (Core). Let P be a set of points, a point p is a core point of P with respect to ɛ and minPts if $\left|{\mathcal{N}}_{\varepsilon }\left(p\right)\right|\ge minPts$ where ${\mathcal{N}}_{\varepsilon }\subset P$.

Definition (Density reachable). A point p is density reachable from a point q with respect of ɛ and minPts within a set of points P if there is a path {p₁, …p_n − 1, q} ⊂ P such that p_i ∈ P are core points of P and p_i+1 ∈ N_ɛ(p_i).

Definition (Density connected). A point p is density connected to a point q with respect to ɛ and minPts if there is a point o such that both, p and q are density reachable from o with respect to ɛ and minPts.

Definition (Cluster). A cluster C ⊂ P with respect to ɛ and minPts is a non empty subset of P that satisfies the following conditions:

• ∀p, q ∈ P: if p ∈ C and q is density connected from p with respect to ɛ and minPts, then q ∈ C.

• ∀p, q ∈ C p is density connected to q with respect to ɛ and minPts.

Definition (Noise). Let C₁, …, C_k be the clusters of a set of points P with respect to ɛ and minPts. The noise is defined as the set of points in P that do not belong to any cluster C_i.

Using the above definitions we can summarize the algorithm to cluster a set of points P as:

• Consider the set P_C = {p₁, …p_m} ⊂ P of core points of P with respect to ɛ and minPts.

• For a point p ∈ P_C take all the density connected points as part of the same cluster, remove these points from P and P_C.

• Repeat the step above until P_C is empty.

• If a point can not be reached from any core point is considered as noise.

In order to fit our aim of detecting structures within structures, we apply DBSCAN recursively to each of the discovered clusters. It is important to consider that the ɛ parameter in DBSCAN is related to the relative point density of noise with respect to clusters, this means that in each recursive application one needs to set an appropriate value for ɛ; preserving the same value across iterations would not uncover density structures within clusters since those structures would, by definition, be of higher density than the encompassing cluster. In López-Ramírez et al. (2018), ɛ is decreased by a constant arbitrary factor on each iteration, in ‘Adaptative DBSCAN’ we propose a novel method for automatically selecting ɛ for each cluster and iteration that overcomes this limitation.

OPTICS

The OPTICS algorithm is fully described in Ankerst et al. (1999), here we will make a brief presentation of its main characteristics, adapted to our research problem. OPTICS works in a similar fashion to DBSCAN and can be seen as an extension of it. The algorithm uses the notion of density-based cluster ordering to extract the corresponding density-based cluster for each point. OPTICS extends the definitions used in DBSCAN by adding:

Definition (Core distance). The core distance (d_core(p, ɛ, minPts)) of a point p with respect to ɛ and minPts, is the smallest distance ɛ′ < ɛ such that p is considered a core point with respect to minPts and ɛ′. If $\left|{\mathcal{N}}_{\varepsilon }\left(p\right)\right|<minPts$ then the core distance is undefined.

Definition (Reachability distance). Let p and o be points in a set P and take ${\mathcal{N}}_{\varepsilon }\left(o\right)$. The reachability-distance (reach − dist_(ɛ,minpts)(p, o)) is undefined if $\left|{\mathcal{N}}_{\varepsilon }\left(o\right)\right|<minPts$, else is max(d_core(p, ɛ, minPts), d(p, o)).

Let P be a set of points, and set values for ɛ and minPts. The OPTICS algorithm can be summarized as following:

• For each point p ∈ P set the reachability-distance value as UNDEFINED (L_reach−dis(p) = UNDEFINED ∀p ∈ D).

• Set a empty processed list L_Pro.

• Set a empty priority queue S

• For each unprocessed point p in P do:

  • Obtain ${\mathcal{N}}_{\varepsilon }\left(p\right)$

  • Mark p as processed (add p to L_Pro)

  • Push p priority queue S

  • If d_core(p, ɛ, minPts) is UNDEFINED move to the next unprocessed point. If not then:

    • Update the queue S based on the rechability-distance by the update function and use ${\mathcal{N}}_{\varepsilon }\left(p\right)$, p, ɛ, and minPts as input parameters.

    • For each q in S do:

      • Obtain ${\mathcal{N}}_{\varepsilon }\left(q\right)$

      • Mark q as processed (add q to L_Pro).

      • Push q in the priority queue S.

      • If the d_core(p, ɛ, minPts) is not UNDEFINED, use the update function with the ${\mathcal{N}}_{\varepsilon }\left(q\right)$, q, ɛ, and minPts as parameters to update S.

• The algorithm expands S until no points can be added.

The update function for a priority queue S, uses as parameters a neighborhood ${\mathcal{N}}_{\varepsilon }\left(p\right)$, the center of the neighborhood p, and the ɛ and minPts values, and is defined as:

• Get the core-distance of p (d_core(p, ɛ, minPts)).

• For all $o\in {\mathcal{N}}_{\varepsilon }\left(p\right)$ if o is not processed, then

  • Define a new reachability-distance as new_reach−dis = max(d_core(p, ɛ, minPts), d(p, o))

  • If L_reach−dis(o) is UNDEFINED (is not in S) then L_reach−dis(o) = new_reach−dis and insert o in S with value new_reach−dis.

  • If L_reach−dis(o) is not UNDEFINED (o is already in S), then if new_reach−dis < L_reach−dis(o) update the position of o in S moving forward with the new value new_reach−dis.

After this, we have reachability-distance values (L_reach−dis:D → ℝ∪UNDEFINED) for every point in P and an ordered queue S with respect to ɛ and minPts. Using this queue, the clusters are extracted using ɛ′ < ɛ distance and minPts by assigning the cluster membership depending on the reachability-distance.

To decide whether a given point is noise or the first computed element in a cluster, is necessary to define that UNDEFINED > ɛ > ɛ′. After this, the assignation of the elements of S is performed as:

• Set the Cluster_ID = NOISE.

• Consider the points p ∈ S.

• If L_reach−dis(p) > ɛ′ then,

  • If d_core(p, ɛ′, minPts) ≤ ɛ′, the value of Cluster_ID is set to the next value, and the cluster for the element p is set as Cluster_ID.

  • If d_core(p, ɛ′, minPts) ≥ ɛ′, the point p is consider as NOISE.

• If L_reach−dis(p) ≤ ɛ′, the cluster for the element p is Cluster_ID.

The order of S guarantees that all the elements in the same cluster are close. The L_reach−dis values allow us to distinguish between clusters, these clusters will depend on the ɛ′ parameter.

In our work we use a recursive application of OPTICS to uncover the hierarchical scale structure, in the same fashion as López-Ramírez et al. (2018).

HDBSCAN

The algorithm Density Based Clustering based on Hierarchical Density Estimates is presented in Campello, Moulavi & Sander (2013). The main intuition behind HDBSCAN is that the most significant clusters are those that are preserved along the hierarchical distribution. Before developing an explanation of HDBSCAN we need the following definitions:

Definition (Core k-distance). Let p be a point in space, the core distance d_corek(p) of k is defined by d_corek(p) = max({d(x, p)|x ∈ k-NN(p) if x ≠ p}) where k-NN(p) is the set of the k-nearest neighbors of p for a specific k ∈ ℕ.

Definition (Core k-distance point). A point p is a Core k-distance point with respect to ɛ and minPts, if the number of elements in the core neighborhood of p (${\mathcal{N}}_{cor{e}_k}\left(p,\varepsilon ,minPts\right)=\left\{x|d_{core}\left(x,p\right)<\varepsilon \right\}$) is greater or equal than minPts. Where ɛ is greater or equal to the core distance of p for a given k value.

Definition (Mutual Reachability Distance). For every p and q the Mutual Reachability Distance is d_mreach(p, q) = max({d_corek(p), d_corek(q), d(p, q)})

Consider the weighted graph with every point as vertices and the mutual reachability distance as weights. The resulting graph will be strongly connected, the idea is to find islands within the graph that will be consider clusters.

To classify the vertices of the graph in different islands its easier to determine which do not belong to the same island. This is done by removing edges with larger weights than a threshold value ɛ, by reducing the threshold the graph will start to disconnect and the connected components will be the islands (clusters) to consider. If a vertex in the graph is disconnected from the graph or the connected component doesn’t have the minPts, the corresponding vertex are considered noise. The result will depend on ɛ.

To reduce computing time HDBSCAN finds a Minimum spanning tree (MST) of the complete graph. The tree is built one edge at a time by adding the edge with the lowest weight that connects the current tree to a vertex not yet in the tree.

Using the MST and adding a self edge to all the vertices in the graph with a distance core value(MST_ext), a dendrogram with a hierarchy is extracted. The HDBSCAN hierarchy is extracted using the following rules (Campello, Moulavi & Sander, 2013):

• For the root of the tree assign all objects the same label (single cluster)

• Iteratively remove all edges from MST_ext in decreasing order of weights (in case of ties, edges must be removed simultaneously):

  • Before each removal, set the dendrogram scale value of the current hierarchical level as the weight of the edge(s) to be removed.

  • After each removal, assign labels to the connected component(s) that contain(s) the end vertex(-ices) of the removed edge(s), to obtain the next hierarchical level: assign a new cluster label to a component if it still has at least one edge, else assign it a null label (noise).

The clusters thus obtained from the dendrogram depend on the selection of a λ parameter based on the estimation of the stability of a cluster. This is done using a probability density function approximated by the k-nearest neighbors.

Once again, instead of relying on the hierarchical structure obtained from HDBSCAN, to discover our structures within structures we apply the algorithm recursively to every cluster.

Adaptative DBSCAN

In López-Ramírez et al. (2018) the authors propose the recursive application of DBSCAN as a suitable algorithm to uncover the hierarchical structure present in spatio-temporal events data sets. The main drawbacks of this proposal are, on the one hand, the need to select appropriate ɛ and minPts values for each recursive application of DBSCAN and, on the other hand, that once this values are selected they are used for every cluster in a given hierarchical level, thus assuming that every cluster has the same intrinsic density properties. In this paper we propose an algorithm to automatically select ɛ values for each cluster, this algorithm draws on methods proposed in the available literature and adapts them to the problem of identifying clusters on hierarchically structured geographical data.

Before explaining our algorithm, let us review how minPts and ɛ values are estimated. For the minPts parameter, the general rule is to select it so minPts ≥ D + 1, where D is the dimension of the data, so in the recursive application of DBSCAN this parameter is fixed. In general, the problem of selecting ɛ is solved heuristically by the following procedure:

• Select appropriate minPts.

• For each point in the data get the kth-nearest neighbors using minPts as the value of K.

• The different distances obtained are sorted smallest to largest (Fig. 4, K-Sorted distance graph).

• Good values of ɛ distance are those where there is a big increment in the distance. This increment will correspond to an increment on the curvature in the plot, this point is called the elbow.

An automatic procedure to select ɛ is presented in Starczewski, Goetzen & Er (2020), where the authors propose a mathematical formulation to obtain the elbow value of the K-Sorted distance graph once minPts is given. This formulation tends to find the highest density clusters present in the data, so when applied recursively to hierarchically structured data it will tend to find the clusters in the deepest hierarchical level and will label as noise points belonging to intermediate hierarchical levels.

To overcome this limitation, we propose an alternative procedure to automatically select ɛ that leads to larger distance values and thus less dense clusters in each recursive application of DBSCAN. First, notice that either on the heuristic described above or in the procedure by Starczewski, Goetzen & Er (2020), the ɛ value depends on the selection of K, larger values of K will in general lead to larger ɛ values, although not in a strictly monotonic way as seen in Fig. 5. In principle it is possible to select K as the maximum number of points in the data set and obtain a suitable ɛ value, the problem is that for large data sets this is not feasible due to computational constraints.

To reduce the computational cost of finding appropriate ɛ values, we observe in Fig. 6 the similarities in the K-Sorted distance graph for a range of K values, it can be seen that the elbow value is similar for close values of K. Using this observation we obtain the K-Sorted distance graph only for a sample K_s ⊂ {1, …, N}, where N is the number of points, instead of every possible value, thus reducing the time and computational cost needed to obtain a suitable ɛ distance value.

To find the elbow value for a given k ∈ K_s, we use the library developed in Satopaa et al. (2011). The algorithm is $\mathcal{O}\left(n^2\right)$ so reducing the points involved in this step is also important. For each k, the algorithm takes as input the N∗k distances. To limit the number of distances passed to the algorithm, we average the N∗k distances in bins of size k, thus we only use N values to calculate the elbow.

Our adaptative DBSCAN algorithm can be summarized as follows:

• Let N be the number of points in the data set.

• For every step k in the range $\left\{minPts,\dots ,\lfloor \frac{N}{10}\rfloor \right\}$ obtain the following:

  • For each point p, obtain the distances to the k nearest neighbors.

  • Sort the distances.

  • Divide the distance interval in bins of size k and obtain the average for each bin.

  • Get the elbow value for all the average bins as in Satopaa et al. (2011).

• Take ɛ as the maximum of all elbow values from the previous step.

This algorithm for calculating a suitable ɛ value for each identified cluster allow us to find clusters with intermediate densities (correspondingly to intermediate hierarchical levels) and, at the same time, do it in a computationally efficient manner.

Global evaluation

To evaluate the performance of the different clustering algorithms as global classification algorithms, we label the obtained clusters using our tree cluster structure. In Fig. 7 we show an example output for the different algorithms. For every level we label points as either belonging to a cluster or as noise. This labeling has the property that when a point is labeled as noise in level n, then for every level n + j with j ≥ 1 the point will be tagged as noise. Thus we can obtain a label for every point as the concatenation of NOISE, SIGNAL tags. For example, in a three level cluster tree, the possible labels are the following: SIGNAL_SIGNAL_SIGNAL, SIGNAL_SIGNAL_NOISE, SIGNAL_NOISE_NOISE, NOISE_NOISE_NOISE.

To evaluate the algorithms, we will use the Normalized Mutual Information (NMI) between the algorithm classification and the ground truth. NMI has the advantage of being independent of the number of classes (Vinh, Epps & Bailey, 2009), thus giving a fair evaluation of the performance of the algorithms across experiments.

Definition (Normalized Mutual Information). NMI is defined as: ${\displaystyle NMI\left(X,Y\right)=\frac{2I\left(X,Y\right)}{H\left(X\right)+H\left(Y\right)}}$ where P(x) is the probability to get the label x, H(X) is the Shannon entropy of the set of labels X and I(X, Y) is the mutual information between the sets of labels X an Y.

Definition (Shannon entropy). The Shannon entropy of a set of labels X is defined as: ${\displaystyle H\left(X\right)=-\sum _{x\in X}P\left(x\right)logP\left(x\right)}$ where P(x) is the probability of label x.

Definition (Mutual information). Let X = {X₁, …, X_l} and Y = {Y₁, …Y_n} be two sets of labels for the same set of points (N), the mutual information between the two set of labels is calculated as: ${\displaystyle I\left(X,Y\right)=\sum _{i=1}^{\left|X\right|}\sum _{j=1}^{\left|Y\right|}P\left(i,j\right)\ast log\left(\frac{P\left(i,j\right)}{P\left(i\right)\ast P^{\prime }\left(j\right)}\right)}$ where P(i) = |X_i|/|N| is the probability of a random point belonging to class X_i, P′(j) = |Y_i|/|N|. And $P\left(i,j\right)=\frac{|X_i\cap Y_j|}{\left|N\right|}$ is the probability that a random point belongs to both classes X_i and Y_j.

Per level evaluation

To reflect the iterative clustering process, we evaluate the classification obtained for each level. In order to carry out this evaluation, we label, for each level, the points as either SIGNAL for points belonging to a cluster or NOISE otherwise. For each level n we include only the points that have the SIGNAL label in n − 1 level since those are the only points seen by the algorithm at each recursive iteration. Using this approach, the per level evaluation corresponds to the clusters obtained on level n and is performed only for the points that correspond to each level.

In general, the labels obtained for each level are unbalanced with more samples of the noise class. Therefore for the evaluation we use the Balanced Accuracy (BA) for binary classification (Brodersen et al., 2010) on each level.

Definition (Balanced accuracy). Let X be the ground truth labels in a point set N and Y the set of predicted labels in the same set. The balanced accuracy for a binary labeling is defined as:

${\displaystyle BA\left(X,Y\right)=\frac{1}{2}\left(\frac{TP}{TP+FN}+\frac{TN}{TN+FP}\right)}$ where TP is the true positive set of labels, TN the true negative set of labels, FN is the false negative set of labels, and FP is the false positive set of labels.

Shape evaluation

For geographic applications such as those described in ‘Detection of crowd activity scales and zones’, it is important to also evaluate the shape of the clusters obtained. Therefore we propose a measure that compares the shapes of the clusters obtained by each algorithm.

Our Similarity Shape Measure (SSM) compares the shapes of the Concave Hulls of each cluster, obtained by the optimal alpha-shape of the point set (Edelsbrunner, 1992; Bernardini & Bajaj, 1997). Thus, for our SSM, each cluster is represented by a polygon.

A simple way to compare polygon shapes is the Jaccard Index, mostly used in computer vision to compare detection algorithms.

Definition. (Jaccard Index) The Jaccard Index or Jaccard Similarity Coefficient between polygons P and Q is defined as: ${\displaystyle Jacc\left(P,Q\right)=\frac{Area\left(P\cap Q\right)}{Area\left(P\cup Q\right)}}$ where Area is the area of the polygon.

The main issue in implementing a shape similarity measure is in determining which polygons to compare in each level. As can be seen in Fig. 7, every algorithm outputs different number and shapes of polygons in each level so there is no straightforward way of assigning corresponding polygons across algorithms (or the ground truth). To overcome this issue, our Similarity Shape Measure compares, using the Jaccard Index, the polygons produced by an algorithm with the polygons in the ground truth, weighting the index by the number of points in each polygon and its corresponding intersections.

Definition (Similarity Shape Measurement). Let ${\mathcal{P}}_i^O=\left\{P_0^O,\dots P_n^O\right\}$ be the set of polygons for the clusters on the i − th level in the ground truth and ${\mathcal{Q}}_i^C=\left\{Q_0^C,\dots Q_m^C\right\}$ the polygons to evaluate for the same level, then the similarity between ${\mathcal{P}}_i^O$ and ${\mathcal{Q}}_i^C$ is:

${\displaystyle SSM\left({\mathcal{P}}_i^O,{\mathcal{Q}}_i^C\right)=\frac{\sum _{P_l^O\in {\mathcal{P}}_i^O{Q}_k^C\in Q_i^C}|P_l^O\cap Q_k^C|\ast Jacc\left(P_l^O,Q_k^C\right)}{\sum _{P\in {\mathcal{P}}_i^O}\left|P\right|+\sum _{Q\in {\mathcal{Q}}_{not}}\left|Q\right|}}$ where ${\mathcal{Q}}_{not}=\left\{Q|Q\cap P_i^O=0̸\right\}\subset {\mathcal{Q}}_i^C$ for all $P_i^O\in {\mathcal{P}}_i^O$.

The SSM takes values in [0, 1], the maximum value corresponds to the case when all polygons that intersect across the evaluation sample and the ground truth satisfy $Jacc\left(P_l^O,Q_k^C\right)=1$, all the points within the polygons $P_l^O$ are in a corresponding polygon $Q_k^C$ and the family ${\mathcal{Q}}_{not}$ is empty. Conversely, SSM takes the value 0 when no polygons intersect across the evaluated algorithm and the ground truth.

If the polygons have a large Jaccard Similarity and the number of points in the intersections is large, then SSM will also be large. On the contrary, if the polygons have low Jacacard similarity then a penalization will occur and the SSM will have a lower value even if the cardinality of the intersection is large. This means that our measure will penalize clusterizations that over or under separate points in different clusters with respect to the ground truth.

SSM will penalize for points that belong to a P_l for some l ∈ {0, …, n}, that are not inside any $Q\in {\mathcal{Q}}_i^C$. Also a penalization will occur for the points that belong to Q_k for some k ∈ {0, …, m} that not belong to any $P\in {\mathcal{P}}_i^O$. This is important to ensure that algorithms that properly classify a grater number of points as signal are not penalized.

Thus our SSM allows for the direct comparison between the polygons in the ground truth and those produced by an arbitrary algorithm, producing a global similarity measure for each level that captures not only how similar the polygons are, but also, how many points are in the most similar polygons.

Application on real-world data

As a use case of the proposed technique we will extract the crowd activity patterns for the geolocated Twitter feed for Central Mexico. The test database consists of all geolocated tweets for 2015-01-20, there are 10,366 tweets for this day. Based on the discussion on ‘Detection of crowd activity scales and zones’ and ‘Clustering Algorithms’ we will only compare the results obtained with Natural Cities, recursive DBSCAN and Adaptative-DBSCAN. In order to qualitatively compare the results of the algorithms we will focus on two tasks: first the detection of the greatest cities within the region and, second, the detection of the Central Business District (CBD) structure of Mexico City, widely described in the literature (Escamilla, Cos & Cárdenas, 2016; Suarez & Delgado, 2009). The idea behind this qualitative comparison is to understand how the different algorithms detect known crowd activity patterns at different scales.

Figure 11 through Fig. 13 show the results for the different clustering algorithms. In the larger hierarchical levels, the figures compare our results with the official delimitation of the metropolitan areas, while in the smaller scales, we show the job to housing ratio as an indicator of potential crowd activity. The first thing to notice is that, although the three algorithms are able to separate the input points into the metropolitan areas, Natural Cities tends to detect the inner cores on the first hierarchical level, while hierarchical DBSCAN and Adaptative-DBSCAN detect larger metropolitan structures. Another interesting finding is that Natural Cities detects the T-shaped pattern of Mexico City’s CBD (Escamilla, Cos & Cárdenas, 2016; Suarez & Delgado, 2009) in two iterations, while hierarchical DBSCAN and Adaptative DBSCAN require 4 and 8 iterations respectively. This means that both DBSCAN based algorithms are detecting intermediate scale structures, while Natural Cities is more aggressively reducing the size of the clusters. Another interesting quality in the clusters detected by Natural Cities in contrast to those detected by the other two algorithms is that Natural Cities tends to detect many small clusters at both levels, qualitatively these small clusters do not seem to correspond with any known structure in the city.

For the deeper levels, both in Hierarchical DBSCAN and Adaptative DBSCAN, the patterns detected seem to closely follow the job to housing ratio, this is a good qualitative indicator that the structures detected correspond with the known geographic activity patterns. To further stress this point, in Fig. 14 we show the different clusters detected with the three algorithms for afternoons (14:00 to 18:00) and evenings (18:00 to 22:00), all figures exhibit more concentrated patterns for the evenings, reflecting people gathered at their workplaces. It is also interesting that all algorithms seem to follow the T-shaped pattern of the CBD and also detect some smaller job centers to the South and West, but Natural Cities detects many smaller clusters that do not seem to correspond with relevant geographic structures.