Networks and long-range mobility in cities: A study of more than one billion taxi trips in New York City

We analyze the massive data set of more than one billion taxi trips in New York City, from January 2009 to December 2015. With these records of seven years, we generate an origin-destination matrix that has information of a vast number of trips. The mobility and flow of taxis can be described as a directed weighted network that connects different zones of high demand for taxis. This network has in and out degrees that follow a stretched exponential and a power law with an exponential cutoff distributions, respectively. Using the origin-destination matrix, we obtain a rank, called "OD rank”, analogous to the page rank of Google, that gives the more relevant places in New York City in terms of taxi trips. We introduced a model that captures the local and global dynamics that agrees with the data. Considering the taxi trips as a proxy of human mobility in cities, it might be possible that the long-range mobility found for New York City would be a general feature in other large cities around the world.


Results
Activity between zones with high demand. We explore taxi trip records taking into account the administrative boundaries including the five boroughs of New York City 36 . As a result, for the seven years studied, we have T 1 148 052 357 = taxi trips (see the Methods section for a detailed description of the datasets explored). In the following, we study this volume of data by using a grid with 500 × 500 square zones with dimensions 100 m × 100 m. Once this grid is defined, we examine the zones contained in the administrative boundaries of New York City. In Fig. 2, we present a map generated with the information of origin and destinations reported in the datasets. For each square zone defined before, we count the number of trips according to the registers of longitude and latitude of the initial and final locations of each trip for taxi registers from 2009 to 2015. The results depicted in Fig. 2 give us a first insight into the global activity of taxis. We can identify a high demand of this service in Manhattan, also the high activity in the John F. Kennedy (JFK) International Airport and how by exploring the origins of the trips we can detect some features of the street network in New York City. On the other hand, we can see in Fig. 2(b) that the destinations are less localized in specific zones observing that in the Bronx, Brooklyn and Queens boroughs the number of taxis arrivals is more uniform in comparison with the origins in Fig. 2(a). This particular feature reveals how taxi transportation manages to permeate almost all the regions of the city.
In Fig. 2, we can identify zones in New York City with low demand for taxis or where only a reduced number of taxis arrives. Even considering the counts in seven years of activity, we can identify zones with dimensions Figure 1. Schematic illustration of mobility as a spatially embedded directed weighted network. We show N = 10 square zones in the plane representing particular regions where agents can start or end a trip; we simulate T 1 000 = trips of agents between these locations. (a) Bar representation of the total number arrivals k (in) and the number of departures k (out) . (b) Diagram of the system expressed as a directed network, we represent with colors the number of trips T ij between sites i and j. In the study of human mobility, this information is expressed as an origin-destination matrix N × N with elements T ij . In particular, the directions of the links are depicted by arrows and self-loops represent the number of trips with the same origin and destination.
100 m × 100 m for which less than 10 3 taxi trips arrive or depart. This is a small number in comparison with the values of zones with a high demand for which we observe more than a million arrivals or departures. Much of these zones are located in Manhattan but also other zones of the city. In what follows, we study the flow of taxis between zones with high demand and we will describe the global spatial dynamics as a directed weighted spatial network. All the zones in our study are defined by a square with dimensions 100 m × 100 m and, for each year, we classify a region as a high activity zone if, in this specific part of the city, the number of arrivals and departures are at least 10 3 . In this way, the minimum number of arrivals at a high activity zone is at least 10 3 trips, and the same rule applies to the number of taxis leaving this region. This limit is reasonable due to the high quantity of trip records explored per year in the complete database. In addition, by using this rule we reduce possible errors produced by the inaccuracy in the origin and destination coordinates. By applying the criteria described before to the taxi trips in 2015, the number of high demand zones for this year is 4 353 and the total number of trips between these zones is = 128 984 657 T that represents a 90.22% of the original database described in the Methods section. We found similar values for the trips from 2009 to 2014. The results for the number of high demand zones N and the total number of trips T are presented in Table 1. Now, we define origin-destination matrices describing the flow of taxis between high-demand zones. In this way, the global dynamics can be explored and treated as a directed weighted network; in particular, a spatial network for which the nodes represent zones of high demand and the links with weights can represent several quantities like the flow of vehicles, the geographical distance between nodes, among other values 28 . trips from taxi trip records between January 2009 to December 2015. In this representation of the data, using only the information of origins and destinations of taxis, we can see in detail the spatial complexity of New York City and how the street network emerges from the large number of trips analyzed.  Table 1. Analysis of the spatial activity of taxi trips in New York City considering zones with a high demand for this service. By using the rule that at least 1 000 = M trips depart and arrive from a zone in a year, we obtain the number N of high demand zones. In addition, we present the total number of trips T between zones and the fraction of the original dataset that each number of trips represents. For the trips analyzed, we show the fraction of local trips with geographical distances d in the interval 0 ≤ d ≤ 1.8 Km and the fraction of long-range movements with d > 1.8 Km.
For each year, we calculate an origin-destination matrix for which the elements T ij represent the number of taxi-trips from zone i to zone j, where i, j = 1, 2, …, N denote the square zones of high demand with dimensions 100 m × 100 m. In addition to the elements of the origin-destination matrix, it is important the in-degree defined as that gives the total number of vehicles arriving at the zone i. We also have the out-degree determined by the relation i N i (out) 1 ∑ = = that counts the total number of trips originated from zone i. On the other hand, to explore the spatial activity is important to have information about the geographical distances between zones. This information is included in a N × N distance matrix D with elements d ij with the geographical distance between i and j. In addition, the degrees in Eqs. 1-2 satisfy: where T is the total number of trips considered in the origin-destination matrix. In Fig. 3, we show the origin-destination matrix and the respective matrix of distances D obtained from taxi trips in 2015. The resulting matrices incorporate the flow of vehicles between N = 4353 high demand zones.
Let us now analyze the statistical properties of the directed weighted network associated with mobility in New York City. In order to do so, we show in Fig. 4 two probability distributions: one associated to the in-degree of the network k i (in) (Fig. 4(a)) and the other one associated with the out-degree of the network k i (out) (Fig. 4(b)). We explore all the in and out-degrees, for seven years, from 2009 to 2015, in the interval 10 3 ≤ k ≤ 10 6 . With the aim of finding the best fit of the aggregated data of mobility for these distributions, we used the tools and procedures described by Clauset et al. (2009) as given in ref. 37 , that are implemented in the powerlaw package library described in references [38][39][40] . In order to decide the best fit and perform a proper statistical analysis, we explore several candidates for the distribution models: power law, power law with an exponential cutoff, exponential, stretched exponential and log-normal.
For the statistical distributions considered, we calculate the Kolmogorov-Smirnov (KS) distance between them in a pairwise fashion. This KS distance gives us a first indicator (goodness of fit) of the proximity of the data and the proposed distribution model. Then, we compare the different distributions via a likelihood ratio test by calculating the log-likelihood function of each one of the selected distributions. The sign of this ratio gives us a criterion to discriminate between distributions (see reference 37 ). After this model selection, the best two fits were the power law with an exponential cutoff (EC), with a probability density 37 : Figure 3. Global activity of taxis between zones of high demand for this service in New York City. We analyze the movement of taxis trips made in 2015 and, from the study of a grid with square zones with 100 m × 100 m, similar to the one presented in Fig. 2, we identify zones of high demand of taxis considering that at least 1 000 trips have departed or arrived from a zone. We found with this criterion N = 4 353 high demand zones. In (a) we present the origin-destination matrix for taxi trips moving between zones, the respective colorbar codifies the trip counts. In (b) we present the geographical distance between origin and destination zones; the values of the distance are represented in the colorbar.
k EC 1 min and the stretched exponential (SE) 37 : where k represents the degree, k min is the minimum value considered in the fit and Γ(x, y) denotes the incomplete gamma function. Notice that both distributions have two parameters, that we will distinguish with a superindex EC for the power law with an exponential cutoff, and with a superindex SE for the stretched exponential; however, we will not indicate these superindexes in Eqs. 4-5 to simplify the notation. We use λ EC and γ EC for the power law with an exponential cutoff and λ SE and β SE for the stretched exponential.
For the in-degrees in Fig. 4(a), the best fit is the stretched exponential with parameters β = . × − 6 730 10 in EC 6 . On the other hand, the same analysis for the out-degrees in Fig. 4(b) concludes that the best fit is the power law with an exponential cutoff with parameters 1 0000000025 out EC γ = . and λ = . × − 6 086 10 out EC 6 ; in addition, for the stretched exponential β = .
0 495 out SE and 6 834 10 out SE 5 λ = . × − . It is surprising that both exponents γ in EC and γ out EC are extremely close to the value one. Thus, both distributions are well described by the power law p(k) ∝ k −1 in some range of in and out degrees.
transition probabilities. All the information in the origin-destination matrix and in the degrees k i (in) and k i (out) allow us to analyze and understand the spatial activity of taxis as a dynamical process in a spatial directed weighted network. In this way, we can describe statistically the global dynamics of taxis in terms of transition probabilities between high demand zones of this service.
The transition probability w i j (out) → between zones i and j is defined in terms of the origin-destination matrix as: With this definition, the transition probabilities satisfy the normalization condition: With the transition probabilities w i j (out) → , we can explore the relationship between the information in the origin-destination matrix and the matrix of distances; these matrices were presented in Fig. 3. Now, to study this connection, we calculate → w i j (out) by using the definition in Eq. 6; for each value, we have the corresponding geographical distance d ij between i and j as an entry in the distance matrix D.  Table 1 for each of the years explored. The results were obtained with normalized counts using logarithmically spaced bins. In both cases, we show with different curves, the power law with an exponential cutoff p EC (k) in Eq. 4 and the stretched exponential fit p SE (k) in Eq. 5. (2020) 10:4022 | https://doi.org/10.1038/s41598-020-60875-w www.nature.com/scientificreports www.nature.com/scientificreports/ In Fig. 5, we depict the logarithm of the transition probability → w i j (out) as a function of the logarithm of the relation d ij /d 0 where d 0 = 1 Km is a reference length. In Fig. 5(a), we consider all the non-null transition probabilities w i j (out) → and distances d ij , for the annual data records of the taxi's activity in 2015; we obtain a distribution of points for distances less than a characteristic value R = 1.8 Km. In contrast, for distances greater than R, the transition probabilities are well described by a power law with an exponential cutoff relation: where continuity for d ij = R requires a = 10 c . Now, to find the best fit, we analyze the pairs  In all these cases, we analyze the non-null transition probabilities → w i j (out) and the geographical distance d ij between zones i and j. We show hexagonally binned twodimensional histograms for the logarithm of w i j (out) → and the logarithm of d ij /d 0 where d 0 = 1 Km is a reference distance. The values codified in the colorbar represent the frequencies denoted as ( ) found in each hexagonal bin. Dashed lines are used as a guide and represent the Scientific RepoRtS | (2020) 10:4022 | https://doi.org/10.1038/s41598-020-60875-w www.nature.com/scientificreports www.nature.com/scientificreports/ probabilities of transition. In a similar way, once defined c, we calculate β in Eq. 8 for d ij > R. In Fig. 6(b), we analyze the probability density of the values β and we identify a peak around β = 0.15 Km −1 .
The piecewise approximations described by the values R = 1.8 Km, c = −3.2 and β = 0.15 Km −1 are represented with dashed lines in Fig. 5. A similar behavior has been detected in the analysis of the transportation network of stations in bicycle sharing systems operating in New York City and Chicago 21 . In these cases, the value R ≈ 1 Km defines local displacements and the long-range dynamics is well described by w d In this way, in bike-sharing systems R is reduced in comparison to our findings for taxi trips; in addition, the long-range spatial activity qualitatively has similar characteristics to those observed in Fig. 5.
→ defined in Eq. 6 allow us to understand human mobility as a dynamical process in a spatial directed weighted network. Well-known results in stochastic processes apply for the transition matrix W (out) 10 . In most cases, origin-destination matrices are non-symmetric; as a consequence, it is convenient to analyze the transition matrix W (out) establishing an analogy with the Google matrix 41 , with a mathematical structure entirely general that applies to any graph or network in any domain 42 . In the following, we explore how by using this connection, the eigenvalues and eigenvectors of W (out) give valuable information to understand the movement of taxis.
The transition matrix W (out) has left and right eigenvectors. Left eigenvectors j Φ → with elements φ j (i) satisfy: are the eigenvalues of the transition matrix. Right and left eigenvectors form an orthonormal base and have the same eigenvalues. On the other hand, the stochastic matrix W (out) fulfills Eq. 7 and, by definition, the elements of T ij satisfy T ij ≥ 0; therefore, W (out) belongs to the class of Perron-Frobenius operators with a possibly degenerate unit eigenvalue λ = 1 and other eigenvalues obeying |λ| ≤ 1 (see 43 for details).
In Fig. 7(a) we plot the eigenvalues of the transition matrix W (out) for taxi trips in New York City in 2015. We use the origin-destination matrix in Fig. 3(a) and the definition in Eq. 6. The results were obtained numerically and, due to the asymmetry of the origin-destination matrix, the eigenvalues are complex numbers. In Fig. 7(a) we show the real and imaginary part of each of the eigenvalues λ i for i = 1, 2, …, N = 4 353. In this analysis, we found that only one eigenvalue satisfies λ = 1, a result that reveals that the directed network associated with the mobility between sites of high demand for taxis is connected. Therefore, the links in the network connect all the zones. This particular result can be interpreted using the terminology of random walks on networks. In this case, the movement of a random walker defined in terms of the transition matrix W (out) is capable to visit any node of the network only by moving on the links, independently of the initial configuration. As we mentioned before, the high connectivity observed in the origin-destination matrix is a consequence of considering high demand zones with a criterion that requires a high number of departures and arrivals in each zone avoiding the emergence of isolated parts. However, the approach developed is general and in other cases, similar spectral analysis of the transition matrix could be an important tool to identify disconnected parts in a transportation system.
In addition to the eigenvalues, the respective eigenvectors of the transition matrix provide valuable information about dynamical processes on networks 10,19 . In particular, the left eigenvector associated with the eigenvalue λ = 1 defines a ranking vector P → ∞ with elements ∞ P i for i = 1, 2, …, N and satisfies P W P In the study of random walks on networks, the vector P → ∞ is the stationary probability distribution. The value ∞ P i gives the probability of a random walker to reach the node i after a large number of steps 19 . In the context of the Google matrix, the vector → ∞ P determines the importance of a node in a network establishing a PageRank of the Web 43 . In the analysis of mobility with a transition matrix W (out) , the vector P → ∞ defines a ranking of the zones used in the definition of the origin-destination matrix. Due to this connection, we call this ranking "OD rank".
In Fig. 7(b), we show the results obtained numerically for the OD rank → ∞ P associated with the eigenvalue λ = 1 of the transition matrix W (out) that describes the taxi's flow in 2015. Our findings in this figure reveal a connection between the OD rank ∞ P i of a zone i and the respective in-degree k i (in) . In a similar way to the findings for the PageRank algorithm for Google, the stationary probability distribution → ∞ P is a measure of the popularity of nodes that is mostly due to the in-degree dependence; in a mean-field approximation the stationary distribution of the PageRank algorithm is given by 44 : where 0 ≤ q ≤ 1. Searching the optimal value q * that minimizes the quadratic error S q P P q ( ) ( ( ) ) , we get for the best fit: Fig. 7(b) we illustrate the approximation given by Eq. 11, for q = 0 and q * = 0.062, obtained for the best fit. However, Eq. 11 is a mean field result and important deviations may appear 10,[44][45][46] . The result given by Eq. 11 makes sense in the description of taxis since the importance of a high demand zone can be defined in terms of the number of taxi trips k (in) that arrive at this specific location. For example, in our schematic illustration presented in Fig. 1(a), now we understand that the bars with the value k (in) determine the importance of the zones.
The transition probability matrix W (out) defined in Eq. 6 captures all the information about the system's global activity. We think that an OD rank of the zones defined as → ∞ P can be a valuable measure in the analysis of different transportation systems and a complement to other types of ranking algorithms introduced to determine location attractiveness incorporating geographic considerations into the PageRank algorithm 47-49 . Random walk strategy. The results obtained before for the relationship of the transition probabilities describing the flow of taxis between zones and the geographical distances separating these locations, suggest that the spatial dynamics can be approximately described by a model with constant transitions to zones in a local Figure 7. Numerical analysis of the eigenvalues and OD rank of the transition matrix W (out) . We analyze the transition probability matrix for the taxi's flow in 2015 with origin-destination matrix presented in Fig. 3 with N = 4 353 high demand zones. In (a), we show the eigenvalues λ of W (out) satisfying Eq. 9. In this way, we have 4 353 values represented in the complex plane with dots; in the inset, we depict the results for the eigenvalues in a region close to the origin, where we observe more eigenvalues with a non-null complex part. In (b) we plot the components ∞ P i of the eigenvector → ∞ P with eigenvalue λ = 1; we represent the numerical values of ∞ P i in terms of the respective degree k i (in) for i = 1, 2, …, N. We also show the values ∞ P q ( ) i obtained with Eq. 11 for q = 0 and the best fit q * = 0.062.

Scientific RepoRtS |
(2020) 10:4022 | https://doi.org/10.1038/s41598-020-60875-w www.nature.com/scientificreports www.nature.com/scientificreports/ neighborhood within a distance R, and a long-range dynamics defined by probabilities of transition proportional The analysis of more than a billion trips reveals a particular emergent pattern in the spatial activity. The movement of taxis between high demand zones can be classified into two types of trips with particular characteristics illustrated in Fig. 8. We have local displacements for which a taxi departs from a high demand site and the probability of moving to another site of high activity is independent of the distance that separates them if they are located at a distance less than a value R. On the other hand, there may also be long-range displacements for which the separation between origin and destinations require distances greater than R. For this type of movements, we find that the probability of having a long-range trip depends on the distance and these particular transitions have characteristics observed in truncated Lévy flights.
In this way, to describe the global activity of the taxi's mobility we use the model: In this model, β and R are positive real parameters. The transition probabilities defined in Eqs. 13 and 14 are illustrated in Fig. 8. The radius R determines a neighborhood around each zone where the trips occur with equal probability to move from the initial site to any of the high demand zones in this region. Therefore, the displacements are independent of the geographical distance between origin and destination. That is, if there are S sites inside a circle of radius R, the probability of going to any of these sites is uniform. Additionally, for places beyond the local neighborhood, for distances greater than R, the transition probability decays as a power law with an exponential cutoff of the distance and is proportional to β In this way, the parameter R defines a characteristic length of the local neighborhood and β controls the probability to have long-range displacements. In particular, in the limit β → ∞ the dynamics becomes local. We introduced a similar model with long-range transitions proportional to α − d ij (α > 0) in reference 8 in the context of human mobility and encounter networks. In this case, the resulting dynamics can be similar to a rank model [50][51][52] and a gravity model 3,[53][54][55] . It is worth mentioning that the inverse of the parameter β in Eq. 14 gives us a characteristic distance; this exponential cutoff takes into account the finite size effect associated with a finite system like New York City.
In our previous analysis in Fig. 5, we found that R ≈ 1.8 Km. This value defines what we understand as a local neighborhood for this transport system. On the other hand, for distances d ij > R, the probability to have a trip to a zone is highly influenced by the geographical distance and this long-range dynamics is determined by the values In the following part, we explore the predictions of this model for the annual global activity of taxi displacements in New York City by using the parameters R = 1.8 Km and β = 0.15 Km −1 found in the analysis of the taxi's flow between high demand zones. In addition to Eqs. 13 and 14, that model the displacement between high demand zones, it is important to consider that these zones have different relevance in the whole dynamics, i.e., a trip can start from different zones with non-uniform probabilities. This fact is well described by the values of the out-degree k i (out) defined in Eq. 2 that gives the number of trips with origin in the zone i. In addition, from the results in Fig. 4, we know that the values k i (out) follow a hierarchical distribution with probabilities that decay as p(k) ∝ k −γ e −λk where k represent the values of the out-degree. This result is observed in the annual datasets from 2009-2015. In this way, we simulate the dynamics of multiple taxis that start from an initial zone chosen randomly with a probability proportional to the values k = that quantify the importance of each zone in the city. Then, Figure 8. A schematic illustration of the mobility of taxis between high demand zones. There are two types of trips from a particular location i: First, to a site j inside a circular region of radius R centered in the location i, the probabilities to have a trip to these zones are constant; and, second, a trip to a zone k outside the circle of radius R. In this case, the probability to have this long-range movement decays as a power law with an exponential cutoff proportional to e d where d ik is the geographical distance between i and k.
a displacement is generated randomly from the origin site to a final zone by using the transition probabilities in Eq. 13; this algorithm is repeated to generate, through Monte Carlo simulations, the same number of displacements between high demand zones, as reported in Table 1.
In Fig. 9 we present the statistical analysis of the taxi's displacements d generated randomly and the real values considering the activity in New York City in 2015. In our simulation, we generate 128 984 657 random displacements following the model in Eq. 13, with β = 0.15 Km −1 and R = 1.8 Km. Our findings show an agreement between the predictions of the model and the real dynamics. However, we observe that the predictions do not agree with the real data around d = 10 Km and d = 20 Km; this is a consequence of the singular dynamics induced by the two airports in New York City. An accurate modeling capturing the effects of these very attractive sites in a city requires modifications to the model explored. This fact is also visible in Fig. 5 where, for distances around 20 Km, we see different values of the transition probability that are not described by a model with long-range trips following − . Finally, we repeat the same procedure to compare the predictions of the model with respect to the real data for taxi's activity from 2009 to 2014. Our results for Monte Carlo simulations are presented in Fig. 10. We observe the same characteristics found in Fig. 9, with a good agreement between model and the data. The number of locations of high demand N and the number of displacements analyzed for each year are reported in detail in Table 1. The results in Table 1 also reveal that in average, in a year, approximately 43% of the trips are local movements for which the geographical distance d ≤ 1.8 Km, the rest of the trips are non-local with d > 1.8 Km.

Discussion
In this research, we explore the massive records of more than one billion taxi-trips in New York City from January 2009 to December 2015. With this dataset of seven years, we generate an origin-destination matrix that has detail information of a vast number of trips. The mobility in New York City can be described as a directed weighted network that connects different zones of high demand for taxis. Each zone is characterized by the number of trips that arrive or depart from it and corresponds to nodes in the network. The arrivals and departures are the in-degrees and out-degrees of the directed network, and the flow gives different weights to the links of this spatial network.
We present a statistical analysis of the travel distance of each trip and found a long-range distribution that is almost the same for each of the seven years studied. On the other hand, the degree distributions, for the in and out degrees are, respectively, well modeled by a stretched exponential and a power law with an exponential cutoff. By defining the transition probabilities between zones, given by the origin-destination matrix and the out-degree, we are able to obtain a rank, called "OD rank", analogous to the page rank of Google. We calculate the spectrum of eigenvalues and the main eigenvector, which is related to the in degree. The components of this eigenvector give the more relevant and attractive places in New York City, in terms of taxi trips.
The dependence of the transition probabilities with the distance between zones is obtained from the dataset, and based on that, we introduce a model that captures the global dynamics of trips. The data and the model describe, for short distances, a local dynamics independent of the spatial distance, and, for large distances, a dynamics that decays with distance as a power law with an exponential cutoff. The data agrees quantitatively with Monte Carlo simulations based on our model. Finally, considering the taxi trips as a proxy of human mobility in cities, it might be possible that the long-range mobility and other features found for New York City would be rather general, and thus we expect a similar behavior in other large cities around the world for which these ideas can be applied as well.

Methods
Dataset description. In this section, we present a global description of the records explored to study the spatial dynamics in New York City. We use data for the activity of taxi trips from January 2009 to December 2015; these datasets are available to the public by the Taxi and Limousine Commission in the New York City open data website 56 . The data available include information for all taxi trips in New York City when the taxis are in service.    Table 2. The records comprise several fields capturing pick-up and drop-off dates and times, pick-up (origin) and drop-off (destination) locations, itemized fares, rate types, payment types and driver-reported passenger counts 56 . Now, to complement the information in Fig. 2, and identify other global characteristics in the taxi's spatial activity, we analyze the geographical distance d between the origin and destinations in each trip calculated from the longitude and latitude coordinates of these locations reported in the database. Here, it is worth mentioning that other types of distances can be implemented; in particular, the distance of the path in the street network connecting origin and destinations. In fact, powerful techniques have been introduced exploring taxi trips in New York City to estimate the driving distance based on the origin and destination coordinates 59 . However, due to the different paths that a taxi can follow to carry out each trip, in the following we use the geographical distance d. In Fig. 11, we present the statistical analysis of the geographical distances d. We depict the frequencies f(d) of the displacements obtained from uniform bin counts with Δd = 500 m for taxi trips. Different markers show the results for the analysis in a year. We can see that the frequencies f(d) maintain the same characteristics from 2009 to 2015, and the statistics reveal three important intervals: the first for d < 1.8 Km with higher values of the frequencies, a second interval for 1.8 Km ≤ d < 20 Km where f(d) gradually decays and finally, for distances around 20 Km, we identify a peak that decays rapidly with the distance; this peak is associated with large displacements from Manhattan to the JFK airport (as a reference, the geographical distance between Times Square and the JFK airport is 20.6 Km). In a similar way, we identify another relative maximum at d = 10 Km: this increase in the frequencies is associated with trips between Manhattan and La Guardia airport (with d = 9.8 Km between Times Square and this airport). These are examples of how important locations can induce long-range dynamics in the taxi's mobility. In this case, the two airports in New York City influence the taxi transportation mode in the whole city. This important feature has been observed in other cities with airports located at the city's periphery (a particular case is reported in 31 ).
In Table 2, we summarize the global information found for the spatial dynamics per year. We present the number of taxi trips analyzed, the average distance, the largest distance traveled as well as the fraction of trips with distances at different intervals. From the information in this table, when we examine the complete records from 2009 to 2015, we observe that 41.33% of the taxi trips have displacements with d less than 1.8 Km, whereas a 57.49% of the trips involve long-range displacements in the interval 1.8 Km ≤ d < 20 Km, and only a 1.18% of the trips have d greater than 20 Km. The average displacement of trips is = . d 3 3 Km and the maximum value We depict the results for the length of the intermediary path and the geographical distance between these intersections; the distances are indicated by the colorbar. In (c) we present the hexagonal bin counts for the geographical distances and the respective length of the intermediary path. We depict a dashed line, with unit slope, that represents the case when the two distances are the same. Clearly, since the intermediary path is always greater or equal than the geographic distance, we only have data in the lower triangle of the figure. We show with a colorbar the frequencies for the values found in each bin. The street map, intersections, and intermediary paths were obtained and analyzed with the OSMnx package 57,58 . (2020) 10:4022 | https://doi.org/10.1038/s41598-020-60875-w www.nature.com/scientificreports www.nature.com/scientificreports/ observed in the records is 51.87 Km. All these quantities give us a first characterization of the spatial activity of the taxi transportation mode.

Geographical and shortest path distances.
In the analysis of the information described before, we use the geographical distance d between origin and destination. This election is based on the fact that we only know the geographical coordinates of origins and destinations for each trip. However, another important quantity to consider is the length of the intermediary path that the vehicle follows on the street network. The information of the street network, the length, and direction of each street and all the intersections, can be obtained from different sources like OpenStreetMap 60 or generated by using specialized algorithms (see for example 59 ). In general, the length of the geographical distance is less or equal than the length of the shortest path between two points in a city. In Fig. 12 we explore this relation for all the intersections in Manhattan's street network. We analyze the information available in OpenStreetMap 60 and the OSMnx Python package 57,58 to generate the street network depicted in Fig. 12(a). From this structure, we obtain the geographical coordinates of 4 409 intersections. In Fig. 12(b) we calculate all the distances between these intersections, taking into account the length of the intermediary path, and the respective geographical distance. The results are presented as matrices for which the entry l, m represents the respective distance between intersections l and m. The two matrices are similar; however, the matrix with intermediary paths is asymmetric since it includes the directions of the streets. In Fig. 12(c) we explore the relation between the two distances by plotting all the values presented in Fig. 12(b). The results reveal that a high fraction of the values is close to a linear relation. Similar results apply for the whole city and, in this way, the main features of the global activity of taxis can be analyzed by using only geographical distances between origin and destination. However, in other contexts, a description of the complete path followed by the vehicle is necessary. See refs. 35,59 for a detailed discussion and models for taxi's mobility at the level of intermediary paths.