Intra-Urban Scaling Properties Examined by Automatically Extracted City Hotspots from Street Data and Nighttime Light Imagery

: A country can be well-comprehended through its core cities. Similarly, we can learn about a city from its hotspots, as they manifest the concentration of urban infrastructures and human activities. Following this philosophy, this paper studies the intra-urban form and function from a complexity science perspective by exploring the power law distribution of hotspot sizes and related socio-economic attributes. To detect hotspots, we rely on spatial clustering of geospatial big data sets, including street data from OpenStreetMap platform and nighttime light (NTL) data from the visible infrared imaging radiometer suite (VIIRS) imagery. Unlike conventional spatial units, which are imposed by governments or authorities (such as census block), the delineation of hotspots is done in a totally bottom-up manner and, more importantly, can help us examine precisely the scaling pattern of urban morphological and functional aspects. This results in two types of urban hotspots—street-based and NTL-based hotspots—being generated across 20 major cities in China. We ﬁnd that Zipf’s law of hotspot sizes (both types) holds remarkably well for each city, as do the city-size distributions at the country level, indicating a statistically self-similar structure of geographic space. We further ﬁnd that the urban scaling law can be effectively detected when using NTL-based hotspots as basic units. Furthermore, the comparison between two types of hotspots enables us to gain in-depth insights of urban planning and urban economic development.


Introduction
As a result of urbanization or the continuous influx of people into cities, the number of worldwide urbanites is predicted to be 6.9 billion by 2050, accounting for 68% of the world's population [1]. The urbanization in China has been unprecedentedly rapid as well in the past few decades [2], reaching 60.6% nationally in 2019 [3]. Consequently, the grasp of city form and function-that is, how cities look and work-has become the key to our sustainable development. Given the circumstances, city-related research has attracted scientists from a variety of subjects and has, inevitably, become cross-disciplinary, including geography, economics, computer science, and physics, etc. To converge these disciplines, scholars have called for a new science of cities in the past few decades, in which they view cities as an organized complexity [4], for studying cities' fractal shapes, complex structures, and nonlinear dynamics (e.g., [5][6][7][8][9][10][11]).
One major aspect of urban complexity is its underlying scaling properties. The scaling pattern of urban entities can be categorized into two perspectives: The power law distribution of a single quantity, such as city sizes (Zipf's law [12]), building heights [13], street lengths [14], and leisure venue densities [15], and the power relationship between two quantities, such as populations versus innovations ( [16,17]) or gross domestic product (GDP) versus street fractality ( [18,19]). This study uses the terms scaling and power law interchangeably. Urban scaling is, to a great extent, a ubiquitous pattern across different measures. Moreover, the theory developed by Bettencourt et al. [16], which is behind the power relationship between urban populations and other socio-economic measures, has been formulated as fundamental laws about cities: Universal scaling law. However, recent studies have shown that the universal scaling law may not work as expected, as the scaling exponent is sensitive to different city boundaries or ineffective urban areas [20,21]. This controversy is likely to be bound with the top-down methods of defining geographic units by governments and authorities, such as administrative city boundaries, census tracts, and some equally partitioned cells, which are essentially for management purposes and hardly consider the scaling pattern of urban morphological and functional entities.
The arrival of geospatial big data has triggered a new paradigm for urban analysis since geospatial big data, such as remote sensing (RS) images and location-based social media data, has the capacity to offer fine-grained, massive-scale geographic information [22]. For instance, nighttime lights (NTL) data, also referred to as RS of human beings and their activities [23], are globally downloadable and can manifest the development of urban and regional areas. OpenStreetMap (OSM), a pioneering volunteered geospatial information platform, provides street data across the globe for probably the first time in human history [24]. Both NTL and OSM data help researchers construct alternative modeling units for spatial analyses at both intercity and intracity levels, and remove the barriers of inter-regional incomparability. The most recent relevant studies are so-called natural cities, referring to the objectively defined cities based on different types of urban elements from the open data, such as building footprints, street nodes, and points of interest (e.g., [25][26][27][28][29]). However, most of these studies take the derived cities as a whole to understand the scaling structure over a region or country, but seldom calibrate a "local" understanding of such spatial configuration at the intracity level.
Thus, the present study attempts to investigate the intra-urban scaling properties through the lens of city hotspots. A city is formed by highly concentrated areas of human settlements or activities within a country extent [30]. Likewise, if we scale down our scope from a country to one of its cities, such concentrations can be regarded as urban hotspots.
With the advance of geographic information system (GIS) technologies, urban hotspots can be delineated more precisely on the support of geospatial big data and bottom-up approaches. The study contributes to the current literature in three aspects. Firstly, we followed the ideas of previous city delineation methods to derive two types of urban hotspots across 20 Chinese cities: Street-based and NTL-based hotspots, from respectively the spatial clustering of individual street nodes and NTL image pixels with the cutoff determined by data's inherent scaling properties (see details in Section 2.2). Secondly, we found that Zipf's law held remarkably well for both street-based and NTL-based hotspot sizes per city, as do the city-size distributions on the national scale. The scaling exponents derived based on NTL-based hotspots were also consistent with the established regimes, implying that NTL-based hotspots can act as better spatial units for urban analysis. Thirdly, we found that the spatial discrepancy between the street-based and NTL-based hotspots can lead us to deep insights on urban planning and development.
The remainder of this paper is organized as follows. Section 2 introduces the data sets and the designed methods for urban hotspot delineation and related scaling analyses. Section 3 presents the maps of the detected hotspots across the top 20 cities in China, as well as the power law metrics of hotspot sizes and associated socio-economic attributes. Section 4 further discusses the intra-urban scaling properties. Section 5 concludes the study and points to future research directions.

Data and Data Processing
We selected 20 well-developed cities in China as study areas and primarily made use of the following three data sets: (1) VIIRS imagery, (2) OSM street network, and (3) socio-economic grid data ( Figure 1a). All data sets are national coverage. The NTL data was obtained from NOAA/NCEI [31]. We chose one monthly image at June 2020, of which the resolution is 15-arc-s (about 500 m at the Equator). We reprojected and cleaned the image to get rid of noises (lit spots) such as burning wildfires and oil drilling, based on the method proposed by Elvidge et al. [32]. The national street network was downloaded from OSM, including 4,419,603 segments from which we extracted 3,172,001 street nodes based on the criterion that a node must intersect with three segments. The socio-economic grid data include the GDP and population from the National Resources and Environment Database of the Chinese Academy of Sciences [33] and environmental grid data include CO 2 emissions from the National Earth System Science Data Center [34]. Raster data sets for GDP, population, and CO 2 were collected in 2010 and had a 1 km resolution. To perform the analysis, we clipped out both the vector and raster data using each of 20 city administrative boundaries, then conducted zonal statistics of cells with socio-economic attributes for each city, which were further joined with city hotspots (Figure 1b).

Urban Hotspot Detection
We adopted the spatial clustering method for urban hotspot calculation and delimitation. As there were two types of data sets (street junction nodes and NTL pixels) to be processed, we applied two rules for cluster detection of each data set: Point proximity and lit pixel adjacency. The threshold (distance between points or pixel value) for clustering was determined by the data's inherent scaling properties uncovered by head/tail breaks and power law detection methods.

Spatial Clustering of Street Nodes and NTL Pixels
Urban hotspots-that is, populated areas in a city-are the basic unit for the analysis in this study. Traditional urban analysis uses pre-defined administrative units provided by local authorities or grids with different resolutions. However, both spatial units cannot represent the merit of "concentration" as they are defined either from a top-down or arbitrary manner. To overcome this issue, we adopted a spatial clustering approach to objectively delimit the boundary of a hotspot from the dense areas of street junctions or lit pixels.
We chose two clustering approaches for each data set. For street junctions, we first computed the triangulated irregular network (TIN) to get junction-junction proximities.
As Figure 2a-c shows, the area of urban hotspot can be directly obtained through the conversion of short TIN edges between points. For NTL images, the first step is to vectorize each raster pixel into a polygonal feature with the light value maintained (Figure 2e), then the hotspot can be derived through grouping the adjacent lit pixels (Figure 2f).
The above procedures can be simply done using any mainstream GIS or RS image processing software (such as ArcGIS and Erdas). The major difficulty lies in identifying the cutoff value for the classification short/long edges and dim/lit pixels across a set of urban areas. In other words, it lacks an objective criterion to make the linkage between the morphological hotspot (the concentration of urban infrastructure) and a set of proximate street junctions or the functional hotspot (the concentration of human activities) and a group of lit pixels. The same issue occurs when delineating the city boundaries (regardless of administrative boundaries) at the country or cross-country level, whereas prior studies (e.g., [35]) have made use of the universal scaling property for finding the effective cutoff value. In a similar spirit, the next section will introduce how to obtain the optimal cutoff value for the accurate delimitation of urban hotspots.

Scaling Analytics for Identifying the Cutoff for Spatial Clustering
A vast body of literature has investigated city-size distributions in different countries. Most of those studies have used the power law model to characterize the uneven spatial distribution of cities, as well as their sizes, such as Zipf's law [12]. Zipf's law states that there is an inverse relationship between the rank and the size of a city. In other words, the largest city is twice as big as the second largest city, etc. Such a statistical distribution would strikingly present the long-tail effect or scaling pattern of far more small cities than large ones. In most cases, the scaling pattern recurs within the power-law distribution and leads to an inherent hierarchy, which can be derived through the head/tail breaks classification scheme. In this study, we change our perspective from a "country-to-city" relationship to "city-to-hotspot" one. In this way, we can borrow the scaling analysis methods (power law detection and head/tail breaks), which were previously used for finding the cutoff value of city demarcation, to delineate hotspots. To start with, we shall first introduce briefly Zipf's Law, power law, and head/tail breaks.
Referring to the size n of each city relative to its rank number r, Zipf's law is denoted by Equation (1): where b usually is equal to 1, indicating that the city size is equal to the reciprocal of its rank. Another way to describe Zipf's law is the Pareto distribution (or power-law, which is a derivative of Pareto distribution) [36]. To do this, it is equivalent to use the inverse function of Equation (1) as r ∼ n − 1 b , where r is further treated as the proportion, Pr, to the whole population by the cumulative distribution function (CDF), and it is relative to how many of the cities are greater than the size, x, is defined as follows: where k > 0. For a specific point of x, the power-law is acquired by the derivative of Pareto distribution by the probability density function (PDF) as: where C is a constant and α = k + 1. In practical terms, the power-law distribution could only be discovered in one part of the whole dataset, where there must be some lower bound denoted as x min . A formal form of the power-law is given as follows proposed by Clauset et al. [37]: With the fixed lower bound x min , the power law exponent α is then derived from the robust maximum likelihood estimation (MLE) method, noted as Equation (5): So far, we can remark that, for detecting Zipf's law, the power law exponent should be two rather than one. Furthermore, a modified Kolmogrov-Smirnov test [37,38], needs to be performed to determine the extent of fitness for the data to an ideal power-law fitted model using the derived x min and α values. Every time we generate 1000 synthetic datasets that follow a perfect power law above x min but have the same non-power-law distribution as the original dataset. Then, we check how many times the maximum difference between each synthetic data and the fitted model are larger than the one between the original dataset and the fitted model, the ratio of number of times to 1000 is the goodness-of-fit index p-value. We set p-value ≥ 0.05 as the acceptance of data being a power law in this study, meaning that at least 50 among the 1000 synthetic datasets are less "power-law-distributed" than the original dataset.
Zipf's law can be used as an effective assessment when performing city demarcations. In other words, if the demarcated city sizes follow Zipf's law, we think that the result is valid. The question then narrows down to how to derive cities whose sizes follow Zipf's law from geospatial datasets, such as the TIN model and NTL imagery ( Figure 3). Here, we introduce the head/tail breaks method [39] to effectively locate the cutoff value. Put simply, data with a power law distribution can be divided into a high percentage in the tail (≥60%) and a low percentage in the head (≤40%) at the arithmetic mean. Therefore, for TIN and image models, the head refers to long TIN edges and light pixels, and the tail refers to short edges and dark pixels. The process then runs recursively for the head part until the head percentage is no longer small (say, ≥40%). During the process, a series of arithmetic means were iteratively computed, naturally forming a scaling hierarchy of the data. The number of mean values, also known as the ht-index [40], can then characterize the tendency of data being power-law-distributed. Namely, the larger the ht-index value, the more likely it is that the data is a power-law. Prior studies have used these nested mean values as cutoffs for extracting the so-called natural cities whose sizes obey Zipf's law at either national or cross-national levels (e.g., [41]). However, the use of those values for hotspot derivation at the city level remains under-researched. The present study would detect urban hotspots through a combination of head/tail breaks for locating the feasible cutoff and MLE method for examining Zipf's law.

Power Function Fitting for Intra-Urban Scaling Law Examination
The examination of urban scaling concerns two perspectives: The power law detection of a single urban indicator (as mentioned in Section 2.2.2) and the power relationship between two types of urban quantities (for example, urban areas versus populations). The latter have been formulated as the universal scaling law [16] for most of the urban indicators, which uses the power function fitting between an urban indicator and the urban population size across cities at time t, denoted as Equation (6): where β is the scaling exponent and k is the constant. The scaling exponent β can be further investigated by means of three categories: The sub-linear (β < 1), linear (β ≈ 1), and super-linear (β > 1) scaling relationships between urban measures [16]. To elaborate, for β < 1, it normally refers to the need of a city's infrastructure scales sub-linearly with its population size due to the economies of scale, whereas the number of a city's innovations and crimes scales super-linearly (β > 1) due to the endogenous social interactions. The regime of β ≈ 1 describes the pattern that the individual demands in a city is proportionate to the urban population size. In this study, we use the detected hotspots as alternative spatial units to reexamine the urban scaling law. To do so, we conduct the power function fitting between urban socio-economic metrics (such as population, GDP, and CO 2 emissions) that are within urban hotspots. To compute the scaling exponent, we first take the logarithms on both axes and adopt the ordinary least-squares linear regression for fitting. The scaling exponent is then the slope of the fitting line.

Derived Urban Hotspots in the Top 20 Chinese Cities
We applied the urban hotspot detection method on street nodes and NTL imagery, respectively, across top 20 Chinese cities, ranked by GDP. To derive the hotspots from the street nodes, we established big TIN models for each city, whose TIN edges range from tens to hundreds of thousands (Table 1). The heavy-tailed distribution statistics were striking for each TIN model, as the average edge length (the mean length of l edge is about 450 m) was classified effectively between short and long TIN edges according to their imbalanced ratios (around 80% versus 20%). The observation of 80/20 division, namely the scaling pattern of far more short TIN edges than long ones, objectively reveals the uneven spatial distribution of street node densities. The delineation of urban hotspots for each city was then conducted by grouping and converting those short edges into many different-sized hotspots. The area of resulting hotspots per city followed well with Zipf's law, as the mean value of 20 cities' power-law exponents was 2.01 (for more details of the basic statistics and related power-law metrics of hotspot size, see Section 3.2). Figure 3 presents the appearance of hotspots across selected cities, clearly showing that a few largest patches were located in the downtown and numerous smaller ones were spaced dispersedly in places other than the city center.  The urban hotspot extraction from NTL data went through experiments with a series of "candidate" mean values along with the head/tail breaking process on each image. To start with, the number of pixels for each image ranged widely, from 9397 (Shenzhen) to 353,344 (Harbin) and, interestingly, also followed the fat-tailed distribution. More specifically, among 20 city NTL images, most of the images (14) contain fewer than 78,045 pixels, some (five) between 78,045 and 141,623 pixels, while only one image has more than 141,623 pixels, resulting in a ht-index value of 3, meaning that there are three hierarchical levels of images regarding the number of pixels. Moreover, the ht-index for the pixel values of each city image was even higher. Figure 4 shows clearly that each image contains far more dark pixels than light ones, and such a scaling pattern recurs at least five times, indicating that there were no fewer than five average lightness values of each image achieved as candidate threshold values for a single city's hotspot delineation (see Appendix A for more details of the head/tail breaks method applied to the pixel values of each city's NTL data). Therefore, for every image we merged the vectorized pixels whose values above each derived candidate thresholds based on head/tail breaks to extract the urban hotspots, ensued with power law detection for each set of the hotspot results. The summary of statistical results for varying thresholds is presented in Table 2, which shows that the optimal cutoff value resided in the third level, since its power-law exponent was closest to 2, leading to hotspots being most akin to the Zipf's law configuration. It should be noted that the average of the cutoff values across 20 cities (33.086) largely echoes the optimal threshold (33.14) based on the VIIRS NTL data in 2013 for Chinese city demarcation [35]. Following the located cutoffs for each image, the layout of extracted urban hotspots exhibited a picture that was overall similar to that from street nodes in terms of the imbalanced spatial distribution from city center to periphery (Figure 4). By comparing Figure 3 with Figure 4, it is clear that two types of patches overlapped, but in varying degrees, with each other, indicating there were similarities and differences between urban physical and functional extents. Here, we applied the intersection over union (IoU) metric to compute the overlapping ratio between two types of hotspots for each city, the average ratio for 20 cities was around 0.27 (see more details in Appendix B). It appeared that inland cities were inclined to have larger ratios, such as Shenyang, Xian, and Zhengzhou had most overlays (around 0.4), whereas coastal cities such as Shenzhen and Qingdao held much less (e.g., only 0.11 for Qingdao). We further opted to map the overlay between two types of hotspots among the top four representative cities in China: Beijing, Shanghai, Guangzhou, and Shenzhen ( Figure 5), whose IoU metrics are all smaller than the average, i.e., 0.26, 0.17, 0.21, 0.18, respectively. Moreover, it is intriguing to note that detailed disparities can be found with respect to the extent of dispersive patches. In other words, with similar power-law exponents (around 2), the sizes of NTL hotspots in top cities seemed to be more even and the spatial distribution were more dispersed than those of street hotspots.

Intra-Urban Scaling Properties Based on Derived Urban Hotspots
We applied the robust power law detection based on the MLE method to two types of hotspots in 20 cities. For each city, we listed the power-law fitting metrics regarding its hotspot areas detected using the cutoffs derived from head/tail breaks (Table 3). We can see that Zipf's law held remarkably well for both types of urban hotspots. As stated, the power-law exponents for street hotspots were centered at 2.01 ± 0.15, while the averaged exponent value for NTL hotspots was slightly smaller, 1.921 ± 0.19, due to the exception of Chengdu (1.46). Most of the p-values were above 0.05 and readers can cross-check the results in Table 3. In addition to the hotspot sizes, we also examined the power law fit of the socio-economic status within the hotspots in the top four cities. As Figure 6 shows, the power-law distribution still holds for GDP, population, and the amount of CO 2 emissions per hotspot, respectively. However, the values of exponents for each city performed slightly differently. Specifically, the exponents of three urban metrics inside the hotspots remained relatively stable with the hotspot size in Guangzhou and Shanghai, but less so in Beijing (up-and-downs around α Area ) and Shenzhen (all smaller than α Area ).  We further investigated how these extracted hotspots worked as cores of each city. Ideally, there should be a disproportionate relationship between hotspot areas and the amount of pertained resources. Consequently, 3% of the city area, constituting either type of hotspot, accommodates, on average, around 15% of GDP, 25% of population, and 20% of CO 2 emissions (Table 4). Extreme cases such as Shenyang, Wuhan, and Kunming showed that derived hotspots could even account for more than 40% of the city's total population or GDP. Such imbalanced ratios enabled us to make use of those urban indicators within the hotspots for exploring the intra-urban scaling law. After correlating the total areas, GDP, and CO 2 emissions with the population, based on two types of hotspots for each city in double logarithm scales, we were intrigued by two findings. Firstly, there were no scaling relationships between the area/GDP/CO 2 emissions and population based on the street hotspots, indicated by the very low R 2 values (below 0.01), while significant scaling relationships existed when using NTL hotspots (R 2 values above 0.4). Secondly, the relationships of area-and CO 2 emissions-population were sub-linear (0.84 and 0.68; Figure  7a,c), whereas the GDP-population relationship was super-linear (1.13; Figure 7b), wherein the corresponding scaling exponent values, computed among the chosen 20 cities, were very consistent with values from the recent study based on 287 Chinese prefecture-level cities [17].

Discussion
Cities have long been treated as complex systems. The formation of cities can be described as a dynamic, self-organized, and nonlinear process of human settlements [5], demonstrating highly-heterogenous patterns in both its spatial and aspatial aspects [42]. The spatial aspect can refer to the fractal urban form and the aspatial aspect can refer to the long-tailed distribution of city-related metrics. However, such heterogeneities cannot be revealed effectively since conventional urban data, formed normally through top-down approaches, lack sufficient geographic scope and granularity. In the current geospatial big data era, we can easily conquer this constraint by acquiring fine-grained open data regarding the city form and function at countrywide coverage. Big data is not only big, but also possesses significant fractal and nonlinear properties [43], based on which we can model and analyze a city in a bottom-up manner. That is, delimiting city boundary at the country level or delineating hotspot area at the city scale by agglomeration of individual-based locations.
By adopting the fractal and nonlinear ways of thinking and doing, the cutoff for hotspot boundary derivation was located effectively. Specifically, drawing the border of hotspots is similar to measuring the length of a coastline-a commonality between the two is that, in reality, there is no ground truth for them. The father of fractal geometry, Benoit Mandelbrot [44], has made it clear that the length of a coastline is immeasurable, while the nonlinearity or scaling property is always measurable. In the present study, we characterized the data's nonlinearity in its inherent scaling hierarchy (by head/tail breaks) and power-law or Zipf's law distribution (by the MLE method), by which we obtained the cutoff guiding the spatial clustering. Taking the NTL image as an example, the nested mean values enable us to quickly classify pixels iteratively into a minority of light ones and a majority of dark ones, without exhausting all pixel values by increasing the threshold one at a time. Accordingly, only a few times of experiments on grouping-light-pixel operations for each city led us to generate hotspot polygons whose sizes follow Zipf's law.
The successfully detected Zipf's law of street-and NTL-based hotspots across 20 cities further strengthen the fractal structure of geographic space. It is well-known that a part of a fractal is similar geometrically or statistically to the whole, termed as self-similarity. Since there has been a good agreement among scholars that Zipf's law holds for cities at the country scale [36,45], such a repeated statistical regularity for hotspots at the city scale in the present study can be considered evidence of the self-similarity of geographic space. The self-similarity across multiple scales makes us connect the system of geographic space with that of biology, where similar power law statistics appear across multiple layers in a human body from organs, to tissues, and further to cells [46,47]. Therefore, we believe that Zipf's law can hold within even smaller sub-units than city hotspots (such as neighborhoods), and thus more refined urban center areas could be further identified with the proposed methods. This certainly warrants further study as long as the data granularity allows.
The detected hotspots in both types constituted only a small part of the city area, but accounted for a considerable portion of the urban population, wealth, and energy. This imbalanced ratio between hotspot sizes and the associated socio-economic statistics sheds light on the fact that not all city areas for people live or perform activities. This is also known as the potential problem of the administrative city boundary for urban analysis [21]. Without an accurate capture of human urban activities, the urban scaling estimations may be subjected to unexpected variations. We also examined the power relationship between selected urban measures within the entire administrative boundary among 20 cities, and failed to achieve expected scaling exponents (small R 2 values or in wrong regimes), similar to the case when using the street-based hotspots. By contrast, through the NTL-based hotspots, the derived scaling relationships of area/GDP/CO 2 to population were consistent with the established regimes (e.g., [17,48]). The obtained scaling exponents, shown in Figure 7, indicated that due to a more concentrated settlement and use of infrastructure, the growth of urban economy paced quicker than that of the population (super-linear regime), while the demands of urban areas and the related energy consumption accelerates slower than the population growth (sub-linear regime). The presence of scaling law further implied that the NTL-based hotspots could work as a new, effective instrument for exploring the system of cities.
The hotspots identified by both street and NTL data, by and large, tally with the locations of central urban areas of these 20 cities in China. As noted, street-based hotspots can represent a city's morphological aspects, whereas NTL-based hotspots can accurately reflect a city's functional aspects. The comparison between the two can give us a comprehensive image of how people utilized the urban space. It is noteworthy that the disparity occurs in their spatial distributions. Given that NTL-based hotspots illustrate the aggregation of human activities, we refer that the NTL-based hotspots better manifest the actual urban populous areas than the street-based hotspots, in the context that the street network constructed or traffic planning normally show a time lag. This discrepancy normally hints the evolution of urban centers. That is, these regions are preferred by humans, but apt to be neglected by the municipal authorities or urban scholars. Thus, the planning authorities should at least pay attention to these regions and other urban infrastructure should be strengthened in order to keep pace with real human needs, as well.
By computing IoU metrics, we are able to find that two types of hotspots have less overlays in coastal cities than in inland cities, while coastal cities in China normally have better economic status. Meanwhile, it is worth mentioning that the NTL-based hotspots are very dispersed in the four headmost metropolises, indicating that well-developed cities tend to exhibit a balanced distribution of human activities. It is further referred that cities with higher economic status shift to a more decentralized structure upon urban autonomous development. On this basis, the governments need to take more measures to promote urban justice (including the even distribution of urban resources, etc.) on the process of urban development.

Conclusions
The ultimate goal of city science is closely related to urban smart growth and sustainable development. In natural and societal phenomena, it has been widely adopted that the scaling pattern and power-law statistics are signs of sustainability [49]. This paper provides an intra-urban perspective to study the underlying scaling structure of urban space through novel spatial units: Urban hotspots, detected from geospatial big data including OSM street data and VIIRS imagery. In contrast to conventional spatial units that were imposed by local authorities, the present study adopted the objectively delineated concentration areas as hotspots using the spatial clustering approach. This is mainly motivated by the instability of urban scaling exponents affected by different cities and its sub-unit demarcations. In sum, we found (1) that Zipf's law also holds strikingly at the intra-urban level; and (2) that NTL-based hotspots can be good proxies for city populous areas, by which the urban scaling relationship can be correctly maintained.
The method for hotspot detection acts as a promising tool and could supplement innovative urban planning toolboxes in the big data era. Despite the strengths of urban hotspot in this work, there is still room for improvement in terms of the following. Firstly, whether the intra-urban scaling law exists in other countries remains to be verified from a global view, in addition to these 20 cities in China. Secondly, it is important to add NTL images before 2020 to check whether and how the intra-urban scaling exponents change or evolve. Further, the updated raster data sets of GDP, population, and CO 2 emissions after 2010 will be combined once they are available, for eliminating possible biases or inaccuracies that occurred due to the difference in data time acquisition. Thirdly, the multiscale effect of scaling analytics (e.g., detecting a more refined spatial unit and related power law statistics) within one city needs to be further conducted. Fourthly, the underlying mechanism of this scaling law has not been revealed yet, concerning policy, landform or demographic traits, etc. Future work will point to these directions. Table A1. The candidate cutoffs for spatial clustering of light pixels for hotspot detection and their corresponding power law metrics. (Note: NA: Not available; for most of cities, the number of scaling hierarchies of NTL image pixel lightness ≥5, meaning that there were at least 5 iteratively-averaged pixel values leading to an imbalanced ratio between dark-to-light pixels per city image. Among those five candidate cutoffs, the third one appeared to be the most suitable for urban hotspot delineation, as the almost related α values of hotspot areas were closest to 2 and with an acceptable p).

Appendix B
This appendix supplements Section 3.1 by presenting the overlapping ratios between street-and NTL-based hotspots among 20 cities. We adopted IoU for assessing how much one type of hotspot overlaps another in a city. The IoU metric between two types of hotspots can be denoted by the following equation: where Area s is the total area of street-based hotspots, Area n is the total area of NTL-based hotspots. The results of IoU for each city is shown in Table A2.