A comparative study on four major cities in Northeastern Thailand using urban land density function

Abstract Urban land density is an important factor to understand how cities expand. An “Inverse S-shape Rule” was implemented for the first time to analyze urban land density in Northeastern Thailand using the four cities Khon Kaen, Udon Thani, Nakhon Phanom, and Nong Khai as study sites. Land density function was tested using different data classification techniques from previous studies. Each city was investigated over two different time periods between 2002 and 2015. Declining pattern characteristics of metropolitan area density outward from city centers can be quantified by fitting the parameters to urban land density functions. An inverse S-shape function was identified as the best data fit. The four selected cities showed conventional density variation for decline in urban land area from city centers to outlying areas. Overall trend indicated that cities became more compact over time since the density differences between the urban core and urban fringe were greater with increasing infilling growth within the urban boundary. All four cities increased in size over time; however, the increasing amount of built-up land in the surrounding rural areas did not follow the same trend in each case. Some functional parameters required careful interpretation because of the linear shape of the city as in the case of Nakhon Phanom. Using highly detailed urban data resulted in lower densities of urban areas compared to the conventional pixel-based classification, and this affected the overall shape of the inverse S-shape function. The fitted parameters and their changing trends indicated that the urban land density function was useful for understanding urban form and urban sprawl in Thailand. Results can be used to develop a specific framework for other cities with similar attributes in the future.


Introduction
In recent years, the world has become more than onehalf urban for the first time in history (54.9% in 2017) and almost one quarter (23.9%) of the population with 500,000 or more now lives in urban areas (Demographia 2017). In 2008, the world demographic population dispersion altered due to urbanization with the registered city dwellers surpassing the rural inhabitants for the first time. Thus, urban planning has become increasingly important, with the growing demand for dwellings within limited land areas and available natural resources. Research regarding the planning of urban land usage has, therefore, intensified. One of the most studied topics involves urban growth modeling, especially for urban sprawl. Urban sprawl refers to the expansion of a city and its suburbs over more and more rural land as inefficient low density development. In built-up urban areas, a basic key parameter for quantifying urban sprawl considers the impervious features. Impervious land surfaces including concrete pavements, buildings, and other fabricated structures can result from uncontrolled development sprawl onto open and agricultural areas.
Database management systems are now employed to monitor and control urban sprawl through geographic information system (GIS) and remote sensing data. The characteristics of urban fragmentation and patch density can be quantified by GIS at landscape level (Sudhira, Ramachandra, and Jagadish 2004). Numerous measurements have been used in the past to quantify urban development functions. Estoque and Murayama (2013) analyzed the dynamics of Baguio City, Philippines to derive meaningful information using a relative entropy method. Pham, Yamaguchi, and Bui (2011) applied the percentage of like adjacency (PLADJ) metric to map urban sprawl in Hanoi, Nagoya, Hartford, and Shanghai, while the patterns of urban sprawl were characterized using the Shannon entropy (Sudhira, Ramachandra, and Jagadish 2004;Jat, Garg, and Khare 2008;Yue, Liu, and Fan 2013). However, the most widely used method in recent studies on urban patterns is the landscape metrics. Landscape patterns are now identified and depicted quantitatively using an extension tool in the ArcGIS program which is simple, convenient, and user-friendly. Useful data can be extracted from thematic maps derived from remote sensing (Herold, Goldstein, and Clarke OPEN ACCESS 2003;Deng et al. 2009). Satellites are now used as transport media for remote sensing equipment to record GIS temporal and spatial data and monitor patterns of urban sprawl (Huang, Lu, and Sellers 2007;Yu and Ng 2007;Huang, Wang, and Budd 2009;Nassar, Blackburn, and Whyatt 2014;Hu et al. 2015).
The study of urban density is imperative for a comprehensive understanding of how metropolises function since it plays a key role in understanding the spatial process of urban growth, which can lead to urban sprawl. A landscape where buildings are highly diversified is referred to as an urban sprawl with increased costs of providing public services and infrastructures. Hence, higher density cities are more sustainable than those with low density since the former conserve land for agriculture and environmental protection (Mohammed, Alshuwaikhat, and Adenle 2016). Concentric partitioning of cities was commonly used in studies on urban land density, with many researchers following a concentric ring pattern approach (Pijanowski and Robinson 2011;Ramachandra, Aithal, and Sanna 2012;Shi et al. 2012;Schneider, Chang, and Paulsen 2015). The negative exponential spatial distribution of population density in Bangkok was quantified by Murakami et al. (2005), who compared results generated by the Clark and Newling models depicting the diverse stages of urbanization. Jiao (2015) postulated a sigmoid function to represent the basic parameters that described urban land density. He applied this concept to compare a city at disparate time points, and also different cities at the same time horizon, and plotted the parameters to quantify and define the measurement of urban sprawl.
In Thailand, urban growth is dominated by the Bangkok urban area, with sprawling growth threatening conservation efforts from poor urban planning. Recently, predictions suggest that economic expansion and industrial growth will be concentrated in the future in the peri-urban areas and they are trending to be the secondary urban centers (Thebpanya and Bhuyan 2015), e.g. in the northeastern region which is a financial trading hub with Laos, China (Yunnan), and Vietnam. The ASEAN Economic Community (AEC), established in 2015 to promote cooperation and development between ten South-East Asian Nations, has greatly benefitted the northeastern region of Thailand (Chia 2013). However, more public infrastructures such as roads, highways, and railways need to be built, and the consequences can lead to increasing urban sprawl in many urban locations in the future, especially in major cities.
Four major cities in Northeastern Thailand were selected as study sites in this study. All the cities are in the early phase of urban growth from AEC trade, which is expected to expand rapidly and may cause urban sprawl. The urban land density function can help to understand urban form and urban sprawl quantitatively. Globally, this function has only been previously tested by Jiao (2015) in China. Therefore, experiments on different geographic locations and diverse urban area size are of interest for global comparison. Data using different classification techniques i.e. detailed digitized urban patches were collected to identify any advantages or disadvantages. Two research questions arose: (1) Does an "Inverse S-shape Rule" land density function help to understand urban form and urban sprawl in the four major cities in Northeastern Thailand? (2) Is the result significantly different from previous studies after applying more detailed urban input data? Land density functions were monitored to characterize the spatial patterns of urban growth and measure urban sprawl based on the Inverse S-shape Rule. This was the first test of the urban land density function developed by Jiao (2015) outside China. Results will benefit researchers working in this specific field, and also people responsible for decision-making regarding urban planning in Thailand.

Study area
Four cities selected for study are located in Northeastern Thailand (as shown in Figure 1). The study region comprised a flat landscape with the latitudes between 14°0′ N and 18°30′N and the longitudes between 101°0′E and 105°0′E on the Khorat Plateau. This area of 160,000 km 2 is bounded by the Mekong River to the north and west, by Cambodia to the southeast, the Sankamphaeng Range south of Nakhon Ratchasima Province, and the Phetchabun Mountains to the west. Four major metropolitan population centers are located within 30 km × 30 km. Borders towns with Laos have sprung up to cater for the international trade in goods and labor from as far afield as Vietnam and Southern China, resulting from the new established AEC. Previously, farming land surrounded the urban centers, but these arable areas are now covered by urban infrastructures. Land prices have spiraled, with ownerships transferred from the local inhabitants to traders and conglomerates outside the region. Large-scale development projects have sprung up, including construction of the Mittraphap Road as a major highway to provide access to the area as a significant expressway for trade advancement. New developments have expanded the economy of Northeastern Thailand and promoted rapid growth of the urban territories.

Defining urban extent
Here, urban land was assumed to have an impervious artificial surface as a built-up environment. This pretext delivers acceptable results for assessments using remote sensing data. The GIS database was created by on-screen digitizing of urban patches (impervious surface) using Google Earth as base map (QuickBird, World View II, III) with zoom level height of 300-500 m (scale 1:4,000). Test sites included two inland cities (Khon Kaen and Udon Thani) and two border cities (Nakhon Phanom and Nong Khai). The test area was based on a 30 km × 30 km boundary and a database was created using comprehensive plan boundaries as the lowest range of the digitizing area (as shown in Figure 2). The urban centers of each city were selected and created as point data. Jiao (2015) used the concept of concentric ring partitioning to define city limits and select urban centers with concentric 1-km buffer rings extending outward (as shown in Figure 3). The final outer ring circumference has to be large enough to encompass the continuous density of the outwardly trending metropolis while excluding small satellite urban microcosms far from the central urban core. Urban land density was calculated by the area of urban land divided by the total land area of each concentric ring partitioning. Water bodies were excluded from where g(x) represents the density function as (x, x + Δx), Δx refers to a radial distance interval (buffer distance), CH Δx records the total land cover changes associated within the area between distance interval (x, x + Δx), and NCH Δx quantifies the total area of land that could be developed in the same distance interval.

Inverse S-shape Rule
The overall urban spatial distribution patterns can be determined by plotting urban land density against the distance from the city center. In general, urban density decreased outward from the city center. The spatial variation of urban land density is described by a modified sigmoid function with an inverse S-shape (Jiao 2015).
(2) g(x) = CH Δx CH Δx + NCH Δx the calculation as they are not used for urban development. Nong Khai and Nakhon Phanom are located by rivers and these large waterbody areas were excluded from the total area of each ring because the difference of urban land density calculation after subtraction was larger than 1% (Jiao 2015). Urban changes over two time intervals were created as dependent variables. The time intervals were different between cities due to the availability of satellite imagery in Google Earth. Details of data used in this study are shown in Table 1. The urban growth map over different time periods for Nakhon Phanom is shown in Figure 4.

The classic urban land density function
City planners use the term 'urban density' to describe populations in selected metropolitan districts. Distance intervals radiating outward from town centers provide a grid to quantify population density. Generally, most cities show a negative exponential density decrease with increasing distances from their centers, which can also be applied to development units such as buildings (Equation(1)) (Goodchild 2000;Brunsdon 2001;Cheng and Masser 2003).
where x represents an individual buffer distance from the city center or central business district (CBD), and λ is the density gradient. The density calculation relating to each buffer distance is shown by Equation (2).
(1)  iterations to refine the parameters. MATLAB 2008a was employed to implement the statistical analysis using a nonlinear regression model to determine the land density function parameters. The Hougen-Watson algorithm was selected as the best choice for linear least squares fitting. The nlinfit functions in Statistic Toolbox were used to estimate the least square parameters.

Measuring the compactness of cities
In order to analyze urban variation, the derivatives of the fitted curves of the land density functions must be examined. The first derivative (denoted by f′(r)) of the plot of the density function was analyzed to measure rates of decrease and variation, and thereby deduce the coordinates of the fastest decrease in the graph of the function curve (Jiao 2015). From the first derivative, it can be inferred that the coordinates of the fastest decreasing point (denoted by p 0 in Figure 5) in the curve of the urban land density function is: where f is the urban land density, r is the distance to the urban center, e is Euler's number, and α, c, and D are constants. The constant α is the parameter that controls the slope of the curve of urban land density function.
A compact city has a large α value while a city with a small α value indicates a low-density or dispersed city. The constant c represents the background value of built-up land density in the hinterland of the city. When c increases, it indicates that construction on land in the surrounding rural area increases due to the development of the entire region. Finally, parameter D is an estimation of the radius of the main urban area of a city and when increases with time for a single city, describes the expansion of an urban area.
Here, a nonlinear least squares method was postulated to derive city land density functions. A non-linear function was fitted to the collected data after successive (3) f (r) = 1 − c 1 + e ((2r∕D)−1) + c As stated by Jiao (2015), k s is useful to compare urban densities among cities of similar size because k s is determined by three parameters in which D is a changeable variable. Thus, a higher D leads to a lower k s . The k s value tends to decrease with urban growth, and large cities tend to have lower k s values. Another index, proposed to avoid bias and represented by k p , can be written as: A compact city has a small k p value. If k p is large, then the city has low density and is highly dispersed. Since the k p value is only related to the constant α, it is a suitable index for characterizing the changes in the degree of dispersion of diverse cities with time. Figure 6 shows the plots of the urban land density functions. There was a decline in the urban land area from city centers to outlying areas. This finding was similar to Jiao (2015) and mirrored conventional understanding of metropolitan density variation. Table 2 indicates that the constant D increased with time in all cities. In terms of city size when implied from the value of constant D, Khon Kaen had the largest range of urban boundary in 2015 with urban radius at 3.36 km. In terms of rate of expansion, the difference of the constant D can be used to estimate the change. For example, during the first period, Nong Khai experienced the highest change in urban growth boundary, increasing city radius by more than double from 0.41 km to 1.04 km (60.6%) with highest annual rate change of 7.6%. The amount of land construction in the surrounding rural area due to the development of the entire region can be speculated by the value of constant c. This constant c is the horizontal asymptote of the function (Equation (3)); as r approaches infinity (i.e. rural area), f(r) approaches c. The value of constant c in all cities increased over time except for Khon Kaen. The graph of urban land density of Khon Kaen (Figure 6(a)) also shows the crossover of the fitted line, indicating a decreasing amount of the c value. The slope of the curve representing urban land density is controlled by the parameter α. Α high value of α indicates rapidly decreasing land density from the city center outward to the surrounding areas, hence representing a compact city (Jiao 2015). Results showed that Nakhon Phanom had the highest average value of α at 1.1. Results also showed that values for constant α increased with time for all cities, except for Nakhon Phanom during the second period, suggesting that in most cities, differences of density the between urban core and urban fringe become greater over time. In other words, cities become

Physical interpretation of the parameters
The second derivative (denoted by f′′(r)) has two extremes which locate where the rate of decrease changes most greatly. The coordinates of the two corresponding points, p 1 and p 2 , on the curve of urban land density, can be calculated by Equations (5) and (6), respectively.
The three points, p 0 , p 1 , and p 2 , can be used to partition the urban area in a concentric manner: the urban core, the inner urban area, the suburban area, and the urban fringe ( Figure 5). The partitions suggest that the urban core with high density of urban land is slowly decreasing while the inner and sub urban area have relatively rapid decreasing rate of urban land density. The urban fringe, however, has slow rate of decreasing in urban land density due to the expanded environment surrounds the main urban area. Since the curves of urban land density for compact cities are steep, but those for sprawling cities decrease slowly. Two points, p 1 and p 2 , are used to define the slope of urban land density function. The slope, denoted by k s , is calculated by: surrounding the core. The higher the value, then the steeper the curve, i.e. a more compact city, while sprawling cities have low k s value. Khon Kaen and Udon Thani had very similar k s values, much lower than Nong Khai and Nakhon Phanom which indicated that both cities had a lower overall urban land density, and possibly experienced more expansive urban growth. A second index is the k p value. This counteracts bias among cities with different sizes. A compact city will have a small k p value. On the other hand, a city with a large k p value is a low-density or dispersed city. Results show that Nakhon more compact. However, in the case of Nakhon Phanom, which had the most linearly spatial structure among the four cities, further discussion of the non monocentric cities is provided in Section 4.3.

Measuring the compactness of cities
The k s and k p values of the four cities are shown in Table  2. The k s value is an indicator that depicts the compactness of a city, which usually has high density in the core area and decreases to a lower value in the areas Figure 6. the fitted graphs of urban land density functions using the 1-km buffer distance. that urban growth happened toward urban cores which made the city's spatial structure denser and more compact. However, this conclusion can only be applied to cities with recorded urban growth direction. In the case of Nakhon Phanom, which had the most linear spatial structure among the four cities, the constant α required careful interpretation. Results of land density function found in Nakhon Phanom showed many similarities to the two linear cities, Lanzhou and Yinchuan, found in the work of Jiao (2015). For example, the values of determinant coefficients (r-square) of the three cities were the lowest among the tested cities. Lanzhou and Yinchuan had the lowest r-square (0.93 and 0.96, respectively) compared to the other cities (0.97-0.99), while Nakhon Phanom had the lowest r-square at 0.96 compared to the other three cities at 0.97-0.99. This indicated that actual data points were less fitted among the other tested cities. Such variations in data might result from asynchronous spatial structure between the linear form of the cities and the buffer partitioning based on monocentric cities. Similar to Lanzhou and Yinchuan, Nakhon Phanom showed a sharp decrease in urban land density from the urban core, which resulted in a high value of the constant α derived from the fitted function. Density outside the urban core could have been underestimated. The circular ring partitioning expands outward asynchronously with linear expansion of the city, thus, the actual city could be more dispersed than the fitted function suggested.

Conclusions
Four major cities located in the northeastern region of Thailand were selected as samples for land density function analysis over two time periods during 2002-2015. Results indicated that data from all the urban areas studied fitted well with a nonlinear least square function, and the variation of urban land density followed an "Inverse S-shape Rule. " Results also suggested that all the studied cities were expanding and becoming more compact overtime. Findings clearly showed that detailed urban area derived from high-resolution images gave lower density values and affected the overall shape of the inverse S-shape function. Urban area should be classified from high-quality Landsat images to achieve accurate results, especially for smaller cities used in this experiment. Quantitative analysis of the inverse S-shape of urban land density implied that specific zoning/regulation in urban land use planning plays an important role in the final results and evolution of parameter values of the density function. More experiments are required to confirm these quantitative findings of the inverse S-shape rule for other cities with different spatial structures such as polycentric cities. Statistical analysis, including the ratio of the growth rate to measure the degree of urban sprawl, could also be assessed in future research.
Phanom was the most compact city overall, while the most dispersed city was Nong Khai. Interestingly, in Nong Khai during the first period which was less than 10 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011), the k p value dropped by more than 100%, including that the city significantly increased in compactness during this time period.

Spatiotemporal change of urban structure
The evolution of land density function parameters at two time points (Table 2) can be discussed in relation to spatiotemporal change of urban structure of the cities as follows: Table 2 indicate that all four cities have expanded in size. The constant D is an estimation of the radius of the main urban area of a city. Jiao (2015) showed that every city had increasing values of D over time. However, unlike the exponential growth found in Chinese cities by Jiao (2015), observed growth rates of D for the four cities in Northeastern Thailand implied urban growth in a more linear fashion. Interestingly, Udon Thani was the only city during the second period that showed an increasing annual rate of expansion. This result can infer that Udon Thani is the most prosperous in terms of future investments.

The constant c
Overall results of constant c in all four cities were less than 0.1.This was consistent with results of 28 major cities in China studied by Jiao (2015). Each city showed an increasing trend of the constant c, except for Khon Kaen, indicating that while construction on land in the surrounding rural area of other cities increased, due to the development of the entire region, Khon Kaen became more densely built up at the urban core area. This could be the result of land use planning over certain areas of urban growth inside the city boundary.

The constant α
Overall results showed a much lower value for α than found by Jiao (2015). This discrepancy may be due to differences in size of the cities chosen for the study. Most cities in China (Jiao 2015) had much larger urban areas and were denser than the four cities in Northeastern Thailand. Another important reason was due to differences in the input data. Jiao (2015) used urban impervious surface data extracted using the Maximum Likelihood algorithm. Our study used an onscreen visual digitization approach, resulting in more detailed urban patches. Many gaps occurred within the urban core area as a result of using this finely detailed data, with reduced differences in density gradient from the urban core to the surroundings which reduced the value of the constant α. Increasing trends of the constant α were found in all cities except for Nakhon Phanom. This suggested