Allometric scaling, size distribution and pattern formation of natural cities

There has been a strong interest in more cogent definitions on economies of scale to reveal general urban growth laws and to develop urban performance metrics. Unstructured data, including satellite images, will provide us with new sources to do so by defining cities as aggregates of human activities. Such a uniform definition on the basis of nighttime light clusters is more consistent and economically meaningful than administrative or official boundaries. In comparison with patterns of city systems found in traditional census data, we examined the allometric scaling, size distribution and fractal geometry of natural cities. It can be concluded from our empirical analysis on regional, country and continental scales that a super-linear scaling between lightness and area with a stable exponent across different low light threshold levels generally holds for natural cities. But Zipf’s Law does not always apply over the whole range of lightness thresholds. Furthermore, we build a model based on the simple geometric matching mechanism to reproduce the self-organized formation process of nighttime light patterns. The statistical properties including allometries, size distributions and fractal geometries generated by our model are in good agreement with empirical evidence. These findings have profound implications for understanding the effects of simple aggregation behaviour in primitive stages of city formation and the urbanization process.


Introduction
F or the first time ever, more than half of the world's population lives in cities, and this proportion continues to grow. By 2030, 6 out of every 10 people will live in cities, and by 2050, this proportion will increase to 7 out of 10 as some 2 billion people move to cities (United Nations Human Settlements Programme, 2010;World Health Organization, 2010). Rapid urbanization and accelerating socioeconomic development have generated global problems from climate change to incipient crises in food, energy and water availability, public health, as well as the global economy. With more people leading urban lives and the number and size of cities expanding everywhere, it is increasingly important to know more quantitatively the way a city organizes itself. Identifying universal patterns in cities is extremely crucial for furthering our understanding of urban dynamics, and may help to manage many contemporary global challenges such as the use of natural resources and the growth of urban poverty.
Cities as typical self-organized systems, exhibit universal macroscopic patterns including allometric scaling, scale-free distribution and fractal geometry (Portugali et al., 2000). In the first place, it has been statistically proved that if all cities follow some proportional growth process, their distribution (at least in the upper tail) will automatically converge to Zipf 's Law (Zipf, 1949;Gabaix, 1999;Córdoba, 2008). However, since the key empirical article in 1980s (Rosen and Resnick, 1980), there have always been disputes whether Zipf 's Law can be observed in population of most cities (or urban agglomerations) on regional, country, continental and even global scales. Contrary to the common belief, from an underlying power law distribution of city sizes and a simple random sampling, one cannot obtain a ranksize rule described in Zipf 's Law for the largest cities (Cristelli et al., 2012). Thus any examination concerning the distribution of city population or areas is required to be consistent and complete in sampling (called "coherence" by Cristelli et al.) before any valid conclusion can be arrived at.
Self-organized cities as consumers of energy and resources, and as producers of artefacts, information and waste, have often been compared with biological entities (Macionis and Parrillo, 2004). In living organisms, various macro variables can be well predicted by the body size since mechanical similarity requires allometric growth (Gould, 1966;West and Brown, 2005), and such scaling laws have also been found for city systems. The exponents concerned fall into three categories: linear associated with individual human needs, sub-linear with infrastructure and super-linear with social interactions (Bettencourt et al., 2007). In contrast, the allometric scaling investigated on different scales, such as the ones among countries or provinces, yield divergent conclusions. For instance, it has been found that population increases sub-linearly with area while GDP increases linearly with population for countries (Zhang and Yu, 2010).
From a morphological perspective, cities can also be viewed as typical fractals in that their structure features statistical selfsimilarity Longley, 1986, 1994). Therefore, their complexity degree can be quantitatively evaluated, in general, using fractal dimension. For city systems, fractal dimension also characterizes the hierarchy of different subcentres or clusters across many scales (Batty, 2008).
Although the aforementioned universal patterns have been extensively documented for city systems, all these scaling exponents (characterizing Zipf 's Law, allometry and fractal properties) are notably sensitive to system description specified by consistently defined boundaries. This is particularly the case for indicators that have a non-linear relationship with population size (Arbesman et al., 2009). However, almost all the primary evidence was gathered with cities delineated, somewhat arbitrarily, by political or geographic boundaries on the basis of administrative, legal or historical status of a particular city. Such a problem of ambiguity in the definition of cities and difference in standards adopted across nations gives rise to incongruous comparisons on large scale. In fact, limited by the boundaries in the sense that city governments have jurisdiction over, traditional census data and statistics are mostly collected in an area which is not economically meaningful (Angel et al., 2005). Many challenges need to be overcome before identifying patterns of behaviours or devising a theory across cities, among which the most pivotal concern is how "cities" are defined.
There have been attempts to construct cities from the bottom up (Makse et al., 1995) on the basis of geographical features of high-quality micro data. A recently introduced algorithm, the City Clustering Algorithm (CCA) (Rozenfeld et al., 2008), defines cities as maximally connected clusters of populated sites identified from gridded data at high resolution (Rozenfeld et al., 2011). In comparison with data collected from conventional "manual" sources such as the traditional population census and household interview, remote sensing data in the form of nighttime light are a more uniform, consistent and independent record for all world land areas (except for high latitudes) (Zhang and Seto, 2011). Therefore, it inspires us to take a step forward and systematically redefine cities from the bottom up in a much more natural and unified way with global composites of temporally stable nighttime light.
On the other hand, numerous city models have been designed to understand the self-organized nature of cities and to elucidate the origins of scaling laws. For example, some early models based on statistical physics, such as DLA (Diffusion Limited Aggregation) model (Batty et al., 1989) and correlated percolation model (Makse et al., 1995(Makse et al., , 1998, successfully manifested some scaling properties and fractal nature resembling those of real city systems. Recent models concentrate more on networks embedded in city systems, including infrastructure networks, street networks, social networks and so on. Confirmed by exponents measured in thousands of cities worldwide at different levels of development, the theoretic framework based on a hierarchical street network model (Bettencourt, 2013) turned out to be powerful in explaining both sub-linear and super-linear scaling relations. Although all these models are capable of evoking our deeper comprehension of the mechanism underlying city self-organization, no single model can reproduce all the significant patterns simultaneously so far. Therefore, there is still an urgent need for another (even simpler) model, which can generate as many as possible statistical patterns found in empirical evidence, to adequately capture the basic rules governing their formation.
Based on the high correlation between nighttime light and economic variables (for example, gross domestic product, urban population and so on) in a statistical framework for both timeseries and cross-sectional approaches, it is found that lightness can be used as a proxy for socioeconomic indicators on country and regional scales, especially for countries with low-quality statistical systems Sutton et al., 2007;Ghosh and Powell, 2010;Chen and Nordhaus, 2011;Henderson et al., 2012). It is necessary to point out here that even after careful calibration, nighttime light is never equivalent to or easily replaceable with GDP or population. In fact, as objective records of humanity's presence on the earth's surface , the distribution of nighttime light has also been used to develop proxies for urbanization (Sutton, 2003), armed conflicts (Agnew et al., 2008;Li et al., 2013), and, more straightforward, the spatial extent of light pollution (Butt, 2012;Bennie et al., 2014). For most places on the earth, stable nighttime light detected (given the preprocessing in "Methods" section) from outer space represents human activities. The brightness of artificial light increases, in general, as the human activities taking ARTICLE PALGRAVE COMMUNICATIONS | DOI: 10.1057DOI: 10. /palcomms.2015 place at corresponding sites get more intense. The intensity can both be driven up by the size of population group participating in the activity, and by the value of goods or services produced at that spot. In this sense, higher activity intensity generates more significant development, which is usually accompanied by anthropogenic landscape modification.
Unstructured data, which are predicted to account for 90% of all data created in the next decade (Gantz and Reinsel, 2011), provide us with relatively untapped sources of insights to define cities as systems of interactions. In this article, we take a step forward in fulfilling the potential of nighttime light data by defining natural cities "from the bottom up" with high-resolution satellite imagery (DMSP-OLS dataset). By this means, cities are defined through a completely data-driven approach with objective, spatially explicit and globally available empirical measurements of human development. Such a straightforward method is therefore supposed to be more economically meaningful than administrative or official definition of cities. We then systematically investigate the size distribution, allometric scaling and fractal dimension of natural cities. As the most panoramic and immediate embodiment of human activity and city growth (Frolking et al., 2013), the gridded imaging data allow us to compare the scaling exponents of natural cities across countries and continents. Furthermore, a dynamic model based on the simple geometric matching mechanism is built to simulate the pattern formation process identified in natural cities. Not only does our model successfully reproduce the allometric scaling and size distribution patterns, but also statistically generates comparable characteristic parameters as is in the empirical data. This is also the first attempt ever to model how nighttime light organizes itself into light clusters and grows over time.
There is no doubt that modern cities are combinations of topdown planning and bottom-up aggregation. However, the political aspect of planning and the political role of planners had not been emphasized until the 1960s and the early 1970s (Pissourios, 2014). In comparison with political discourses of urban intervention, bottom-up definition of city systems will help us get closer to the essence of human settlements. Our findings have profound implications for understanding the effects of simple aggregation behaviour in primitive stages of city formation and the urbanization process.

Methods
Material. The global satellite images of nighttime light used in our study are collected by the Operational Linescan System (OLS) of the US Air Force Defense Meteorological Satellite Program (DMSP) and archived at NOAA National Geophysical Data Center (NGDC). The flagship product of this dataset (Version 4) is the stable lights, which is available to download at http://ngdc.noaa.gov/eog/ dmsp/downloadV4composites.html. It is an annual cloud-free composite of the average digital brightness values for the detected nighttime light, altered to remove ephemeral lights and background noise . The image is 30-arcsecond gridded and spans from − 180 to 180 degrees longitude and from − 65 to 75 degrees latitude. The digital number (DN) values of the nighttime light range from 1 to 63, while 0 represents the identified and eliminated background noise and 255 represents an area where no cloud-free observation has been collected. In addition, although sunlit data, moonlight, glare, observations containing clouds and lighting features from the aurora are excluded from the DMSP nighttime stable lights dataset, gas flares are not. Therefore, we used the global gas flare map generated by NGDC  to identify and remove gas flares, to reduce the possibility of mistaking them for urbanized areas. This map can be downloaded at http://ngdc.noaa.gov/eog/interest/gas_flares_countries_shapefiles.html.
We obtained the administrative boundary data from GADM (the Database of Global Administrative Areas, available online at http://gadm.org/) to extract our regions of interest (ROI) from the global nighttime light images. Since we need to qualitatively study the evolution of global nighttime light in nearly two decades on a large geographic scale, it is necessary to inter-calibrate ) and re-project (Imhoff et al., 1997) each image correspondingly (Cauwels et al., 2014). In the examination of basic scaling properties of nighttime light clusters, we extracted regions with a vast territory and a huge population like China and Contiguous United States (CONUS) as our ROI. By using ArcGIS, we re-projected the nighttime light images of CONUS and China into Lambert conformal conic projection. For detailed comparison between nighttime light agglomerations detected and simulation results, we narrowed our scope down to part of the south central CONUS in which clusters including saturated lightness only makes up a negligible proportion. Due to the existence of upper limit in DN, there are some problems with respect to simple integration or summing of lights as a proxy measure of light emission from a given area. With the purpose of avoiding the problems of saturation in the DMSP images, a non-linear relationship between population and areal extent was used by Sutton to create a proxy measure of GDP (Sutton et al., 2007). For simplification, when to estimate the value of allometric scaling exponents, instead of referring to other data sources and indirect relationships, we only fitted the data points representing natural cities without saturated lightness. This can be validated by the fact that natural cities containing saturated pixels only account for a very small percentage.
Simulation. We use python 2.7 plus numpy and scipy packages to conduct all the simulations.

Results
Definition of natural cities. In order to define cities in a more natural and unified way, we explore the nighttime light images by varying the low light thresholds designating the spatial extent of contiguous development. That is, pixels whose digital number (DN) of lightness are below the threshold (DN th ) will be set as unlit, which means setting the digital number values of corresponding pixels to be 0 (Cauwels et al., 2014). We apply "burning algorithm" (Stauffer and Aharony, 1994) (see Supplementary Section 1 online) to identify nighttime light clusters representing natural cities, and require connected pixels to have lightness above the given threshold. In other words, for each stable light cluster identified, a natural city is defined at the same spot, whose area is the extent of land surface represented by all pixels in that cluster. Regardless of whether night lights depict population, economic activity, land cover modification or merely stable emitted light, natural cities are defined with explicit metric for anthropogenic development on global scales in this approach. Therefore, the adopted definition of natural cities are dependent on the threshold DN th , and higher DN th generally corresponds to more intense (or more developed) cutoffs while lower thresholds extend the lighted area to include less intensively developed periurban or suburban fringes and agricultural areas. We also suggest a specific threshold to define the natural cities according to the statistical properties which will be discussed in the following sections (see Fig. 1).
In this work, we set out to understand cities as a system of interactions, in other words, an aggregate of human activities. Viewing nighttime light as the embodiment of human activities (not necessarily human presence), the outline of maximal connected nighttime light clusters identified by clustering algorithm given specific threshold are defined to be the boundaries of natural cities. By using global composites of temporally stable nighttime light, we redefined cities from the bottom up by setting their contours in a more consistent and unified way. After delineating the geographic boundaries, the population and areas of corresponding cities can be calculated by adding up the population counts and land coverage within its boundary. Other socioeconomic indicators of each natural cities can also be calculated by the sum of related measures of socioeconomic activities that have happened within each boundary. Technically, gridded datasets of demographic and socioeconomic statistics are handy for such calculation.
Size distribution of natural cities. The distribution of city size follows Zipf 's Law (Zipf, 1949), that is, the size of the large cities decays with its descending rank r in a power law r − α with an exponent α. It has already been pointed out that the size (that is, area) of nighttime light clusters also obeys Zipf 's Law (Small et al., 2011;Jiang et al., 2014). Since the process from fitting the distribution to validating the hypothesis that the data is distributed according to a power law has been covered rigorously in the work of Clauset et al., we followed technical details of their maximum likelihood methods (Clauset et al., 2009) for preliminary test (see Supplementary Section 2 and Supplementary  Fig. S11 online). Furthermore, we find in this study that only the upper tail of size distribution curve can be well described by a power law, whereas the whole spectrum can be better depicted by the discrete generalized beta distribution (DGBD) function (Martínez-Mekler et al., 2009) (see Fig. 2).
in which r max is the total number of natural cities in the region under examination, r is the rank, and S(r) is the size of the r − th largest connected cluster measured by the area it covers as is in the image after preprocessing (more details in Methods). A, b, and α are parameters to be estimated, in which exponent b characterizes the exponential head. The conventional Zipf 's Law (with a power law exponent α) can be recovered by setting b = 0 (see Supplementary Section 2 online). In comparison with only one exponent at which city size decays linearly in the log of rank, estimation on the basis of DGBD provides us with an opportunity to scrutinize and evaluate our model across the whole spectrum of size distribution. In fact, we take into account the shape of the head of the size distribution curve in our model design. After all, if the only testable prediction a model makes is about the shape of the tail, it can be intrinsically hard to check its design. China is selected as the representative country to illustrate the size distribution of natural cities, shown in Fig. 2(a) (see Supplementary Section 4 online). As the low light threshold is increased from dim to bright, the overall pattern changes for two reasons: (1) the spatial extent (area) of each spatially contiguous cluster on developed land shrinks, and (2) large clusters connected by dimmer regions fragment as dimmest lights are exceeded by the threshold (Small et al., 2011). This results in increasing the curvature of the rank-size distributions and decreasing the dominance of the largest cities.
By varying the threshold value DN th , we systematically investigate the size distributions of the United States and China (see Fig. 3  transition of the power law exponent with decreasing DN th is analogous to the phase transition phenomenon observed in percolation models (Stanley et al., 1999). In fact, anthropogenic land cover development, by its nature, is consistent with the notion of a spatial phase transition from more isolated to more interconnected development patterns (Small et al., 2011). Thereafter, different from the "natural cities" defined by Jiang et al. according to "head/tail breaks" (Jiang et al., 2014), we are especially interested in the definition with DN th set to be 4.
Allometric scaling of natural cities. Previous studies have proved that many socioeconomic variables, such as GDP, number of patents and total salary, can be predicted by the size of a city which is measured by the population (Bettencourt et al., 2007). In natural cities specified by nighttime light clusters, we study the most important allometry between total nighttime lightness which is the embodiment of human activities and city size quantified by the areas of corresponding clusters. For a region we considered, such as Contiguous United States (CONUS) and China, let L i denote the total nighttime lightness in natural city i, and S i the size (that is, total area) of it, then we suggest that the following scaling form holds: Parameter β is the allometric exponent which characterizes the increment speed of lightness relative to spatial expansion. Therefore, larger β indicates more remarkable scale effect of the natural cities in regarding regions. The allometric exponent is relatively insensitive to the threshold compared to the Zipf 's exponent. We select the United States (CONUS region) as an example to illustrate the allometric relationship in Fig. 2(b). Exponents across different thresholds are largely concentrated in the interval [1.11 − 1.14], fairly close to the allometric scaling between economic activity and city population previously demonstrated in extensive research (Bettencourt et al., 2007). By taking the average of the allometric scaling from 1992 to 2009 under each threshold (DN th ), it is shown that the exponents of both the United States and China generally decrease as DN th varies from dim to bright while remain above 1 (see Fig. 3). It is interesting to observe from Fig. 3(b) that at the low light threshold DN th = 4, both the United States and China attain almost the same value which, more importantly, also approximates their inflection points.
The fractal nature of natural cities. Cities are classic examples of fractals in that their form reflects a statistical self-similarity or hierarchy of clusters (Batty and Longley, 1994). We use the standard box-counting method to calculate all the fractal dimensions of natural cities larger than 60 × 60 km 2 in the North America (for top 40 largest cities in the United States, see Supplementary  Table S1 and Fig. S6-10 online), and plot against the cluster size in Fig. 4(a). The trend that the fractal dimensions increase with cluster size is quite clear in Fig. 4(b). Albeit not influential on the dimension computed, saturation effect is visualized here by distinguishing clusters including no pixels at the upper limit of lightness with grey marker in Fig. 4

(a).
Comparison with MSA and CCA. It is a long-standing problem that different definitions of spatial units based on administrative boundaries give rise to inconsistent conclusions at different scales (Rozenfeld et al., 2008). In comparison with our definition of natural cities, administrative boundaries of MSAs provided by the US Census Bureau (available online at http://www.census.gov/ geo/maps-data/data/cbf/cbf_msa.html) are constructed manually based on the subjective judgement. Since they are only available for the most populated cities in the United States, any study of human agglomerations based on MSA is limited within the subset at the upper tail of size distribution. Moreover, due to the fact that many MSAs are constituted by aggregating small disconnected clusters, this definition framework overestimates the area of small agglomerations (see Fig. 5 (a)) by including large unoccupied areas (Oliveira et al., 2014). As a consequence, super-linear scaling exponents associated with interactions are supposed to be smaller for MSAs than natural cities.
While MSA may delineate similar areas as natural cities (given proper DN th ) for large and developed agglomerations (see Fig. 5 (c)), its origin for statistical purpose makes MSAs left behind by the real dynamics of city systems that is taking place. Natural cities that come from the mergence of previously separate agglomerations tend to be designated into different MSAs as before (see Fig. 5 (b)). Recently, as a revision of conventional MSA, a two-step methodology is developed which allows for a dynamic, economically meaningful definition of cities (Arcaute et al., 2015). This revision greatly improves the consistency of definition by overcoming constraints of historical dependence and by supporting implementation in almost all countries. However, since the unit of agglomeration for this algorithm must be division or district already defined for administrative or political purposes, the realizations of such definition actually vary among different countries.
Methodologically, the definition of natural cities is more similar to CCA than MSA. The uniqueness of the idea behind natural cities lies in our emphasis on human activities over human presence. To demonstrate how the definition based on nighttime light will affect the allometric exponents, we examined the pivotal allometric relationship between area and population (see Supplementary Section 12 and Fig. S13 online). Specifically, given a relatively low DN th (for example, DN th = 4), natural cities are neither as overestimated as conventional MSAs, nor as actually populated as cities defined by CCA. Places significantly influenced by human activities yet not highly populated are taken into account as parts of natural cities, which can explain possible differences in scalings compared with CCA.

Modelling the pattern formation in nighttime light clusters.
Most human activities scale super-linearly with city size, for the number of potential face-to-face interactions increases more than proportionately (Batty, 2011). In essence, our definition of natural cities based on nighttime light clusters views city systems as sets of interactions which enable construction, production and creation. We build a model from this viewpoint to reproduce the patterns recognized in the empirical nighttime light clusters. Consider a W × W lattice as the correspondence of the interested region. At the beginning, this is an empty place. At each time step t, this region can either (1) generate a new cluster (city) with a probability ϵ; or (2) expand one of the existing cluster (city) with the probability (1 − ϵ).
If rule (1) is activated, a random position X ∈ {(x, y)|xoW and yoW} is selected as the place where a new cluster forms. Then m (uniformly distributed in [1, m 0 ]) new settlements represented by points are generated around X in a distance less than r to become the first residents of this new city. If rule (2) is activated, then a new settlement p is generated and will randomly select an available place X in the region to settle down. Here, X is available if there is at least one existing settlement within the r radius of X (see Fig. 6). We call this rule as the geometric matching mechanism which is more important than other rules. In such scenario, a new point can settle down at any available place no matter how many existing settlements are already there.
However, due to the basic needs of living materials and the scarcity of natural resources, densely populated area is unattractive, even not liveable, for newcomers. In real life, the crowding effect should also be taken into consideration. Assuming that i can settle down at the place X with a probability Π(μ(X)) = μ(X) − γ which depends on the local density of settlements, in which μ(X) is the total number of existing settlements within the r radius of X. γ is a non-negative free parameter called crowding coefficient in our model. Therefore, if the place X is very crowded, it is highly probable that i may fail to settle down there.
In this way, new settlements join this region sequentially. At each time step, we recognize clusters by identifying all connected settlements as one unique cluster. Here, settlements i and j are connected if their distance is less than r. According to (2), the cluster will naturally be expanded through random selection, while larger clusters may have higher probability to be selected since they can provide more available positions. When a cluster expands, it may collide with another one and spontaneously merge into a new huge cluster. To simulate the formation of natural cities in our ROI, the model can technically generate settlements according to the following steps: Step 1: Generate a random number ξ; Step 2: If ξoϵ, then pick up any place X in W × W region to generate m new settlements around X (build up a new cluster); Step 3: Otherwise Step 3.1: Pick up a random place X in all the available places (close to existing settlements); Step 3.2: Generate a random number ψ, and if ψoμ(X) − γ , then generate a new settlement at X (expand an existing cluster).
Step 4: Identifying all clusters of settlements.
Step 5: Go to Step 1 until the total number of clusters in the region is equals to the number of clusters in the empirical data.
We suppose that each settlement can interact with its connected neighbours within r radius, and the intensity of a b socioeconomic interaction (equal to the number of its neighbours) is proportional to the lightness at that place. Then, the total lightness of cluster s is calculated as the following equation: in which C s denotes the set of settlements in the cluster s, O j is the set of all neighbours of settlement j within r radius. Similarly, we define A s the size of the cluster s as the total number of lattices covered by these settlements. One lattice is covered by cluster s if there is at least one settlement belonging to s in this lattice. Therefore, we can study the allometric scaling relationship between L s and A s as well as the distribution of A s . To sum up, there are five parameters in the model: W, r, ϵ, γ, m 0 , only three of which (ϵ, γ, m 0 ) are free. W = 1000 is set to the size of the region considered, and r is fixed to be 1. γ is interpreted as the sensitiveness toward the crowdedness of the region. It can adjust the allometric scaling exponent between A s and L s . ϵ models the relative rate of frequency between expanding a city and building a new one. This parameter can take effect on the slope of the size (that is, area) distribution curve. As the maximum settlements to be built up in a new cluster, m 0 can influence the exponential head of the size distribution curve.
Actually, when we set ϵ = 0 and m 0 = 1, the model is analytic solvable (see Supplementary Section 7 online). The allometric scaling exponent is: Therefore, the scaling exponent varies in the interval [1, 2] which is consistent with our empirical observations. When W → ∞, ϵ ≈ 0 but is non-zero, and m 0 is small, the growth process of clusters can be understood as a sub-linear preferential attachment (Yule-Simon) process (Krapivsky et al., 2000) because the clusters are approximately independent each other  when W → ∞, and the growth of an existing cluster is proportional to its size. Therefore, the size distribution of clusters can be solved as the following formula when t → ∞ (see Supplementary Section 8 online), where P(A)dA is the probability of a cluster with size in [A, A+dA]. c is the constant coefficient of the scaling relation between A s and N s , the total number of settlements in cluster s. η = (2+2γ)/(3+2γ)o1 is the exponent of the same scaling relation. ζ is the solution of the equation . This distribution resembles a power law if η is close to 1, and the exponent is not analytic solvable (Krapivsky et al., 2000).
Different from observed dynamics of clusters identified with CCA from 1981 to 1991, which exhibit all five types of changes (that is, no change, expansion, reduction, division and merge) in cluster shape (Rozenfeld et al., 2008), our model exhibits only three types (that is, no change, expansion and merge) in computer simulations. Since clusters may collide each other and merge together over time, they are generally not independent. Therefore, the size distribution and allometric scaling exponents may not follow previous analytic values exactly. We implemented a series of numeric experiments to generate 921 clusters, namely a multi-city system containing 921 natural cities, in comparison with a nighttime light satellite image containing the same number of natural cities identified (denote the corresponding area as our region of interest, or ROI for short; see detailed information of ROI in Supplementary Section 3 online). By tuning γ to 1.5 and ϵ to 0.03, we can obtain simultaneously the exponents of the allometry and Zipf 's distribution comparable with, if not exactly equivalent to, those from natural cities within the multi-city system in our region of interest (see Fig. 7(c) and (d)).
Each cluster, namely every single natural city in the multi-city system generated on the basis of simple rules of interaction and land development, also manifest self-similarity across scales. Thus the extent to which they fill space can be measured by their fractal dimension. Instead of computing fractal dimension for each cluster, we study the largest cluster (in most cases, this is the cluster that first comes into being in simulation) in different time steps. Such tracking evaluation is meaningful since the growth processes of all other clusters can just be viewed as duplication of the course through which the largest cluster expands. We also find a stable positive correlation between its size and fractal dimension (see Fig. 4(b)) in general. Since all clusters expand following the same rule, our model predicts such a positive correlation between area and fractal dimension of clusters.  PALGRAVE COMMUNICATIONS | DOI: 10.1057/palcomms.2015 This prediction is validated by empirical data on regional, country and continental (see Fig. 4(a)) scale. However, it is also noteworthy that the fractal dimensions from simulation are in general larger than that from empirical results. As can be directly perceived through Fig. 7(a) and (b), clusters generated by our model are more compact in comparison with real natural cities.

ARTICLE
In short, our model exploits the simple mechanisms to reproduce all the observed empirical patterns of natural cities. Among all the rules, the geometric matching (rule 2) is the most crucial one-it indicates that the connection to and thus interaction with neighbours is of fundamental importance.

Discussion
Previous sections present the size distribution, allometric scaling, and fractal geometry of natural cities in two countries. However, the regions we focus on in this study of natural cities is not limited by the country border lines, but can be extended to continental scale. It is interesting to make a comparison of the exponents among six continents (Africa, Asia, Australia, Europe, South America and North America) when DN th = 4 (detailed comparison under all low light thresholds in Supplementary  Fig. S5 online). The characterizing parameters of distribution α listed in Table 1 show that the size of natural cities in Europe and Asia are relatively more heterogeneous, in other words the differences between big cities and small ones are quite significant there. By contrast, two continents in the Southern Hemisphere have the smallest exponents indicating that the size variances of natural cities there are relatively small. As regards to allometric scaling, Australia has an extraordinarily large exponent while the exponents of all other continents are approximately 1.13. It is thus strikingly illustrated that the socioeconomic scale effect is comparatively higher in Australia. In fact, Oceania regions (for example, Australia/New Zealand) is the most urbanized part of the world, with 90% of its population urbanized in as early as 2000, while barely 24% were urbanized in Melanesia which is the least urbanized region at that time (United Nations Population Division, 2002).
In this article, we exert ourselves to exploit the potential of nighttime light for a better understanding of city systems. The correlation between the spatial distribution of nighttime light and economic activity has been utilized in several studies for estimating economic activity at the national and sub-national levels (Doll et al., 2000;Elvidge et al., 2001;Doll et al., 2006;Sutton et al., 2007;Henderson et al., 2012;Zhou et al., 2014). The underlying expectation, despite the researchers' ultimate goals in such literature (for example, to predict or to evaluate), is largely to improve resolution and accuracy of census statistics effectively. To this end, nighttime light functions well as worldwide and cross-calibrated supplementary information, which is free of any possible interference from political interests or limitation of national statistical system quality. However, what deserves more attention is the new perspective on city systems provided by this dataset, which goes beyond constraints from administrative or official definition of urban area. And this may bridge the gap between empirical probes into cities and mechanism underpinning emergence and growth of cities in the real world. In essence, process of land cover development occurs by combination of isolated nucleation and intensification of previously developed land in human history. Overall growth of city system just follows a progression analogous to the spatial growth that results from lowering light thresholds in the definition of natural cities.
But it is also noteworthy that DMSP nighttime light data is far from perfect in terms of depicting city growth, mainly due to the problem of overglow and saturation. Although nighttime light brightness and population density are highly correlated, the relationship is less direct for dimmest lights detected. By setting proper threshold, we actually determine the standard of "urbanization" based on intensity of human activities instead of population counts. Therefore, the low levels of light associated with the peripheral overglow of larger settlement is not a problem to our final end. Since these overglowing regions generally still encompass areas with anthropogenic modification, nighttime light can act as a suitable proxy for intensity of human activity or level of socioeconomic development.
There are reasonable concerns that decreases of brightness (especially regional decreases in Europe) may not be caused by decline in human activity intensity there, but by increases in efficiency of municipal lighting. Belgium, for example, where almost 100% illumination of the country's highways was indeed able to be seen from space with a telescopic lens, switched off lighting in the central reservations of many motorways by early 2014. Despite widespread increases in the brightness of nighttime light across Europe over the last 15 years, there also exist regionally significant decreases attributable to economic or industrial decline (Bennie et al., 2014), as has been noted for some countries of the former Soviet Union and Eastern Europe (Elvidge et al., 2005). Improved environmental awareness can lead to cost saving, but more prosaically, it is usually inspired by the need of cost reduction in the first place. Though subtle, the links between reductions of light emissions to both declines in industry and to increased municipal lighting efficiency are still possible to be disentangled given careful calibration based on case studies.
This article proposes a model to simulate the growth of natural cities and reproduce all observed properties with very encouraging results. Except for the fractal dimension, all other scaling exponents observed can be precisely generated by the model simultaneously. This is the first time ever empirical allometric scaling, size distribution, and fractal geometry have been obtained from one single model. In this model, the geometric matching mechanism is the most important one compared with other rules. However, since the modern transportation tools have actually twisted the geographical space of cities, the geometric matching in the Euclidean space is not exact. It also explains why our model generates clusters which are much more compact than real natural cities. Therefore, the geometric matching mechanism in a twisted space by taking into consideration the influence of ships, planes, and land vehicles et al. will later be explored for future studies.
Our model is no doubt a highly simplified and abstract sketch of real cities. Many factors, such as economic, social, and environmental ingredients are not considered. However, the simple design in the current model enables us to introduce Notes: b and α are parameters estimated from DGBD regression of the size distributions. β is the allometric exponent between lightness and area. All values are estimated by setting DN th = 4 in 2009. more realistic but complicated rules and parameters. We believe that the extended models can ultimately be applied to city planning in the future. However, the success of this model at least proves that natural cities, on a macro scale, are just like other complex systems following some ubiquitous selforganization rules and can be adequately captured by simple yet elegant models. In this regard, our work has made its own contribution to the development of urban theory.