Preserved layout features embedded in road network development

Historical maps predate the digital maps obtained from satellite images, and are thus useful for extending the temporal scope of the study. The quality of the maps gives a glimpse as to the degree of scientific and artistic competency of the population during the period they are created [40]. Needless to say, hand-drawn historical maps, while correctly depicting the structure of the road network, lack the geographical accuracy required for quantitative analyses. As such, a historical map should first be subjected to proper georeferencing before digitization and data extraction.

In this work, we use historical maps of the city of Manila at various points in its historical development. The years and the data sources are presented in table 1. To provide a good representation of the different periods of development, we selected maps from the Spanish colonial period (1875 and 1894), the American colonial period (1908 and 1945), and the postwar republic (1991, and the contemporary 2019 map).

The raw historical map data are obtained by digitally scanning an actual physical copy of the map using the Niji scanner with a resolution of 600 dpi. Meanwhile, the contemporary map is extracted from OpenStreetMap (OSM) [41] imported through QGIS. The OSM data is accurate as of January 2019. For our purpose, this contemporary road network will serve as the baseline reference to evaluate the networks extracted from the maps.

The historical maps, now being handled as digital images, are georeferenced with respect to the contemporary GIS map from OSM using the OpenLayer plug-in in QGIS. First, we create a convex hull outlining the borders of the city to limit the extent of the geographical matching. Afterwards, we establish ground control points (GCP) that are deemed to be corresponding locations in the historical and the OSM map. Whenever possible, road network features that are present in all historical maps are included in the GCP to ensure the consistency and accuracy. Using the WGS 84 Pseudo Mercator as the reference coordinate system, we select either a thin plate spline or polynomial 3 transformation and nearest-neighbour resampling method for each historical map to warp the planar image onto real space. The resulting digital image in geotiff format has every point corresponding to a real geographical location denoted by the longitude ϕ and latitude θ coordinate. An illustration of the process is presented in figure 1.

Zoom In Zoom Out Reset image size

Using the georeferenced historical maps, shapefiles of the Manila road network are drawn in QGIS. The roads are manually traced by following the centerlines of the hand-drawn roads from the map. While a combination of several techniques from image processing and machine learning can be used to make this process automated, we opted to manually trace the roads in QGIS for two reasons: (1) to be able to follow the tracks of roads by their names; and (2) to ensure substantial accuracy, especially in regions where the roads are no longer extant. Needless to say, the process also has drawbacks, chief of which is the fact that the egregious systematic errors of the cartographer can never be corrected.

The road network shapefiles are then extracted for the nodes in QGIS. Majority of the nodes mark intersections of two or more distinct roads, or bends for winding roads. Additionally, some nodes can be found along a stroke of road. This happens when long isolated portions of the road have been segmented to denote important points, such as locations of national monuments or local markers, or in places where the road widens or constricts. There are also cases when a stroke of road consists of two distinctly named streets; in such a case, the node marks the boundary between the two.

The road network formation is a complex self-organizing process [27, 30], one that forms fractals [10, 29, 42] and other morpohological patterns [20–22] in space. In quantifying the evolution of the road network, one can view it as a collection of individual roads r, as named by the city. In this case, each road can be measured based on its geographical length L_r. The entire network for each map representing the different periods can be characterized by: (i) the total road length L_sum = ∑_rL_r; (ii) the average road length ⟨L⟩ = L/N, where N is the number of roads obtained in the map over that period; (iii) and the road length area density, ρ_L = L_sum/A, where A is the area of the convex hull bounded by the existing roads over the given map period.

Because of its structure, the road networks can also be understood using the techniques developed in complex network theory [43]. Mathematically, graphs are constructs made up of nodes whose interactions are modeled by a connection or edge drawn between them. For this work, we create the primal representation of the road network wherein an individual road is a path in the network, the geographical intersections or turning points are the nodes i, and the physical portions of the road between intersections i and j are the edges e_ij. With this, we can think of the network at any given map period as a collection of nodes and edges, and we can track: (i) the number of nodes n_i and edges n_e; and (ii) the node and edge densities over A, p_i = n_i/A and P_e = n_e/A, respectively.

Apart from viewing it as a simple collection of elements, we represent the road network as a weighted graph wherein each edge has a weight value w_ij corresponding to the geographical length of e_ij. This introduces real-world traffic implications and node-to-node relationships. In traversing the network, the path of minimum distance required to reach j from i is called the shortest path length l_ij, which is reported in terms of the number of node-to-node jumps.

Additionally, representing the road network using graphs allows us to utilize the tools of complex network theory for characterizing the state of the network over time. Traditionally, the importance of each node in the network [39] are quantified via different measures of centrality, with each having different interpretations based on the specific aspects of the network connectivity being measured. For our purpose, we track the closeness and betweenness centralities, two measures that are related to the physical layout of the nodes in space. The closeness centrality ψ is given by

$$\begin{equation*}{\psi }_{i}={\left[\sum {l}_{ij}\right]}^{-1}\end{equation*}$$

for all j ≠ i; this metric gives the structural center of the network, as the node with highest closeness centrality has the minimum shortest paths, ${l}_{ij}$, to all the other nodes. On the other hand, the betweenness centrality is given by

$$\begin{equation*}{\delta }_{i}=\sum \left[{\sigma }_{hj}\left(i\right)/{\sigma }_{hj}\right]\end{equation*}$$

for all h ≠ j ≠ i. It measures how many times the node i has appeared, σ_hj(i), in the collection of all shortest paths in the network, σ_hj. The high-δ nodes, therefore, may be thought of as bridges where efficient connections between peripheral nodes are passing through. In the subsequent analyses, the structural differences between the historical and contemporary maps can therefore be quantified using these metrics.

Figure 2 shows the digitized maps from various years (colored roads) laid over each other and on top of the contemporary road network structure in the region from OSM (gray roads). A closer inspection of the maps reveal road patterns that are not directly coinciding. In very limited cases, we find instances of road translations by 75 m from contemporary locations. As mentioned previously, the georeferencing procedure is subject to possible errors due to manual cartography, resulting in minor deviations in road network patterns. Additionally, the postwar reconstruction resulted in actual geographical changes in some roads, which explains some differences between pre-war and post-war layouts. We find that these changes have minimal effect on the characterization, especially in computing for the topological properties.

Zoom In Zoom Out Reset image size

Due to the geographical constraint (the Manila Bay to the left), the spatial growth had generally extended eastward across the observation period. Interestingly, however, we find the development of new features that coincide with key events in the history of the city. At the bottom panels of figure 2, we present the successive maps, with the older map in black and the later map retaining its color, to be able to observe the changes that transpired. The 1875 map clearly shows the road network plan for Manila developed by the Spanish since the late 16th century: the core city, near the opening of the river to the bay, is grid-like and walled. This is called the Intramuros ('within the walls') district, and is still extant and preserved.

The period between 1875 and 1894 correspond to the waning periods of Spanish colonization of the Philippines, before the Philippine revolution of 1896. As such, we find very minimal growth between this period, with just a few pockets of grid-like developments sprouting in 1894 from an already stable road network layout. The period between 1894 and 1908, in contrast, can be characterized to be very dynamic, as reflected in the changes in the physical network layout. Interestingly, this is associated with the corresponding socio-political conditions at the time: Spain ceded control of the Philippines to the United States in 1898, thereby beginning the American colonization and, with it, massive construction in the capital. In 1904, in an attempt to resolve the status of the friars, the United States bought the major portions of the friar's holdings of 1,660 km², half of which was in the vicinity of Manila [44]. This largely increased the territorial coverage of the city available for infrastructure growth. These major efforts in building new and longer roads expanded the road network span but lessened the density of roads built. Merely ten years after the Americans took over, the road network structure grew significantly. In the figure, we find the development of what would now be the South Port of Manila to the west of Intramuros, and other grid-like developments southward in the 1908 map. Incidentally, this map also showed exploratory mechanisms towards the east, as part of the efforts to connect the capital to the nearby provinces.

Before the end of World War II in 1945, we observe a heavily homogenized space in the Manila map. The grid-like portions in the northeast coincided with the construction of the planned Quezon City that started in 1939, which was designed to be the new capital. The boundary between these two cities, containing residential communities of 'New Manila', are heavily densified to recover uniformly sized blocks, inspired by the utilitarian network layouts typical of other American urban areas such as Manhattan in New York and Chicago [45]. Between the 1945 and 1991 maps, not much changed in the network layout, save for the road network development at the North Port of Manila to the northwest of Intramuros. However, between the 1991 and 2019 maps, we find major road network growth in the northwestern part of the city due to reclamation efforts [46]. Incidentally, the development of these regions correspond to the period of rapid economic development of the Philippines, which resulted in the rise of infrastructure spending by the government and private entities. These areas now host the bigger port operations and large business and mixed-used establishments.

Apart from the qualitative description provided by the spatial mapping, the metrical and topological properties for all the years are plotted in figure 3. The top panels denote the raw values of the metrical properties (a) N and (b) L_sum, and the topological properties (c) n_i and (d) n_e. As expected, the road network grows in spatial extent and topological (network) complexity in time. Comparing the spatial and network growth rates of Manila with other cities in literature such as Groane, Italy [27] and Paris, France [30] reveals a faster growth rate for Manila over similar periods, most likely due to the constrained geographical limits. In particular, the heavy spikes observed between 1991 and the contemporary data for the number of roads (figure 3(a)), number of nodes (figure 3(c)), and number of edges (figure 3(d)) signify the extent and rates of the reclamation and the road-building during this relatively short period. However, the total length of roads (figure 3(b)) is expected to show a tapering off, as the space cannot be entirely converted into roads. These results suggest the type of development over the last few decades: new city areas that are added by reclamation are developed to have multiple short roads that divide the space into regular blocks. These punctuations of highly-planned development over the predominantly self-organized growth are also easily captured by the visual map comparisons from figure 2.

Zoom In Zoom Out Reset image size

To further investigate the growth, we obtain the effective spatial area A spanned by the roads within the convex hull of the Manila map for the particular period in consideration. The bottom panels of figure 3 are measures of spatial density (i.e. raw values divided by the effective area being covered) of the metrical and topological properties over the years being considered. If the expansion is uniform, we expect to see flat regions in the spatial density plots. In figures 3(e)–(h), however, we find generally increasing trends, especially at the later years, signifying a densified development.

Despite the generally increasing trends in the density values, interestingly, both the road count density N/A (figure 3(e)) and road length density ρ_L (figure 3(f)) show dips in 1908. These signify exploration; the road network increased in longer and curved roads following the terrain instead of short, densified roads forming blocks. The decrease in N/A in figure 3(e) is associated with the fact that the Americans extended some pre-existing eastward roads, increasing the extent of the convex hull without significantly increasing the number of roads. On the other hand, the decrease in ρ_L for 1908 suggests that the exploratory mechanisms dominate over the densification, resulting in large swaths of areas without roads in the east, as can be observed in figure 2 between the 1894 and 1908 maps. Incidentally, a similar decrease in ρ_L is observed for the contemporary map (figure 3(f)). This is attributed to the fact that the city has grown in area due to reclamation, and that the newly formed roads are well-planned, resulting in an increase in the convex hull area that is not commensurate with the increase in total road length.

Finally, the ρ_i (figure 3(g)) and ρ_e (figure 3(h)) count densities are found to follow mostly the trends of the actual counts in (c) and (d) respectively, signifying that the road network complexity is independent of the actual area spanned by the roads in space. This, again, is illustrative of the densification process; in many instances, the increase in the intersections (nodes) and the corresponding street segments (edges) happens in the interstitial spaces between the larger blocks that are already present in the network and less likely to happen for peripheral areas not yet reached by previous roads.

Despite the various events that have happened throughout its historical past, and even in the absence of a centralized planning body in modern periods, we observed statistical regularities in the Manila road network over time. Figure 4(a) shows that the road lengths follow unimodal fat-tailed distributions P(L) throughout the years. The distributions are highly skewed and non-normal, implying that the creation of roads is not simply a random process. Similar trends have been observed for other data sets such as London, United Kingdom [28], Amiens, France [21], and even for the global road network (GRN) [31], which the researchers opted to model with different fitting functions. For London [28], the trend is fitted by a log-normal behavior. On the other hand, the length distribution for Amiens [21] was fitted with a combination of an exponential and a powerlaw. The choice of fitting function indicates that there is no pre-existing typical road length scale; road creation in cities is a multiplicative process of cutting through bigger blocks until a considerable homogenized state [27]. Interestingly, the GRN [31] is found to collapse into a universal form upon rescaling by the characteristic quantities of the road network, regardless of the type and function of road. These studies strengthen the case for identifying the road network as a complex system that produce emergent statistical features.

Zoom In Zoom Out Reset image size

In the case of Manila, we observe generally unimodal distributions, as shown by figure 4(a). The unimodality of the distributions suggests a similar mechanism for the road creation even in the absence of a strong, centralized planning and strict zoning rules. In the case of Amiens, a bimodality is observed due to local modifications [21]. However, one feature of the distributions shown in figure 4(a) that is deemed to be more important are the fat tails that appear to follow similar trends. The persistence of the fat-tailed behavior for all the years considered implies that very long roads were deliberately created during the exploration stages that served as main trunks from which shorter local roads emerge during the later densification. And, throughout the evolution of the network, these longest roads are rarely, if ever, being cut into shorter roads, thus preserving the statistics at the tails.

Because of their particular significance, we first focused on an approximate description of the tails of P(L). In figure 4(b), we rescaled the distributions by their correponding mean to remove the variability for all the years. Interestingly, the tails of the distributions collapsed under a single curve, while the short road regimes did not. The tails of the distribution all scale as P(L) ≈ L^−ζ, where ζ = 2.90. This is consistent with previous findings for the 2016 network of Manila [22]. This observed invariance is a result of the consistency of the hierarchy of roads in the network [20]. Note that the Manila road network grew organically over the years; it is impossible to have had a continuous, coherent plan that spanned the entire period being observed here, due to the various historical conditions of the city. The collapse into a universal scaling with exponent ζ is an emergent property that is manifested in the temporal domain of the road network evolution. Incidentally, previous works have shown that ζ is also a universal scaling exponent of the P(L) tails for all the cities in the Manila conurbation [22].

The shorter-than-mean regimes in P(L) give us an insight into the degree of the densification process. As shown in figures 4(a) and (b), the fraction of the very short roads in the 2019 P(L) significantly increased compared to the other years, even upon rescaling with the mean. This suggests that the Manila city government has fairly recently embarked on a major road-building process with very densified (i.e. short roads) geometries, consistent with the sudden spikes in figures 3(a) and (b) and observed in figure 2.

To further highlight the extent of the creation of very short roads, we use the rescaling procedure utilized by Strano et al [31] for the GRN. The universal scaling function F(⋅) may be obtained by plotting L^γ ⋅ P(L) with L/⟨L⟩^α, where γ = 1/(2 − α) ensures normalization. In figure 4(c), the case of α = 1 (and, thus, γ = 1) is presented, i.e. rescaling by the simple mean road length ⟨L⟩ of the data. As expected, all the distributions collapse into the universal scaling function upon the proper rescaling. Noticeably, the contemporary road network shows higher regimes at the very short road length scales, that deviates significantly from the collapsed plots for the other maps. This appears to indicate the nature (i.e. very short roads) and extent (i.e. significantly higher preponderance) of the recent densification of Manila, which, by comparing with figure 2, is highly planned.

The densification is not confined to the reclaimed areas, however. To probe the road growth within the preexisting road layout, we plot in figure 4(d) the contemporary road network, and mark the convex hull of the 1875 network, the oldest historical map that we have. We establish an arbitrary boundary between the regions within and outside the 1875 convex hull and present the conditional road length distributions in figure 4(d). A quick comparison between the original road length distribution of 1875 and that of the modern network within the 1875 bounds reveals an order-of-magnitude shift in the mode towards the shorter L values, and an increase in the actual probability density P(L) at this new modal value. This shift quantifies the within-city densification, which, along with the planned development at the reclaimed areas, contributed to the short-L regimes in the P(L) distribution of the 2019 network.

We also looked at the distributions of the closeness ψ and betweenness centralities δ of the networks from the different periods and present them in figure 5. As stated previously, the ψ and δ values are the most associated with physical space among all centrality measures. Previous works have indicated that their distributions (particularly, that of the δ) show statistical consistency across various city data sets [33].

Zoom In Zoom Out Reset image size

The P(ψ) distributions (figure 5(a)) are expected to shift to smaller values simply because of its definition as an inverse of the spatial extent; the average shortest distance increases with continued spatial growth. Upon a simple rescaling with the mean (figure 5(b)) or by the scaling function (figure 5(c)), the characteristic modes appear to align near the value of the mean, or at ψ/⟨ψ⟩ ≈ 1. One noteworthy feature of the contemporary map that is observed in figure 5(b) but is not readily apparent from figure 5(c) is the large proportion of ψ → 0 relative to the other periods. This adds a detail into the picture of the densification process described earlier. The sudden surge in the number of nodes with very small ψ suggests that the densification phase—already visualized in figure 2, and found to be due to very short roads in figure 4(c)—happened at the periphery, significantly far away from the other pre-existing nodes. This, of course, is expected of the reclamation areas.

Figure 5(d), on the other hand, shows the P(δ) across the different periods, with a general shape that resembles the one observed previously for other road networks and for random planar graphs [33]. A simple rescaling by the mean in figure 5(e) significantly aligned the tails of the distributions, while it can be shown that a rescaling by the number of nodes n_i results in the matching of the secondary peaks [33]. When the rescaling is done to reveal the relative frequency distributions (figure 5(f)), the contemporary network reveals a significant proportion of very small δ nodes, which may be attributed to the fact that the planned development at the reclamation areas suddenly increased the number of nodes with low importance in terms of bridging the others efficiently.

These results paint a picture of continuing development of the road network despite the geophysical constraints faced by Manila. As observed previously, the densification of Manila, as measured by the significant increase in the P(L) of very short roads (figures 4(a)–(c)) happened from (i) within the existing road boundaries (figure 4(e)); and (ii) at the periphery (figure 5(b)). The fact that the P(δ) also shows an excess of very small values (figure 5(f)) suggests that the peripheral densification dominated the statistics: if the development were from within, the P(δ) would only rearrange the most important bridging nodes and not produce regimes having very small importance.

After discussing the evolution of the road network through the historical map, we focus our attention to the contemporary map. Since road networks are the oldest extant elements in the urban landscape, we expect that the historical development of the city must be embedded in their architecture, and can be extracted automatically from the current network.

In figure 6, we present the locations of the nodes with centralities above the percentile Q (black points, middle and bottom columns) in the 2019 map. To provide context, the historical road networks (black lines, top column) are also presented, along with their convex hull areas (colored overlays). We observed stark differences in the spatial arrangements of the highest ψ and highest δ nodes.

Zoom In Zoom Out Reset image size

The spatial configuration of the highest ψ nodes (middle panels) is found to mimic the layout of the historical road networks. This indicates that the development of the road network preserves the historical core; new pathways in the network aim to connect new locations to these central locations. For the case of Manila, this historical core is found to be around the Intramuros area, the oldest district in the city. Additionally, the spatial arrangement of the highest ψ nodes can also indicate sudden, well-planned developments within the city. When the cluster of nodes with the highest closeness centralities deviates from the area defined by the core historical portions of the network, the directionality of the shift should point to the regions of sudden, densified construction that did not undergo the self-organized process of expansion and gradual block homogenization. For example, in the Manila road network, the high-ψ clusters show a slight bias to the northeast, indicated by the grid-like regions that are outside of the convex hull of the historical maps. This, of course, is introduced by the planned construction of Quezon city, as discussed previously.

The automatic identification of the historical core regions from just the high-closeness-centrality node clusters in the contemporary map has practical significance, as it gives an easy window for highlighting the risk-prone areas in terms of urban development. The oldest quarters in the city contain heritage structures that are oftentimes preserved, making them inflexible to modern urban development. More importantly, these districts are oftentimes created without regard for the eventual rapid pace of urbanization, making them prone to traffic congestion or urban decay. For example, the ψ > Q90 clusters in Manila are found to correspond to the Binondo district, one of the oldest Chinatowns in the world and one that organically grew to accommodate foot and horse traffic.

In contrast, the spatial profile of the highest δ nodes of the contemporary map (bottom panels) reveal no such comparable patterns as the historical structure of the road network. The spatial arrangement of connectivity is expected to vary radically as new developments are being added into the network. However, it has been suggested that these high-δ nodes indicate possible locations of traffic congestions [33]. In figure 7, we obtained traffic information in the form of annual average daily traffic (AADT) from the data of actual car counts by the Metri Manila Development Authority (MMDA), Philippines from 2012–2018 (figure 7(a)) [47] and a model projection of volume vs capacity for 2030 by the Japan International Cooperation Agency (JICA) [48]. Superposing the traffic data with the δ > Q90 nodes from figure 6 reveals that the traffic chokepoints indeed lie along the paths of the highest-δ nodes.

Zoom In Zoom Out Reset image size