Is Spatial Resolution Critical in Urbanization Velocity Analysis ? Investigations in the Pearl River Delta

Grid-based urbanization velocity analysis of remote sensing imagery is used to measure urban growth rates. However, it remains unclear how critical the spatial resolution of the imagery is to such grid-based approaches. This research therefore investigated how urbanization velocity estimates respond to different spatial resolutions, as determined by the grid sizes used. Landsat satellite images of the Pearl River Delta (PRD) in China from the years 2000, 2005, 2010 and 2015 were hierarchically aggregated using different grid sizes. Statistical analyses of urbanization velocity derived using different spatial resolutions (or grid sizes) were used to investigate the relationships between socio-economic indicators and the velocity of urbanization for 27 large cities in PRD. The results revealed that those cities with above-average urbanization velocities remain unaffected by the spatial resolution (or grid-size), and the relationships between urbanization velocities and socio-economic indicators are independent of spatial resolution (or grid sizes) used. Moreover, empirical variogram models, the local variance model, and the geographical variance model all indicated that coarse resolution version (480 m) of Landsat images based on aggregated pixel yielded more appropriate results than the original fine resolution version (30 m), when identifying the characteristics of spatial autocorrelation and spatial structure variability of urbanization patterns and processes. The results conclude that the most appropriate spatial resolution for investigations into urbanization velocities is not always the highest resolution. The resulting patterns of urbanization velocities at different spatial resolutions can be used as a basis for studying the spatial heterogeneity of other datasets with variable spatial resolutions, especially for evaluating the capability of a multi-resolution dataset in reflecting spatial structure and spatial autocorrelation features in an urban environment.


Introduction
The rate of urbanization all over the world is quite alarming, with the proportion of the world's population living in urban environments projected to reach 66% by 2050 [1].However, the understanding of urbanization is primarily based on population figures obtained from the United Nations; these statistics do not include any information on the population spatial distributions, or evolution metrics of built-up areas within urban environments [2][3][4].
Urbanization, urban growth, urban sprawl and urban expansion are different concepts that have caused much confusion in multidimensional urban systems analysis.For a differentiation, it can be noted that urbanization can be viewed as a characteristic of the population as a particular kind of land use and land cover, as well as a characteristic of social and economic processes and interactions affecting both population and land [5,6].Urban growth mainly refers to an increase in urban population size, independent of rural population [7].Urban sprawl is treated as a process that focuses on describing pattern of land-use in an urbanized area through eight distinct dimensions: density, continuity, concentration, compactness, centrality, nuclearity, diversity, and proximity [5,8].There is no specific definition of urban expansion-this concept is commonly used to describe urban population, physical expansion, quality of urban layout, land and housing regulation and so on [9,10].In most of real world situations, these terms cannot be clearly separated, since urbanization, urban growth, urban sprawl and urban expansion are highly interlinked.However, it is important to realize that the huge growth of urban population may cause uncontrolled urban growth, resulting in urban sprawl and urban expansion.Urbanization may also result from and contribute to urban growth, urban sprawl and urban expansion [5,9].
Although the complexity of these four terms and their ambiguous and partially overlapping meanings make it difficult reach a consensus of a distinct urban phenomenon analysis, a variety of urbanization parameters from the standpoint of the built environment have been proposed to describe urbanization trends [11][12][13][14].For instance, some of these used specific landscape spatial metrics to characterize the configuration and composition of urban areas (e.g., [15][16][17][18][19]).While many such parameters have been identified, it is generally difficult to distinguish between those that are useful and those that are not.Some studies have measured the accessibility within each city on the basis of gravity transportation modelling (e.g., [12,20,21]), but the acquisition of these transport network datasets is problematic.Generally, there appears to be no consensus between those investigating urban landscapes on the parameters to use for urbanization velocity evaluation.This research therefore provides below a succinct review of the most widely accepted and commonly used parameters for measuring urbanization in order to provide a descriptive framework that can be used for measuring urbanization velocity [8].
A number of studies have also shown that remote sensing data and associated techniques are advantageous for classifying, monitoring and analyzing urbanization and its development over time at a range of scales, with an emphasis on mapping large areas at a time [22][23][24].Urbanization velocity (also called urban expansion speed) is defined as the annual growth rate of urban area within a period.It indicates the absolute differences (in terms of footprints) of urban areas within a certain time period [25].The measurement of urbanization using remote sensing imagery has been widely used for mapping, quantification, and documentation of the extent, growth rates, and percentage change in urban areas [26][27][28], and can also be used to predict possible future urban growth [29,30].Especially the grid-based urbanization velocity analysis, which typically involves the use of a "grid-based moving window" on remote sensing imagery to detect the spatial gradient changes of grid-based land use and land cover [31][32][33], nighttime luminosity changes [30,34], population changes [35], or temperature changes [36] through time, were popularly defined to describe the urbanization process.Grid-based urbanization velocity analyses do not require extensive auxiliary data to obtain spatial-temporal urban growth [27,37], it is able to avoid some of the redundancies of estimations that are caused by many landscape parameters [8,38,39] when attempting to describe quantitatively the human settlement patterns [37,40,41].It aims to quantify the local urbanization velocity across a landscape [32,42,43].Spatial and temporal grid-based gradient information can therefore provide consistent, spatially explicit parameters that can be used to record the expansion of urban areas by estimating the speed and direction of urban growth [28,44].
Grid-based urbanization velocity analyses are, however, sensitive to the size of the "grid" (or the spatial resolution) used to compute spatial and temporal characteristic of the urbanization process [32,45,46].Grid-based urbanization velocity analyses can, in theory, be carried out at any scale and descriptions of urban growth may therefore need to be tied to a specific scale.There is little knowledge about the effects of spatial resolution on urban velocity analyses [47].At very high spatial resolutions, transitional patches of human settlement tend to be complex, while coarser resolutions tend to smooth out the effects of urbanization.By contrast, although a coarse resolution urbanization velocity map may show more instances of the same urban expansion features than a high resolution map, it may demonstrate the agglomeration of land cover types and reveal greater urbanization evolution process that are not easily detected from the high resolution map [48].The challenge is therefore to establish a meaningful and useful multi-resolution urbanization velocity model to investigate the spatial resolution issue on grid-based urbanization velocity analyses.
The first objective of this research was to investigate the effect of imagery spatial resolution on grid-based urbanization velocity analysis, which required an investigation into the most appropriate size to use for the "moving window" (expressed as the grid size, which is an effective surrogate for spatial resolution).Since the widespread use of remote sensing data has generally increased interest in studying the capability issues relating to spatial resolution of images, the second objective was thus to investigate whether or not the available fine resolution remote sensing data for grid-based urbanization velocity analyses in the case study of Pearl River Delta are appropriate for the spatial resolution at which the investigated processes operate, or the spatial resolution at which decisions are required.The third objective was to take the grid-based urbanization velocity as the independent variable, and test the inherent characteristics of spatial autocorrelations and spatial structure heterogeneities.These characteristics are useful to effectively validate the results of grid-based urbanization velocity analyses at different spatial resolution.
Increasing numbers of flexible variance methods are providing effective ways of identifying spatial resolution thresholds and spatial resolution domain problems.Frequently used methods for dealing with different spatial resolutions (referred to as multi-scale methods) include local variance analysis [49,50], geographical variance analysis [51,52], semivariance analysis [53,54], multifractal analysis [55,56], wavelet transform analysis [57,58], and Fourier transform analysis [59,60].All of these methods are able to quantify landscape characteristics of different spatial resolutions in their mathematical formulations and procedures, but their selection heavily depends on the nature of the data and the objectives of the investigation [53].
In this study, the grid-based urbanization velocity analysis is based on spatio-temporal changes within an urban area compared to neighboring areas, using a moving window.This research systematically organizes the moving window (i.e., spatial resolution or grid size of remote sensing data) into a hierarchical, grid-based nested dataset.Through a thorough and succinct literature analysis [50,51,53,[61][62][63][64][65], the empirical variogram model, local variance model and the geographical variance model, were selected to investigate the urbanization velocity analysis results at different spatial resolutions.They show the advantages of analyzing the spatial autocorrelation and spatial structural heterogeneities for hierarchical grid size.Meanwhile, these three methods have rarely been used to identify spatial resolution thresholds and domains within urbanization velocity analyses.This research also can quantify and validate their capability in investigating the influence that the hierarchical grid size (i.e., spatial resolution) has on quantitative descriptions of urban growth.The investigation of grid size (i.e., spatial resolution) effects on urbanization velocity analysis was set out in the Pearl River Delta study area between the years of 2000 and 2015.

Study Area
The Pearl River Delta (PRD) in southeastern China extends from 21 • N to 24 • N and from 112 • E to 115 • E (as shown in Figure 1).The twenty seven cities of the PRD are home to a population of more than 55 million.According to the World Bank Group (2015), the PRD has become the largest "megaregion" in the world in terms of both surface area and population.Long-term monitoring by Taubenböck et al. [23] has revealed that its spatial extent in 2011 was 13.14 times greater than in 1975, making this megaregion among the most dynamic areas in the world, even outperforming the two-dimensional spatial growth rates of China's other megaregions such as Shanghai (about 6 times) or Beijing (7.5 times) [66].A comprehensive evaluation of multi-resolution urbanization velocity of this emerging megacity in the PRD, which is the objective of this study, is expected to make a valuable contribution to ensuring the sustainable development of the area.
Remote Sens. 2017, 9, 80 4 of 25 or Beijing (7.5 times) [66].A comprehensive evaluation of multi-resolution urbanization velocity of this emerging megacity in the PRD, which is the objective of this study, is expected to make a valuable contribution to ensuring the sustainable development of the area.

Data Collection
The imageries used for the land use and land cover (LULC) classification were at a resolution of 30 meters (hereafter referred to as 30 m) and obtained from Landsat Surface Reflection products (Landsast 4/5-TM; 2000, 2005, 2010) and Landsat 8 archive products (Landsat 8 -OLI/TIRS; 2015).As shown in Table 1, eight Landsat images can cover the whole areas of PRD in each year.The time period for each year was defined to be from July and October, which revealed the most abundant leafage and minimum cloud cover (<0.1%).Since PRD has a subtropical climate and the spectral information of the typical vegetation areas do not show dramatic changes during the period July to October, there are moderate seasonal variations of the NDVI values in this area for the time period

Data Collection
The imageries used for the land use and land cover (LULC) classification were at a resolution of 30 meters (hereafter referred to as 30 m) and obtained from Landsat Surface Reflection products (Landsast 4/5-TM; 2000, 2005, 2010) and Landsat 8 archive products (Landsat 8 -OLI/TIRS; 2015).As shown in Table 1, eight Landsat images can cover the whole areas of PRD in each year.The time period for each year was defined to be from July and October, which revealed the most abundant leafage and minimum cloud cover (<0.1%).Since PRD has a subtropical climate and the spectral information of the typical vegetation areas do not show dramatic changes during the period July to October, there are moderate seasonal variations of the NDVI values in this area for the time period chosen.The one-kilometer LULC data for 2000, 2005, 2010 collected from the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences (RESDC) served as auxiliary data for the land use classification through generating coarse resolution masks of LULC classification in PRD.
High resolution imagery data (<5 m), which were available online from Google Earth, and urban footprint data with a spatial resolution of 30 m, obtained from the German Aerospace Center ( [66]), were used as reference datasets to provide appropriate ground truth information to validate LULC classification over the PRD area.
The original census data were obtained from the Guangdong Statistical Yearbooks for 2000, 2005, 2010, and 2014 (http://www.gdstats.gov.cn/tjsj/gdtjnj/).Census data from these four years (which are available for 27 original city-scale administrative units) were used to obtain the socio-economic indicators, including GDP, GDP per capita, Total value of retail sales, Total number of Job, Income per capita, Deposit balance, Population.This research employed these socio-economic indicators as the criteria used to investigate the effects of different spatial resolutions in urbanization velocity analyses.

Workflow
A hierarchical grid-based urbanization velocity framework (Figure 2) was established based on Landsat data with variable sized moving windows to investigate how the size of the grid size influences urbanization velocity analysis results.The variable resolution grid-based urbanization velocity model is described below, followed by a description of the methods used to analyze the spatial correlation, spatial structure and spatial autocorrelation characteristics of the urbanization velocity analysis results at each grid size (e.g., spatial resolution).Finally, all these models were combined to determine the optimum spatial resolution to use in grid-based urbanization velocity analysis.

Workflow
A hierarchical grid-based urbanization velocity framework (Figure 2) was established based on Landsat data with variable sized moving windows to investigate how the size of the grid size influences urbanization velocity analysis results.The variable resolution grid-based urbanization velocity model is described below, followed by a description of the methods used to analyze the spatial correlation, spatial structure and spatial autocorrelation characteristics of the urbanization velocity analysis results at each grid size (e.g., spatial resolution).Finally, all these models were combined to determine the optimum spatial resolution to use in grid-based urbanization velocity analysis.

Land Use and Land Cover Classification
The atmospherically and geometrically corrected Landsat-4/5 TM surface reflection products were used to facilitate the land surface changes directly.The Landsat-8 archive images were converted to the surface reflection products through the equations provided by the US Geological Survey [67].
Despite the fact that some influences from satellite drift/changeover, incomplete calibration, correction loss, images monthly variation, as well as atmospheric effects may still exist, these influences are typically much smaller compared with those caused by environmental drivers.Hence, in this study, such potentially existing differences between images were ignored.
Object-based image analysis (OBIA) was employed to identify the LULC classification based on Landsat 4/5 TM and Landsat-8.When compared with the traditional per-pixel methods in image classification, OBIA approaches have already demonstrated some advantages in pattern recognition, such as exhibiting useful shape, context and texture information for urban phenomenon analyses [68,69].This research aims to explore the pattern characteristics of built-up areas under the urbanization process.Here, OBIA is an alternative to dividing the built-up areas in Landsat images into multi-scale built-up objects.The pixels in each semantic object share relatively similar spatial, textual and contextual information to avoid "salt and pepper" phenomena that commonly appear in per-pixel analyses.
A multi-level segmentation/and classification OBIA framework was designed within the eCognition Developer 9 software in order to create spatially contiguous and homogeneous objects.The most appropriate segmentation scale (scale = 15) was determined by use of the scale estimation tool ESP2 [70].Then the fuzzy classification approach which was developed by Wei et al. [69] for Landsat images was used to classify the data into two classes: built-up areas and non-built-up areas.Fuzzy ruleset-based classifications have the advantages of addressing the uncertainty in the Landsat surface reflection products values.It uses n-dimensional tuples of memberships to describe the degree of class assignments, making the ruleset transferable to all the four years of images [69,71,72].Hence, fuzzy sets were employed to handle most sources of vagueness in the spatial-temporal LULC classification of Landsat surface reflection products.
As a first step, the LULC was classified into five classes: urban, rural, vegetation areas (excluding forest), forest and water following the approach of Wei et al. [69].However, this research mainly focuses on measuring the annual growth rate of built-up areas within a period.Therefore, in a second step the LULC classification was simplified through the following three-step process: (1) The 1-km LULC cover data for 2000, 2005, 2010 provided by RESDC were used to generate coarse resolution masks of built-up areas, which served as auxiliary data for the following steps.(2) The built-up areas were determined by using the normalized difference vegetation index (NDVI, the fuzzy membership ranged from 0.24 to 0.38), the normalized difference water index (NDWI, the fuzzy membership ranged from 0.04 to 0.08), the specific leaf area vegetation index (SLAVI, the fuzzy membership ranged from 0.53 to 0.58), the land and water masks index (LWMI, the fuzzy membership ranged from 105 to 200), and the normalized difference moisture index (NDMI, the fuzzy membership ranged from 0.03 to 0.05).The membership functions for each feature were combined by the fuzzy-logic operator.And (3) all other areas were assigned as non-built-up areas.Before carrying out the fuzzy classification, the eight images were mosaicked for each year.

Urbanization Velocity Estimation at Different Spatial Resolutions
As described above, the velocity of urbanization is defined as annual growth rate of built-up areas growth within a period.It indicates the absolute differences in built-up areas footprints within a certain time period.The annual growth rate of built-up areas is based on the unit of a grid, i.e., a raster cell.Therefore, when using multiple scales, such grid-based urbanization velocity estimation compares the built-up densities changed in each unit with its neighborhood areas within a moving window over a specific time periods.It is closely akin to the index of instantaneous horizontal velocity of temperature change and nighttime-light change proposed by Loarie et al. [46] and Ma et al. [30].
As shown in Equations ( 1)-( 4), the grid-based urbanization velocity estimations were based on the areas classified as "built-up" in the 30 m spatial resolution Landsat images from 2000, 2005, 2010 and 2015.
V urban = T trend G spatial (1) (2) where T trend is the temporal trend of images from 2000 to 2015, and t ave is the mean value of the year.G spatial is the mean annual density of build-up areas at the location of each grid from 2000 to 2015.L i stands for the pixel value of built-up areas density within the image.G spatial_i stands for the spatial gradient of built-up areas density at the location of each grid, which is estimated by the average maximum technique [73] based on the 3 × 3 grid cell neighborhood (as shown in Figure 3 and Equations ( 3) and ( 4).For the urbanization velocity (UV) estimates processing at the boundary of PRD, the boundary pixels were assigned as the null values, and the author only took into consideration those moving windows that did not contain any null values.
As shown in Equations ( 1)-( 4), the grid-based urbanization velocity estimations were based on the areas classified as "built-up" in the 30 m spatial resolution Landsat images from 2000, 2005, 2010 and 2015.
_ 8) where Ttrend is the temporal trend of images from 2000 to 2015, and tave is the mean value of the year.Gspatial is the mean annual density of build-up areas at the location of each grid from 2000 to 2015.Li stands for the pixel value of built-up areas density within the image.Gspatial_i stands for the spatial gradient of built-up areas density at the location of each grid, which is estimated by the average maximum technique [73] based on the 3 × 3 grid cell neighborhood (as shown in Figure 3 and Equations ( 3) and ( 4).For the urbanization velocity (UV) estimates processing at the boundary of PRD, the boundary pixels were assigned as the null values, and the author only took into consideration those moving windows that did not contain any null values.probably unrealistic for the individual sensor case, they may be statistically valid for coarse resolution of remote sensors [62].The objective was to investigate the basic effect of different spatial resolutions on grid-based urbanization velocity estimation, rather than to mimic the response characteristics of a particular remote sensing sensor.This approach avoids the controversies over the exact definition of spatial resolution and how it should be measured [62,74].

Statistical Test for the Dependence of Urbanization Velocity on Spatial Resolution
How to effectively validate the scale dependent/independent characteristics of urbanization velocity (UV) estimations across multiple grid sizes (i.e., different spatial resolutions) is the biggest challenge.This research therefore approached the UV analysis for different spatial resolutions from a statistical perspective, while focusing on concepts of spatial dependence on spatial resolution and correlation.The scale-dependence of the UV analysis results was tested by investigating the Pearson correlation coefficients between the UV analysis results and the corresponding city-based socio-economic indicators.All 27 of the large cities in the Pearl River Delta are individual municipalities as well as self-contained economic areas.It is therefore appropriate to calculate the geometric mean of the UV estimates on a city-scale.The annual growth rates of socio-economic indicators for the 27 cities of the Pearl River Delta are fixed values and not influenced by the spatial resolution.They can therefore also be used to establish the relationships with multi-level urbanization velocity results to explore the scale influence.Seven typical socio-economic indicators (including GDP, GDP per capita, Total value of retail sales, Total number of jobs, Income per capita, Deposit balance, Population) from the Chinese statistical yearbooks were selected to do a correlation analysis with UV estimates.These socio-economic indicators represent specific urbanization processes at specific instants in time and can therefore be used to investigate the spatial resolution issues in urban growth analysis.

Quantitative Descriptions of Spatial Resolution Domains and Thresholds
Characterizing a landscape's spatial variability helps to identify a suitable resolution for capturing the spatial grid size of the surface properties and for optimizing the field data collection [53,75].Through a thorough and succinct literature review [50,51,53,54,61,62,[76][77][78], this research divided the spatial variability analyses of multi-resolutions UV estimates into two groups: Spatial autocorrelation analysis, which is a linear combination of empirical variogram models used to calculate the structural parameters of the semivariogram from the spatial autocorrelation perspective to evaluate the differences in spatial variability at different spatial resolutions and to characterize the impact of the spatial resolution.
Spatial structure analysis, which involves using the geographical variance and local variance methods to quantify the level of variability in spatial structures of different grid sizes within a landscape.The mean extent of the grid size with the highest level of variability is defined as the "proper spatial resolution" to capture the spatial variability of the landscape pattern.

Empirical Variogram Method
Empirical variogram models are powerful geostatistical tools for understanding the loss of spatial autocorrelation and spatial variability as spatial resolution decreases [49,79].The nature of a variogram is also an important property characterizing the spatial continuity and regularity of the data at different resolutions, based on a probabilistic approach considering the image as one of many possible realizations of a second-order stationary random function expressed as follows: where γ(h) is the modeled variogram, δ 2 is the variogram sill (i.e., overall image variance), b k is the proportion of overall image variance, δ 2 related to each range γ k , γ k is the variogram range related to each elementary variogram function g k .This use of theoretical variogram models has been fully described by [53,54].
To detect the multi-scale spatial autocorrelation of the UV analysis results, the exponential, spherical, Gaussian, Matern, and nugget effect models were tested to identify the best fits of the experimental variogram models for each scale.The sill δ 2 and the variogram range γ k were used as spatial variability parameters for quantifying the overall spatial variance of multi-scale datasets [54].The attributes of land use classification results with highly skewed data distribution may present problems in variogram calculation.Hence, the normal test was performed before establishing the variogram models.The variogram interpretations were based on certain lag distance, the normal test were thus based on the stratified random samples.The Shapiro-Wilk normality test showed that p-values of all the datasets were larger than 0.05.
Empirical variograms for different orientations have been reported by Garrigues et al. [54] to showed little variation, so this research only compared the omnidirectional variograms that represented simultaneously the horizontal, vertical, and two diagonal directions of the multi-resolution pixel grids in this study.The nugget effect in the variogram model relates to either uncorrelated spatial structures or randomly spatial structure features that are evident at spatial scales smaller than the pixel size.

Local Variance Method
Woodcock and Strahler [62] proposed a local variance method to find an appropriate resolution for understanding the spatial structure of multi-resolution images.Local variance is calculated as the mean value of the standard deviation of a moving 3 × 3 pixel window over an entire image.Graphs of local variance in images with different spatial resolutions can be treated as a function of the resolution in order to describe the relationship between the sizes of the objects in the images and the spatial resolution.
The results of the urbanization velocity analyses based on moving windows of different sizes can be considered to represent the original image data that have been degraded to several levels of lower spatial resolutions.The graphs of local variance can be used to determine the scales at which the objects in the built-up areas most closely approximate those in the surrounding areas within a complex urban system.If, for instance, the spatial resolution is considerably finer than the size of the objects in urban growth areas within the image, most of the objects in the images will correlate well with neighbors' objects and the local variance will be low.If, however, the sizes of objects in urban growth areas approximate the pixel size (i.e., the spatial resolution), then the likelihood that they will resemble neighbors' objects is reduced and the local variance will increase.As the pixel size (i.e., the spatial resolution) increases further until one pixel may contain a number of objects, the local variance will start to decrease.The peak of the graph therefore indicates that the point at which the size of the moving window (or the spatial resolution) matches the size of the objects in the areas of urban growth [51,62].The objective was to investigate and determine the grid size (i.e., spatial resolution) for moving windows that can capture a major proportion of the urban growth.

Geographic Variance Method
The geographic variance (GV) method proposed by Moellering and Tobler [52] has the potential to consistently and logically detect and describe landscape objects of different sizes.The GV procedure allows the relative spatial variability to be measured, as well as the independent contribution of spatial variability made by each grid size (i.e., the spatial resolution) to a nested hierarchy.The multi-resolution urbanization velocity analysis results can be treated as regularly grouped hierarchical data to which GV analysis can be applied and in which the total variance of UV estimates within the PRD can be partitioned hierarchically at different grid sizes.Moellering and Tobler [52] derived the scale variance components for a 3-level (∂, β, γ) hierarchy in which the total variation of the system is expressed as the total sum of squares.As shown in Equation 6, γ represent the original spatial resolution of the image, ∂ represent the grid group aggregated from γ level, β represent the grid groups aggregated from γ level.

SS Total
where I is the number of ∂ level groups, J i is the number of β level groups in each ith ∂ group, and K ij is the number of γ level groups in each ijth β level group.X ijk stands for the UV estimates at each pixel at the γ level, and N is the total number of X ijk observations.n i is the number of observations in the ith region, and n ij is the number of observations in the jth state of the ith region.This method has been fully described by [52].
When plotting the scale variance of the hierarchical levels of grid-size dataset, a peak would imply that high variability occurs at the corresponding grid size (i.e., the spatial resolution).A peak would also reflect the relative contribution of that particular grid size (i.e., the spatial resolution) to the total variability of the hierarchical dataset [78,80].

Land Use/Land Cover Accuracy
Two methods were employed to assess the accuracy of the LULC classification.In the first accuracy assessment method, the sample design tool in ARCgis 10.3.1 was used to randomly generate 2000 reference samples for each image.Then the reference data were assigned the attributes based on the high resolution SPOT images from the online Google Earth portal.These high resolution historical remote sensing maps can provide appropriate ground truth information to validate our LULC classification results.Taking the year of 2010 as an example (as shown in Figure 4), 18.5% of the samples were attributed as built-up areas, and the rest samples were assigned as non-built-up areas.The same procedure was employed to the other three years.Then all the samples were used to carry out the combined analysis of the 30 m resolution LULC classification results.Finally, the overall accuracy of each year was obtained through their respective error matrixes (as shown in Table 2).The overall accuracies of the four images yielded values of 92.35%, 89.60%, 90.45%, and 87.95%.These classification accuracies are clearly above a generally assumed satisfactory classification overall accuracy of 85%.
In the second accuracy assessment method, two years (the years 2000 and 2005) of urban footprint data in the Pearl River Delta provided by Taubenböck et al. [23,24] were used to validate the built-up areas classification.These urban footprint data were originally derived from the combination of high resolution TerraSAR-X strip maps images (3 m resolution) in combination with Landsat images.Taubenböck et al. [23,24] mainly selected random distribution points separately in different cities around the Pearl River Delta (e.g., Shenzhen, Guangzhou, Dongguan, Foshan, Zhaoqing, Zhongshan, Huizhou, Jiangmen) by visual comparison to available VHR Ikonos and QuickBird data as well as topographic reference information from these cities.Their overall accuracies range from 80% to 90%.Thus, these urban footprint maps prove high and stable accuracy values on the regional scale, and they are powerful to demonstrate the spatial connectivity between cities based on the continuity of built-up areas patterns [53,78].The results show that the proportions of overlay areas between urban footprint and built-up areas in this research were 90.3% in 2000 and 86.5% in 2010.This high consistency of the built-up areas and the urban footprint information strongly supports the urbanization velocity detection based on the patterns of built-up areas.

Multi-Scale Urbanization Velocity
Overall mean geometric velocities of urbanization in the whole PRD area ranged from 2.10 sq•km•year −1 to 3.67 sq•km•year −1 across the multiple moving window sizes.As shown in Figure 5, twenty-seven major cities in the PRD have undergone varying degrees of urbanization velocity.Based on the geometrical interval theory in Arcgis 10.3.1 [81], the urbanization velocity can be categorized into four levels: very low (≤0), low (0-5), median (5-30), high (≥30).More specifically, nine cities (Huadu, Guangzhou, Nanhai, Dongguan, Panyu, Foshan, Shunde, Zhongshan, Shenzhen) revealed above-average urbanization velocities across multiple grid sizes (Figure 6).As the grid sizes increased to Level 16, the urbanization velocity for the city of Zengcheng gradually increased and eventually reached the average urbanization velocity for all 27 cities.All of these cities gradually expanded, becoming interconnected and coalescing to form increasingly continuous built-up areas between 2000 and 2015 through con-urbanization.The urban growth velocity analyses (shown in Figures 5 and 6) revealed that urbanization velocities varied for each city, but that the grid size had no statistically significantly influence on the magnitudes of urbanization velocities determined for these 27 cities.

Multi-Scale Urbanization Velocity
Overall mean geometric velocities of urbanization in the whole PRD area ranged from 2.10 sq•km•year −1 to 3.67 sq•km•year −1 across the multiple moving window sizes.As shown in Figure 5, twenty-seven major cities in the PRD have undergone varying degrees of urbanization velocity.Based on the geometrical interval theory in Arcgis 10.3.1 [81], the urbanization velocity can be categorized into four levels: very low (≤0), low (0-5), median , high (≥30).More specifically, nine cities (Huadu, Guangzhou, Nanhai, Dongguan, Panyu, Foshan, Shunde, Zhongshan, Shenzhen) revealed above-average urbanization velocities across multiple grid sizes (Figure 6).As the grid sizes increased to Level 16, the urbanization velocity for the city of Zengcheng gradually increased and eventually reached the average urbanization velocity for all 27 cities.All of these cities gradually expanded, becoming interconnected and coalescing to form increasingly continuous built-up areas between 2000 and 2015 through con-urbanization.The urban growth velocity analyses (shown in Figures 5 and 6) revealed that urbanization velocities varied for each city, but that the grid size had no statistically significantly influence on the magnitudes of urbanization velocities determined for these 27 cities.

Lack of Scale-Dependence in Urbanization Velocity (UV) Estimations
The Pearson correlation analyses indicate that there are only three socio-economic indicators, (including GDP, total value of retail sales, total number of jobs) which revealed a moderate relationship with urbanization velocity results based on the 30 m resolution Landsat images (Level 1).Therefore, these three indicators were employed to carry out the correlation analyses with other spatial resolutions for the UV estimates.Table 3 shows in great detail that 8 out of 10 relationships between UV estimates and socio-economic indicators using different spatial resolutions yielded moderate correlations (>0.396), with the exceptions being for two UV analysis results at Level_256 resolution (including totally 727 pixels) and Level_512 resolution (including totally 182 pixels).Although the missing pixels caused by the coarse resolution at these two low spatial resolution levels influenced the correlation coefficients between UV estimates and socio-economic data, the mutual relationships between UV estimates and socio-economic indicators are generally stronger than the influence of spatial resolution.It therefore appears to be reasonable to conclude that the inter-relationships between the city-based UV analysis results and socio-economic indicators are independent of the spatial resolution.

Lack of Scale-Dependence in Urbanization Velocity (UV) Estimations
The Pearson correlation analyses indicate that there are only three socio-economic indicators, (including GDP, total value of retail sales, total number of jobs) which revealed a moderate relationship with urbanization velocity results based on the 30 m resolution Landsat images (Level 1).Therefore, these three indicators were employed to carry out the correlation analyses with other spatial resolutions for the UV estimates.Table 3 shows in great detail that 8 out of 10 relationships between UV estimates and socio-economic indicators using different spatial resolutions yielded moderate correlations (>0.396), with the exceptions being for two UV analysis results at Level_256 resolution (including totally 727 pixels) and Level_512 resolution (including totally 182 pixels).Although the missing pixels caused by the coarse resolution at these two low spatial resolution levels influenced the correlation coefficients between UV estimates and socio-economic data, the mutual relationships between UV estimates and socio-economic indicators are generally stronger than the influence of spatial resolution.It therefore appears to be reasonable to conclude that the inter-relationships between the city-based UV analysis results and socio-economic indicators are independent of the spatial resolution.

Spatial Autocorrelation and Spatial Variability of UV Analysis Results Based on Different Spatial Resolutions
The empirical variogram models are shown in Figure 7 and the respective model parameters can be found in Table 4. Previous studies illustrate that the semivariance value in the empirical omnidirectional variograms models consistently decreases as nominal pixel size increases [54].However, from Figure 7 it can be seen that the semivariance value in the omnidirectional variogram models does not decrease consistently as the nominal pixel size increases, which appears to contradict the general hypothesis and findings to date.In fact the semivariances alone are of limited use for characterizing the impact that the spatial resolution has on the results of land use and land cover analyses [53].This research therefore evaluated the differences in spatial autocorrelation and spatial variability between UV estimates obtained using different spatial resolutions by comparing different variogram model parameters (including the variogram sill ∂ 2 , the range, and the rate of spatial variability), in order to estimate the overall image variance.The UV estimates at 10 spatial resolution levels have been generated using linear regionalization variogram models in Figure 7.The models' sill ∂ 2 generally increases as the window size (i.e., spatial resolution) increases from Level 1 resolution to Level 8 resolution, but then decreases as the window size increases further from Level 8 resolution to Level 512 resolution.The empirical variogram models are flat, with nugget effects at Level 16 resolution and Level 32 resolution.This means that there were no correlations between UV analysis results and the spatial extent of the image at these two levels of spatial resolution, indicating spatially independent variance of UV estimates.The high variability of the sill ∂ 2 at Level 8 resolution is explained by the intrinsic variability in the structure of built-up areas (e.g., variations in the density and connectivity of built-up areas, etc.), as well as vegetation areas that contrast with the built-up areas.In this case, the UV analysis results at Level 8 resolution were better able to represent the intrinsic variability of urbanization in the Pearl River Delta.
In the empirical variogram model, the sill is an indicator of the overall spatial variance of the data.If the sill is not reached before a maximum distance of image, it implies that the spatial extent of the image is not sufficiently large to encompass the low frequency variation in the data.The maximum distance is usually defined as one-third of the full spatial extent of the image subsets [54].In this research, most of the variograms reach a sill within one-third of the spatial extent of the image, except for the Level 4 resolution and the Level 64 resolution.The characteristics of the spatial structures of UV estimates at Level 4 and Level 64 are consequently very uncertain, since the spatial extents of the images are too small to detect the spatial variability of UV analysis results.Moreover, the variograms from Level 128 resolution through to Level 512 resolution increase at a very short spatial extent in the image, indicating that the UV analysis results at these three spatial resolutions are roughly structured and spatial structures of UV analysis results have dimensions similar to the spatial resolution of the image.In order to further demonstrate the rate of spatial variability, the authors assumed that the variogram derived from a 30 m resolution moving window approximates the actual spatial variability of the landscape.The decreasing or increasing sill rates relative to the sill at 30 m resolution then correspond to the loss or gain in the spatial variability of UV estimates as the pixel size increases.As shown in Table 4 (column 5), intra-pixel spatial heterogeneity of UV estimates increases rapidly with increasing pixel size until this size reaches the Level 8 resolution.From this point on spatial heterogeneity of UV estimates gradually decreases with lower resolution.UV analysis results at the Level 8 resolution indicate an overall spatial variability that is 3.077 times greater than at the original resolution, and most of the UV estimates' spatial variability is lost at the Level 512 resolution.The rate of decay of sill served as a reference for evaluating spatial information decay in multi-resolution UV analysis results.
Generally speaking, the parameters of variogram sill ∂ 2 , the range of the variogram model, and the rate of spatial variability all characterize the important variability of UV analysis results at different spatial resolutions.The variogram models highlight that at Level 8 resolution (i.e., at 240 m spatial resolution) the intra-pixel spatial heterogeneity is fully explained by this spatial resolution.It is therefore a suitable resolution to capture the major part of the variability in the spatial structural of UV analyses results in the Pearl River Delta.

Local Variance and Geographical Variance Of the UV Analysis Results at Different Spatial Resolutions
Local variance as a function of the grid size used for the UV estimation is presented graphically in Figure 8.The graph begins with low local variance at the initial 30 m resolution of the image.Due to the dominance of urban agglomeration processes in the PRD over the past 15 years, most of the pixels have similarly high UV estimates.The standard deviation of a 3 × 3 window based on the original resolution is therefore very low.As the grid size increases, the number of pixels with a high UV estimates decreases and the likelihood that surrounding pixels will be similar also decreases.In this situation the local variance increases.This trend continues until a peak in local variance is observed at Level 16 (i.e., at about 480 m coarse resolution).As the resolution increases beyond this peak, pixels tend to become aggregated to the extent that they again become more similar and the local variance decreases.Level 16 therefore indicates a moving window size (or spatial resolution) that matches the sizes of urban growth units at which a large proportion of urban growth processes occur.
For an additional city-level analysis, the whole of the PRD was further partitioned into city-based scenes.Their corresponding graphs of local variance were investigated in order to identify the moving window sizes that best matched urban growth units in the different cities.As shown in Figure 9, for the cities in the central urban agglomeration areas such as Guangzhou, Foshan, and Dongguan, the peaks in local variance were observed at low levels (ranging from Level 2 to Level 4), or at 60 m to 120 m fine resolution.In contrast, cities in the outer areas of the PRD yielded local variance peaks at higher levels (Level 8 to Level 32) or at coarser resolutions (480 m to 960 m).This means that the most appropriate spatial resolution for analyzing urbanization velocity may be different in the core areas (where urbanization was already high at the beginning of the 15 year study) than the outer areas, where major urbanization processes started later.
The geographical variance derived using the GV method is shown graphically in Figure 10.The cumulative scale variation can be seen to reach a peak at Level 16 resolution, which is consistent with the tendency of the local variance analysis described previously.Since all of the large cities are known to be self-contained economic geographic areas undergoing rapid urban agglomeration with other cities, the very large geographical variance value at Level 16 resolution can be interpreted as indicating that the urban growth units in the Pearl River Delta area generally form large spatial patterns (Level 16 spatial resolution equal to 480 m coarse resolution).This may therefore be the appropriate spatial resolution at which to examine variations in urbanization velocity within the whole PRD area.For an additional city-level analysis, the whole of the PRD was further partitioned into city-based scenes.Their corresponding graphs of local variance were investigated in order to identify the moving window sizes that best matched urban growth units in the different cities.As shown in Figure 9, for the cities in the central urban agglomeration areas such as Guangzhou, Foshan, and Dongguan, the peaks in local variance were observed at low levels (ranging from Level 2 to Level 4), or at 60 m to 120 m fine resolution.In contrast, cities in the outer areas of the PRD yielded local variance peaks at higher levels (Level 8 to Level 32) or at coarser resolutions (480 m to 960 m).This means that the most appropriate spatial resolution for analyzing urbanization velocity may be different in the core areas (where urbanization was already high at the beginning of the 15 year study) than the outer areas, where major urbanization processes started later.Nevertheless, a city-based geographical variance analysis revealed spatial structure differences within the PRD. Figure 11 shows that the peak variances in the central part of the PRD are detected mainly at Level 1 to Level 4 resolutions.This area is also the main urban corridor of the most rapid urban sprawl [23].The cities in the outer areas of the PRD yielded the largest geographical variances, detected mainly at Level 8 to Level 32 resolutions.These results indicate the importance of choosing the appropriate resolution to use for different cities or groups of cities.The geographical variance derived using the GV method is shown graphically in Figure 10.The cumulative scale variation can be seen to reach a peak at Level 16 resolution, which is consistent with the tendency of the local variance analysis described previously.Since all of the large cities are known to be self-contained economic geographic areas undergoing rapid urban agglomeration with other cities, the very large geographical variance value at Level 16 resolution can be interpreted as indicating that the urban growth units in the Pearl River Delta area generally form large spatial patterns (Level 16 spatial resolution equal to 480 m coarse resolution).This may therefore be the appropriate spatial resolution at which to examine variations in urbanization velocity within the whole PRD area.Nevertheless, a city-based geographical variance analysis revealed spatial structure differences within the PRD. Figure 11 shows that the peak variances in the central part of the PRD are detected mainly at Level 1 to Level 4 resolutions.This area is also the main urban corridor of the most rapid urban sprawl [23].The cities in the outer areas of the PRD yielded the largest geographical variances, detected mainly at Level 8 to Level 32 resolutions.These results indicate the importance of choosing the appropriate resolution to use for different cities or groups of cities.Nevertheless, a city-based geographical variance analysis revealed spatial structure differences within the PRD. Figure 11 shows that the peak variances in the central part of the PRD are detected mainly at Level 1 to Level 4 resolutions.This area is also the main urban corridor of the most rapid urban sprawl [23].The cities in the outer areas of the PRD yielded the largest geographical variances, detected mainly at Level 8 to Level 32 resolutions.These results indicate the importance of choosing the appropriate resolution to use for different cities or groups of cities.

Discussion
Our research has primarily addressed the spatial resolution issue when measuring urbanization velocities and urbanization patterns.This research mainly investigated the effect that multiple grid sizes have when defining the geographical extent of urban growth through modeling their spatial autocorrelation and spatial variability effects at different resolutions.

Scale-Independence in the City-Based Urbanization Velocity Analysis Results
Statistical analyses (including the Pearson correlation analysis) were used to investigate urbanization velocity analyses at different spatial resolutions for scale-dependence or independence from a mathematical perspective.Ten different window sizes (or spatial resolutions) showed moderately strong correlations between socio-economic indicators and urbanization velocity, yet remained statistically independent of spatial resolution.These results confirm that the GDP, Total value of retail sales and Total number of jobs indicators can be very useful for measuring the rate and extent of urban growth [82].However, it is worth to note that although the integration of urbanization information and data from socio-economic sources is vital for accurate definition of political urban area units by urban planners, these three indicators only may not fully explain the driving factors of urbanization process.Other typical socio-economic factors, including GDP per capita, Income per capita, Deposit balance and Population, did not reveal significance in respect to the multi-level urbanization velocity results.While the three socio-economic indicators selected cannot comprehensively represent the driving factors of urbanization, the associations between the urbanization velocity results with other urban development indicators, such as transportation network, infrastructure facilities, etc., will be taken into consideration in future studies.The multiple spatial resolution urban velocity analyses in this research have also shown that those cities with above-average urbanization velocities are always conjunct with the core area of the Pearl River Delta and the extent of central agglomeration areas was not significantly affected by changes in spatial resolution, challenging the notorious modifiable areal unit problem (MAUP).The urban

Discussion
Our research has primarily addressed the spatial resolution issue when measuring urbanization velocities and urbanization patterns.This research mainly investigated the effect that multiple grid sizes have when defining the geographical extent of urban growth through modeling their spatial autocorrelation and spatial variability effects at different resolutions.

Scale-Independence in the City-Based Urbanization Velocity Analysis Results
Statistical analyses (including the Pearson correlation analysis) were used to investigate urbanization velocity analyses at different spatial resolutions for scale-dependence or independence from a mathematical perspective.Ten different window sizes (or spatial resolutions) showed moderately strong correlations between socio-economic indicators and urbanization velocity, yet remained statistically independent of spatial resolution.These results confirm that the GDP, Total value of retail sales and Total number of jobs indicators can be very useful for measuring the rate and extent of urban growth [82].However, it is worth to note that although the integration of urbanization information and data from socio-economic sources is vital for accurate definition of political urban area units by urban planners, these three indicators only may not fully explain the driving factors of urbanization process.Other typical socio-economic factors, including GDP per capita, Income per capita, Deposit balance and Population, did not reveal significance in respect to the multi-level urbanization velocity results.While the three socio-economic indicators selected cannot comprehensively represent the driving factors of urbanization, the associations between the urbanization velocity results with other urban development indicators, such as transportation network, infrastructure facilities, etc., will be taken into consideration in future studies.The multiple spatial resolution urban velocity analyses in this research have also shown that those cities with above-average urbanization velocities are always conjunct with the core area of the Pearl River Delta and the extent of central agglomeration areas was not significantly affected by changes in spatial resolution, challenging the notorious modifiable areal unit problem (MAUP).The urban agglomeration process in the PRD consequently seems to be unaffected by spatial resolutions, being evident over a wide range of spatial resolution (with resolution (pixel sizes) ranging from 30 m resolution to 15,360 m resolution).

Spatial Variances of the Multi-Resolution Urbanization Velocity Analysis Results
This research outlined three complementary approaches (empirical variogram, local variance, and scale variance) that can be used to represent spatial heterogeneity of the urbanization velocity analysis results.The spatial heterogeneity of urbanization velocity analysis results can be considered to be essential for validating the spatial extent of urbanization and for the process of choosing a suitable modeling scheme.It may be dependent on the spatial size of the surface structures and may be useful for optimizing field programs [53,78].
The empirical variogram models showed that the spatial autocorrelations varied with the grid size (i.e., spatial resolution), as a result of the different rates of loss of spatial variability for different spatial resolutions.The results show that a relatively coarse image resolution (such as the 240 m or 480 m resolutions) is able to capture more spatial variability in urbanization velocity than the original fine spatial resolution (30 m), as well as exhibiting a higher degree of autocorrelation.This research therefore concludes that relatively coarse resolutions (such as 240 m or 480 m) are already able to delineate urbanization velocity in the PRD.The degree of spatial autocorrelation and spatial variability of urbanization velocity estimates in the PRD can be used as references for other studies, although not in absolute terms, as the positional accuracy, spatial consistency and content efficiency for urbanization processes may be specific to the particular study area.
The results have confirmed the theory that the level of spatial autocorrelation in the variogram model is higher at coarse spatial resolutions (i.e., 3840 m resolution (Level 128) to 15,360 m resolution (Level 512) in this research; see also [53,77]), deviating significantly from UV analysis results at the original smaller scale.Nevertheless, at some fine and coarse spatial resolutions (i.e., 120 m resolution (Level 4) to 1920 m resolution (Level 64), in this research) variogram models failed to detect spatial autocorrelation in the model results when the landscape variability exceeded the modelling scale.This means that the variogram method is constrained by the image resolution.Some previous researchers have stated that descriptions of spatial autocorrelation and variability in an entire study area may not be sufficient to quantify the local spatial variability in sub-areas [49,54].Meanwhile, this research only compared the omnidirectional variograms that represented simultaneously the horizontal, vertical, and two diagonal directions of the multi-resolution pixel grids in this study, the pattern of urbanization velocity estimates that might all growth in a certain direction might also influence the spatial autocorrelations heterogeneities.It is therefore still necessary to consider strategies for combining different orientations of empirical variogram models with the specific spatial extent of an image subset, in order to quantify local variability and local spatial autocorrelation.Since analyses using different spatial resolutions have generally been gaining in popularity, this research also suggests that the variogram model should be included as an important component of the multi-scale spatial data investigations.
The results of the local variance and geographical variance methods are also particularly pertinent to understanding the "resolution of action" of the UV analysis results.The "resolution of action" is assumed to allow the identification of dimensions for particular patterns, which can subsequently facilitate an understanding of the underlying processes [51,78].When this research used moving windows of different grid sizes to obtain a hierarchical nested urban sprawl structure for the PRD, the graphs for local variance and geographical variance demonstrated similar scaling trends of the spatial structure characteristics.In other words, with the spatial resolution that yielded the maximum spatial variability in the geographical variance model coincided with the spatial resolution that yielded the highest likelihood of the neighborhood being similar in the local variance model.Graphs of the local variance and geographical variance models consequently implied that the within-class variance of urbanization velocity is an understandable, and even predictable, effect.Coarser spatial resolutions (i.e., 480 m) were better than the original resolution for delineating the urbanization patterns and processes within the PRD.
The empirical variogram model has proven useful in studying the spatial autocorrelation and optimal spatial resolution of remote sensing data in ecosystem [83][84][85][86].The local variance model is commonly used in defining an optimal size of geographical forest patterns [87], ecological biodiversity [88,89] and LULC patterns [61,62].The geographical variance model has also been employed to measure the scope and conceptual context of neighborhood [90][91][92] and landscape ecology [80].Nevertheless, these three methods are rarely tested in the grid-based urbanization spatial pattern analysis.There have been very few discussions as to whether any such grid-based urbanization patterns might depend on resolutions of observation.This research is the authors' first attempt to embed these experimental or measured variations models to quantify heterogeneous spatial variation in urbanization process.The components of three spatial variance models (i.e., plotted as a function of spatial resolution) are robust for describing hierarchical spatial structure features and spatial autocorrelation characteristics of urbanization patterns in a simple computation way and their results are easy to interpret.
Generally speaking, all three methods highlighted by this case study show that the influence of spatial resolution on spatial variance in urbanization velocity analyses at different resolutions is generally predictable: as the spatial resolution decreases (or the size of the moving window increases), the spatial variability and spatial autocorrelation gradually increase to a peak and then decrease again as the resolution decreases further.Although the spatial variability reaches a peak at a Level 8 resolution in the variogram model and at a Level 16 resolution in both the local variance and geographical variance models, all these three models (including the variogram model, the local variance model and geographical variance model) are necessary in order to systematically evaluate the properties of data regularization that shape grid-based urbanization phenomena.Depending on the characteristics of spatial autocorrelation and spatial structure variability in the multi-resolution hierarchical nested dataset, this research identified the coarse resolution 240 m (Level 8) spatial resolution to be appropriate for describing the variability in spatial structure and spatial autocorrelation of the urbanization velocity analysis results in PRD.This level was a suitable "resolution of action" for characterizing grid-based urbanization velocity within the PRD.Nevertheless, urban land types in different cities demonstrate different land use and land cover objects.The resolution variance diagnostic in this research would be insufficient to resolve the resolution selection issue in other areas.It is clearly demonstrated that it is difficult to identify changes in urban agglomeration processes using only one specific spatial resolution for different cities.The applicability of these spatial variance models should be further examined.

Conclusions and Outlook
Remote sensing data are widely believed to be constrained by the minimum spatial resolution of the sensor.This research has, however, shown that a resolution of 30 m in Landsat satellite data was more than adequate for a multi-scale urbanization velocity analysis over the whole of the PRD.Coarser resolution satellite data (>30 m) clearly presents new possibilities for multi-resolution urban growth analysis, particularly with regard to differentiating urban structures within cities.Nevertheless, fine to coarse resolution imagery data are indispensable for long-term observations of spatial distributions and fluctuations in built-up areas, and enable planners to track urbanization processes in a spatially and temporally explicit manner.Compared to other ways of grid-based urbanization velocity analysis utilizing built-up area information, which are based on using the variations of vegetation cover [33,93,94], impervious surface [95], night-time luminosity data [96], and land surface temperature data [97], are somewhat limited by the time-consuming nature of image classification.Spatial resolution can be used to refer both to the magnitude of a study (e.g., its geographic extent) and also to the degree of detail (e.g., its spatial variance) [77].Concepts of spatial autocorrelation and intra-pixel heterogeneity within the structure of multi-resolution data were thus used in this research as a basis for understanding the effects of spatial variability in a grid-based hierarchy.With the advantages of being simple, well grounded in theory, and largely compatible with multi-scale grid models, the empirical variogram, local variance and geographical variance models used in this research have provided an understanding of the nature and cause of spatial variability in satellite imagery.Nevertheless, investigations into the effects of scale on irregularly grouped hierarchical data now have access to a growing arsenal of techniques.Among them this research encountered various kinds of scale-based decomposition of variation.Further research is required into developing systematic procedures for extrapolating information from one scale to another.
Spatial resolution has emerged as a critical issue in the description of the hierarchical structure of complex urban systems.Remote sensing data and geographic technologies will undoubtedly advance the development of research into multi-scale urban geography.When exploring multi-resolution issues in urbanization velocity analyses, this research found that urbanization velocity analysis across multiple resolutions may yield accurate representations, but the specifications of feature characteristics may differ for different scales.Meanwhile, relationships between socio-economic indicators and urbanization velocity that are independent of scale may provide an initial insight into likely quantitative relationships between the spatial structures within urban landscapes and the social characteristics of urbanization.This may aid in the understanding of how social and economic factors, as well as population growth, might influence urban expansion and density within the PRD.The usefulness of grid-based urbanization velocity analysis may be critically affected by the nature of the grid size in the original remotes sensors.Nevertheless, the remote sensing-based conceptual distinctions and methodological guidelines regarding spatial resolution in this research can help resolve "the spatial resolution question" in urban geography.
chosen.The one-kilometer LULC data for 2000, 2005, 2010 collected from the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences (RESDC) served as auxiliary data for the land use classification through generating coarse resolution masks of LULC classification in PRD.High resolution imagery data (<5 m), which were available online from Google Earth, and urban footprint data with a spatial resolution of 30 m, obtained from the German Aerospace Center ([66]), were used as reference datasets to provide appropriate ground truth information to validate LULC classification over the PRD area.The original census data were obtained from the Guangdong Statistical Yearbooks for 2000, 2005, 2010, and 2014 (http://www.gdstats.gov.cn/tjsj/gdtjnj/).Census data from these four years (which are available for 27 original city-scale administrative units) were used to obtain the socio-economic indicators, including GDP, GDP per capita, Total value of retail sales, Total number of Job, Income per capita, Deposit balance, Population.This research employed these socio-economic indicators as the criteria used to investigate the effects of different spatial resolutions in urbanization velocity analyses.

Figure 2 .
Figure 2. Schematic evaluation of urbanization velocity analysis for different spatial resolutions.Figure 2. Schematic evaluation of urbanization velocity analysis for different spatial resolutions.

Figure 2 .
Figure 2. Schematic evaluation of urbanization velocity analysis for different spatial resolutions.Figure 2. Schematic evaluation of urbanization velocity analysis for different spatial resolutions.

6 L 1 L 2 L 3 L 4 L 5 L 6 L 7 L 8 L 3 L 1 L 2 L 5 L 8 L 7 LFigure 3 .
Figure 3.The hierarchy of moving windows that were employed when investigating urbanization velocity at different resolutions.

Figure 3 .
Figure 3.The hierarchy of moving windows that were employed when investigating urbanization velocity at different resolutions.

Figure 4 .
Figure 4. Location of samples selected in the PRD region in 2010.

Figure 4 .
Figure 4. Location of samples selected in the PRD region in 2010.

Figure 5 .
Figure 5. Urban growth velocities in the Pearl River Delta for the years 2000 to 2015.

Figure 5 .
Figure 5. Urban growth velocities in the Pearl River Delta for the years 2000 to 2015.

Figure 6 .
Figure 6.Major cities in the PRD with above-average urbanization velocities for the years 2000 to 2015.

Figure 6 .
Figure 6.Major cities in the PRD with above-average urbanization velocities for the years 2000 to 2015.

Figure 8 .
Figure 8. Local variance graph showing the multi-scale structure of UV analysis results (see text for explanations of the non-linearity of the curve).

Figure 8 .
Figure 8. Local variance graph showing the multi-scale structure of UV analysis results (see text for explanations of the non-linearity of the curve).

Figure 9 .
Figure 9. City-based analyses yield maximum local variances for each city in the Pearl River Delta and differentiate a core area from an outer zone.

Figure 10 .
Figure 10.Scale variance graph showing the multi-scale structure of UV analysis results yielding a peak for Level 16 (480 m resolution).See text for explanation.

Figure 9 .
Figure 9. City-based analyses yield maximum local variances for each city in the Pearl River Delta and differentiate a core area from an outer zone.

Figure 9 .
Figure 9. City-based analyses yield maximum local variances for each city in the Pearl River Delta and differentiate a core area from an outer zone.

Figure 10 .
Figure 10.Scale variance graph showing the multi-scale structure of UV analysis results yielding a peak for Level 16 (480 m resolution).See text for explanation.

Figure 10 .
Figure 10.Scale variance graph showing the multi-scale structure of UV analysis results yielding a peak for Level 16 (480 m resolution).See text for explanation.

Figure 11 .
Figure 11.The peaks of scale variance in different cities of the Pearl River Delta.

Figure 11 .
Figure 11.The peaks of scale variance in different cities of the Pearl River Delta.

Table 2 .
The overall accuracy of land use and land cover (LULC) classification in the Pearl River Delta for the years 2000 to 2015.

Table 3 .
Relationships between mean city-based urbanization velocity (UV) analysis results and socio-economic indicators across multiple scales.

Table 3 .
Relationships between mean city-based urbanization velocity (UV) analysis results and socio-economic indicators across multiple scales.

Table 4 .
Parameters of the fitted variogram models.

Table 4 .
Parameters of the fitted variogram models.