Rank-size Distribution of Cities and Municipalities in Bangladesh

: This paper examines and updates the rank-size distribution of cities and municipalities in Bangladesh between 1990 and 2019 based on two criteria: (1) built-up urban areas; and (2) population. The distribution of built-up urban areas and population are compared to provide a robust theoretical underpinning of Zipf’s law for future urban developmental planning framework. The data on built-up urban areas is extracted from land cover classification using Google Earth Engine and the population data is obtained from the decennial censuses. The comparison of the conformity to Zipf’s law indicated contradictory results. While a greater proportion of the population has been increasingly concentrated in the smaller and midsized cities over the last three decades, built-up urban areas, on the other hand, have been mostly clustered in two largest cities— Dhaka and Chittagong—accounting for 50 to nearly 60 percent of the total built-up urban areas. These results shed light on the magnitude of continued spatial inequalities in urban development amongst cities and municipalities in Bangladesh despite there being an overall increase of evenness in the distribution of population over time. These results imply an unsustainable rate of urban expansion in Bangladesh and reinforce the need for the exploration of policies and regulations targeted at guiding the rate and direction of evenness in urban expansion.


Introduction
Bangladesh has experienced significant growth in all segments of the community, accompanied by built-up urban areas since the 1970s [1,2]. Less than nine percent of the population in Bangladesh was urbanized in 1974, but today 37% of the country's population live in urban areas, and this amount is expected to rise to 57% by 2050 [2]. New infrastructural developments have arisen to provide housing, roads, schools, services and parks for this growing population, but mostly with minimal regulation and planning. As a result, there have been haphazard urban buildups at the expense of productive agricultural land [3,4], vegetated land [5] and fragile wetlands [6], which are often encroached upon illegally. The uncontrolled and unrestricted urbanization lead to the transformation of natural landscapes into built-up urban areas, and their impact on the local environment and overall ecosystem of the region-raised questions about the sustainability and livability of cities in Bangladesh.
While urbanization has generally been portrayed as a sign of economic growth, with all others environmental adverse impacts being considered as a growing pain, this concept has been challenged by Bloom et al. [7]. There is no evidence of close links between the level of urbanization and the rates of economic growth in developing countries like Bangladesh. Even though the GDP per capita has been increasing in the country for the last three decades or so, and Bangladesh may soon be considered as a middle income country [2], there has been a lack of subsequent socio-economic improvement in the distribution of benefits across the geography and the population [8][9][10][11][12][13]. A well-balanced and well-informed urbanization is unquestionably key to maximizing the benefits and enhancing the quality of life for everyone while minimizing environmental degradation and other potential adverse impacts of a growing number of city dwellers [2].
The urban system is complex. The complexities within the urban system are deeply rooted in the regularity that is seen in the ordering of the urban areas in a variety of spatial scales, regardless of the region [14,15]. This consistency, presented by the hierarchy of urban systems, reflects the economic functions and the economies of scale that drive their growth and can further inform the pressure of the resources needed to serve these demands in a given area [14,16]. Therefore, the application of theories on city size, while underexplored within the context of shaping of urban policy, can provide a much needed additional dimension when exploring the questions of the sustainability and livability of cities, as well as constructing related policies [17].
The primary goal of this study is to examine and update the rank-size distribution of urbanized areas of Bangladesh based on the most recent data available (between 1990 and 2019) by measuring built-up urban areas and population. Rank-size rule distributions can be used simultaneously to support policies targeting smart development and the equitable distribution of resources needed in urban areas of Bangladesh. Previous studies [18][19][20] have examined the rank-size distribution of urbanized areas for Bangladesh using different data approaches-population census data [18,19] or satellite data to define natural cities derived through the identification of areas of human activities and urban agglomeration [20]. Hence, a direct comparison between rank-size distributions of these studies has been difficult, particularly due to their use of different urban units for the analyses. In this study, to fill that gap, we establish a common urban unit based on the existing administrative boundaries and compare the conformity of the distribution of built-up urban areas and the distribution of population within these urban units in Bangladesh to Zipf's law. Due to the lack of more recent up-to-date land cover data spanning the entire Bangladesh for different time periods, we use Google Earth Engine (GEE) to create a land cover classification and extract the built-up urban areas within the urban units and census data to extract the population living in these urban units. A direct comparison between the distribution of the two metrics provides valuable insight on the expansion of the urban footprint in reference to the increase in urban population and the sustainability of future urban development of the region. Understanding this urban distribution and the impetus behind these distributions is particularly important for developing countries such as Bangladesh in order to contribute towards a wide range of policy debates, including the optimal use of economic resources based on the urban hierarchy [21].

Zipf's Law
The hierarchical relations of city sizes and urban systems have been previously examined in the works of Christaller [22] through the center place theory, that looks at the size and distribution of settlements, and Lösch [23], who extended the central place theory by integrating economic and mathematical components in the theory. Similarly, Zipf's law is one of such theories examining the relation between the hierarchy of cities and their respective sizes [24]. Zipf's law, also known as the rank-size rule, is one of the prominent methodologies that have been implemented in the study of urban geography and examination of the city size distribution and growth of cities [25]. The initial concept for Zipf's law was put forward as a special case of rank-size law by Auerbach, a German geographer, in 1913, who indicated that the distribution of the city sizes within a given region follows a Pareto distribution [26]. This general rank-size law was revised and refined most notably by Zipf, hence being named Zipf's law [19]. While Zipf's law is frequently used in reference to the size of the cities, it was implemented by Zipf [27] to examine word frequency and has been used in a wide range of topics, including the size and income of companies [28], the expression of genes [29], the use of internet [30] and a variety of natural and physical phenomena [31].
Mathematically, the Pareto distribution can be represented by the equation: In the logarithmic form as: where i represents the city taken into consideration, S i represents the population of the city i and R i represents the ranking of the city i. A is the constant, with b representing the Pareto exponent [24]. Equation (1) was restated by Zipf as Equation (3). Particularly when taken in terms of city-size states, the absolute size of the city is inversely proportional to the rank-size of the city. In other words, in the law states that are based on the ranking of the cities, the largest city will be roughly two times the size of the second biggest city, three times the size of the third biggest city and so on, or "the size of the nth city is approximately one-nth of the size of the largest city" [32]. This relationship between the rank-size relation is also commonly referred to in the literature as the "Power law" [33]. Mathematically: where i represents the city taken into consideration, S i represents the population of the city i and R i represents the ranking of the city i. K is the constant, and q represents the Zipf's exponent [24]. This rank-size can be expressed in the logarithmic form as: Here, if the value of q = 1, the rank-size distribution is said to follow the Pareto distribution and satisfy Zipf's law. This suggests that within the urban system "the concentration ratio of population distribution and the dispersion ratio of population are completely equal" [34]. If the value of q < 1, it suggests the population is relatively evenly distributed, with the urban population being more concentrated in the medium to smaller sized cities compared to the larger cities. Finally, if the value of q > 1, it suggests the population is unevenly distributed, with bigger cities hosting a larger portion of the population as compared to smaller and medium sized cities [34].
Additionally, Zipf's law, which is the static form, can be expressed dynamically as Gibrat's law, which indicates that the growth of the city does not depend on the initial size of the city [35]. The formal definition of Gibrat's law states that the mean and variance of the distribution function are independent of the initial size of the city; instead, they depend on external functions such as the industrial activities within the city [25,36]. Modica et al. [15] indicates that previous studies by Champernowne [37] and Simon [38] have argued that that rank-size distribution arises if Gibrat's law is satisfied. Gabaix [25] further suggested that once the drivers of city growth satisfy Gibrat's law, the distribution will converge to Zipf's distribution.
A generalization of Gibrat's law proposed by Córdoba [39] suggested an allowance for city sizes to affect the variance of growth for the cities but not their mean growth rate considering Gibrat's law [15]. This relaxing of Gibrat's law on non-proportional variances allowed for the examination of a divergence in Zipf's coefficient (Equation (3)) away from 1 [39]. Modica et al. [15] outlined the relationship between Zipf's law and Gibrat's law as follows: if Zipf1's law holds (i.e., Zipf's coefficient is 1), for Gibrat's law to apply, the mean and the variance should not depend on the size. If the rank-size distribution is more unevenly distributed (as given by Zipf's coefficient > 1), for Gibrat's law to apply, the mean needs to be independent of the size but the variance does not, indicating a greater volatility in growth for smaller cities. Finally, if the rank-size distribution is more evenly distributed (as given by Zipf's coefficient < 1), for Gibrat's law to apply, the mean needs to be independent of the size but the variance does not, thus indicating a greater volatility in growth for larger cities as compared to smaller cities.
Mathematically, Gibrat's law can be expressed as [15,40,41]: Here, i represents the city taken into consideration, t represents the time, S represents the size of the city, ε i represents the independently and identically distributed error term and, β represents the coefficient. If the value of β = 1, Gibrat's law holds. If the value of β > 1, the size of the city is seen to diverge from the mean, and the growth is expected to be greater in larger cities. If the value of β < 1, the size of the city is seen to converge to the mean, and the growth expected is smaller for larger cities [15].

Literature Examining Zipf's Law in Developed Countries
There have been a large number of empirical studies in both developed and developing countries using Zipf's law to examine the distribution of urban areas and urban hierarchies. For the US, frequently cited studies-such as those by Krugman [42], using population data from 130 metropolitan areas, Gabaix [25], using population data for 1990 form the 135 largest metropolitan statistical areas, and Berry & Okulicz-Kozaryn [40], using a regionally-integrated economic area with a population greater than 500,000-have concluded that Zipf's law holds for these urban regions. Internationally, Giesen and Südekum [43] used population data of large cities in Germany and concluded that the city rank-size distribution followed the Zipfian power law, while Budde and Neumann [44] used a varying definition of urban region within Germany based on the population density across 1 km-square grid. These results indicated a moderate but systematic deviation away from Zipf's law when areas with lower density were included in the definition of urban regions [44]. Rastvortseva and Manaeva [45] used 2014 population data to examine the rank-size of the Federal Districts in Russia and concluded that the Russian territory conformed to Zipf's law. However, for individual Federal Districts in Russia, the conformity to Zipf's law depended on the size of these Federal Districts. On the contrary, for countries such as Canada, using the 152 largest urban areas [46], and Poland, using urban communities as their urban units [47], the distribution of population was seen as not conforming to Zipf's law and was observed to be skewed towards the larger cities.

Literature Examining Zipf's Law in Developing Countries
The application of Zipf's law in developing countries has largely indicated a non-conformity of the rank-size distribution of urban areas. The results from the study by Gangopadhyay & Basu [48] for India, using census data from 1941 to 2011 for cities, indicated that for the upper tail of the cities, with a population greater than 120,000, the city size distribution followed Zipf's law. However, other studies from India by Luckstead & Devadoss [49], examining the rank-size distribution of cities from 1950 to 2010, also using census data, indicated that the cities exhibited a log-normal distribution instead of a Pareto distribution, particularly after the economic liberalization of India in 1980, resulting in an increased rate of urbanization. Arshad et al. [50] used population data of five census years, from 1951 to 1998, to examine the city rank-size distribution at the national and provincial levels of Pakistan. The results of this study concluded that the city rank-size distribution followed Zipf's law at the provincial level, but not when these cities were examined at a national scale [50]. Tuholske et al. [51] used OpenStreeMap data augmented with synthetic gridded population data to examine the change rank-size distribution of 4750 individual urban settlements across the continent of Africa from 2000 to 2015. The result resonated with the findings from a previous study by Jiang et al. [20], which indicated a noticeable deviation from the Pareto distribution and concluded that the African urban settlements did not follow Zipf's law [51].
In other cases, the conformity of urban areas within the given areas of interest to Zipf's law has been less concrete. For example, the examination of Zipf's law in Malaysian cities was performed using population census data for five points in time, from 1957 to 2007, by Soo [52]. The results indicated a rejection of Zipf's law for all years, except 1957, when the full population sample was considered. Over the years, the distribution of the city sizes moved away from the Pareto distribution and towards a greater degree of unevenness [52]. Alternatively, when a truncated sample of only larger cities was considered, this trend was reversed, with an increased evenness in the distribution of population among the larger cities [52]. For cities in Morocco, the adherence to Zipf's law was examined based on Sustainability 2020, 12, 4643 5 of 26 three truncated population census data at three different years, from 1982 to 2004 [53]. The results from this study concluded that even though the Pareto coefficient for the cities at these levels was not statistically different from 1, the urban system for Morocco tended towards a more balanced distribution over time [53].
In the case of China, the empirical validity of Zipf's law for Chinese cities using census data from 1990 and 2000 has been examined by Gangopadhya & Basu [48,54] at different levels of sample size truncation, at different time frames and using various definitions of urban areas. The results indicated that the city rank-size distribution of population in China followed Zipf's law, but only when the biggest cities in terms of population were considered, showing an increase in the Pareto coefficient over time. This increase in the Pareto coefficient was further confirmed by Peng [55] and Ziqin [56], both suggesting an increase in the concentration of population in the smaller and medium sized cities, leading to an evenness in the distribution of population in the biggest cities. Q. Huang et al. [57] used naturally delineated cities or natural cities based on nighttime light data to identify their urban footprint and examine its conformity to Zipf's law [57]. The results indicated that the Pareto exponent showed a gradual increase from 1992 to 2008, showing an evenness in the distribution of urbanized areas [57].
Examinations of Zipf's law have also been conducted at the regional and global scales. These large-scale investigations of the rank-size distribution of cities have included studies by Rosen and Resnick [21] using population data from 1970 for cities in 44 countries. The results from the study calculated the average Pareto coefficient to be 1.136, not conforming to Zipf's law [21]. More recently, a large-scale investigation of Zipf's law was conducted by Soo [19] for 73 countries using population census data and based on the definition of urban areas as cities and urban agglomerations. The results derived for cities defined as urban areas was consistent with the previous study by Rosen and Resnick [21] and rejected Zipf's law. The results differed when the definition of urban areas as urban agglomerations was used. Furthermore, using the Hills estimator instead of the Ordinary Least Square (OLS) method, the results were seen to be reversed [19]. A meta-study conducted by Nitsch [58] using 29 studies and 515 Pareto estimates from a wide range of territories and time period concluded that Pareto estimates significantly deviated and were higher than 1. This suggested that the cities around the world, on average, were seen to be more evenly distributed. Nighttime light data was used by Jiang et al. [20] to determine if Zipf's law holds at a global scale, including a total of 30,000 cities around the world. This research concluded that Zipf's law held well for natural cities at a global scale and at a continental scale (except Africa), but showed considerable variation at the country scale and during different time frames.
Overall, Arshad et al. [24], surveying the literature, divided the findings from the empirical studies related Zipf's law and rank-size distribution into five categories: 1. Studies that show urban areas within the selected area of interest conforming to Zipf's law; 2. Studies that reject the distribution of urban areas within the given areas of interest conforming Zipf's law; 3. Studies that have shown a mixed result when examining the distribution of urban units within those following Zipf's law; 4. Studies that show the distribution of the urban areas moving away from Zipf's law as the country experiences urbanization; 5. Studies that have shown that city size distribution moves away from Zipf's law over time.

Literature Examining Zipf's Law for Bangladesh
There have been several studies examining Zipf's law for Bangladesh (Table 1). Winidowa [18] used national census data for Bangladesh from 1961, 1974 and 1981 for the analysis. In this study, areas identified as towns in the census were taken as the primary unit areas of study and included 77, 108 and 451 regions that were defined as towns for the years 1961, 1974 and 1981, respectively [18]. The results indicated that the rank-size distribution for the towns in Bangladesh increasingly did not conform to Zipf's law, showing a gradual increase in the Pareto coefficient over the years, going from 1.0296 in 1961 to 1.1526 in 1974 and to 1.2654 in 1981. This indicated bigger towns in Bangladesh hosting a greater proportion of the population. Winidowa [18] describes these deviations away from Pareto distribution as a function of "political divisions, the low degree of urbanization and a weak economic development level of the country" (p. 200). In other global studies associated with Zipf's law, the Pareto Sustainability 2020, 12, 4643 6 of 26 coefficient reported for Bangladesh has shown a large degree of variation. Soo [19], using population data from 1991 for 79 cities in Bangladesh, calculated the Pareto coefficient at 1.0914, not significantly different from 1; while using urban agglomeration data, with a total of 43 agglomerations, the Pareto coefficient was calculated at 0.8068 and was shown to be significantly different from 1 [19]. Similarly, D'Costa [59] examined the rank-size distribution of urban centers in Bangladesh from 1901 to 1981 using population data. The results indicated a largely even or "deconcentrated" distribution of population for a majority of the time frame, with an increase in primacy from 1961 to 1971. This was largely attributed to Pakistan's policies of developing major cities such as Dhaka and Chittagong, formerly belonging to East Pakistan, following its independence from India [59]. Natural cities derived based on nightlight data with areas greater than 10 sq. km.
The rank-size distribution of natural cities showed an uneven distribution, with a Pareto coefficient greater than 1 for all three years Nishiyama et al. [60], also using urban agglomeration data from 1991, but with a total of 140 agglomerations, calculated the Pareto coefficient at 1.27 and found the value to be significantly different from 1 at one percent. In a more recent study, Jiang et al. [20] used nighttime light data to derive natural cities for 1992, 2001 and 2010, with natural cities being defined as those with areas greater than 10 square kilometers. The result indicated a Pareto coefficient of 2.1 for 1992, 2.2 for 2001 and 2.1 for 2010 [20]. Currently, however, there is no study examining the city rank-size distribution using the most recent data for Bangladesh. The primary goal of this study is to examine the conformity of built-up urban areas and the population in these urban areas to Zipf's law and compare the rank-size distribution of the two.

Some Considerations Regarding Zipf's Law
This brief review of literature highlights a number of different considerations associated with Zipf's law and the rank-size distribution. These considerations follow an empirical point of view as well as a theoretical standpoint. One of the reoccurring issues while examining Zipf's law has been its sensitivity to the changes in the definition of the urban unit commonly referred to as Modifiable Aerial Unit Problem [61]. In other words, depending on the definition of the urban units based on which Zipf's law is being examined, there may be a change in the results. Historically, for Zipf's law, the size of the city has been estimated using the population within the administrative boundaries of the city [44]. Although the use of administrative boundaries as urban unit areas has been viewed as misleading and contributing to inconsistencies, particularly due to the political and social influence in the assignment of these boundaries. Rather than the boundaries acting as a demarcation for built-up urban areas and human activities [62], they serve as a good starting point. Other studies, such as by Rauch [63] and Reggiani & Nijkamp [64], have indicated that an urban agglomeration, which consists of the central city along with its hinterland, could provide the best representation for urban rank-size distribution, rather than only including the central cities. Previously, there has been no one uniform standardized definition of cities and urban areas [65,66]. The definition followed arbitrary rules, depending on the country [67]. However, more recently, advancements in geospatial technology and big data have made it possible to define cities based on the concept of natural cities through a variety of processes, such as the bottoms up approach, based on building footprints [67] and nighttime light data [20,62].
The difference in the identification of urban unit areas considered in the analysis has a significant impact on the estimates of Zipf's exponent derived from it. The impact of change in the definition of urban areas on the conformity to Zipf's law has been highlighted in a number of different studies, including Soo [52] and Nitsch [58], where the use of agglomeration data compared to population census data yielded two contrasting results. For the US, studies by Krugman [42] and Gabaix [25], using the US metropolitan level data, indicated its adherence to the Pareto distribution, but using census-defined places from the 2000 US census, Eeckhout [41] showed that data at this scale deviated from the Pareto distribution. Additionally, Berry & Okulicz-Kozaryn [40] indicated that the specific clustering of urban areas into a megalopolitan urban-economic region allowed the US city size to strictly follow the Zipfian distribution. More recently, an increasing number of studies have applied the concept of natural cities and are using remotely sensed images to investigate the conformity of city sizes to Zipf's law. The application of remotely sensed images to Zipf's law has included their use for urban area detection by land cover, land use classification [68,69], nighttime light [57,70] as well as urban surface heat [71].
The sizes of the samples used to analyze Zipf's law were also shown to have a significant influence on the results [72]. Using population data from French cities, Guérin-Pace [72] showed an overall decrease in inequality in the distribution of population when samples of cities were truncated based on a threshold population size that resulted in the selection of cities at the top of the hierarchy. Similarly, using population data for Chinese cities, Gangopadhyay and Basu [54] and Huang et al. [57] found that a higher amount of truncation of tail-end data resulted in an increase in evenness in the rank-size distribution of cities. Peng [55], using a rolling sample regression for Chinese cities, also concluded that an increase in the amount of truncation in the data led to a decrease in the Pareto coefficient, thus indicating an increased evenness in the rank-size distribution. Eeckhout [41] further stated that for the power law to apply to cities, it is essential to introduce a truncation point into the analysis. Arshad et al. [50], in the review paper, further noted that all the studies that have supported Zipf's law have worked with a truncated sample and data from the upper tail of the distribution, but also acknowledged that the choice of truncation points by the researchers introduces arbitrariness into the sample.
Skeptics such as Ausloos and Cerqueti [73] have criticized Zipf's law as being inadequate to describe city rank-size and have suggested an alternative distribution based on arguments from theoretical physics. Benguigui and Blumenfeld-Lieberthal [74] have also suggested that Zipf's law should not be used as a universal law and have called for an abandonment of the paradigm of Zipf's law. Furthermore, methodological changes suggesting a reevaluation of Zipf's law, away from Equation (3) and OLS regression, and the use of Rank-1 2 [36] and the implementation of recursive regression [75] have been suggested. Additionally, there has also been much discussion related to the use of alternate distribution, such as a log-normal distribution [41,75,76] and the double Pareto log-normal distribution [43,77], better fitting the rank-size rule as compared to the Pareto distribution. These deviations away from the rank-size, seen as a result of the definition of urban areas and the truncation of data, along with the simplicity in the formula used to derive the distribution, have been some of the major sources of skepticism and controversy surrounding Zipf's law [72].
In reality, the Pareto coefficient or Zipf's exponent value corresponding to 1 is rarely achieved [72]. However, rather than examining whether Zipf's law is accepted or rejected for a particular urban system, questions regarding how well or how poorly the theory fits the urban system should be evaluated [36]. Guérin-Pace [72] calls for change in the way in which city rank-size is interpreted and emphasizes the need to consider the deviations away from the Pareto distribution as a way to examine the "dynamic process of evolution of a city size and of the urban system" (p. 561) over time. Similarly, with respect to urbanization, Wu et al. [62] describe Zipf's law as proving a universal degree of measure describing the urbanity of a system. Furthermore, Guérin-Pace [72] underlines the need to investigate the explanations for the observed deviation away from the Pareto distribution. These deviations in the distribution away from Zipf's law can be attributed to a myriad of underlying distortions in urban systems, such as institutional factors, economic factors, the allocation of resources or even mere accidental factors [50].

Identification of Urban Areas for Bangladesh
The definition of urban units has a big influence on the overall outcome when examining Zipf's coefficient [21]. There have been many calls by scholars for a consistent definition, allowing for the delineation of urban areas and urban population [65,66]. Urban areas have been identified based on a variety of definitions, which have included those based on administrative regions, population size, population density, population characteristics or a combination of these criteria [65,66]. An example of those based on economic activities is the requirement that the majority of the population be employed in non-agricultural activities [78]. Other definitions of urban areas have been based on administrative functions, facilities and infrastructure such as schools and municipal buildings within a geographic area [78,79]. More recently, urban areas have been defined based on the concept of natural cities, derived through the identification of areas of human activities and urban agglomeration [20,80,81].
For this study, we derived the definition of urban units using a two-step process by combing traditional administrative boundary data with the concept of natural cities [20,80]. Traditionally, the urban administrative hierarchy structure in Bangladesh has been divided into four categories-namely megacities, Statistical Metropolitan Areas, Pourashavas or municipalities and other urban areas (OUA) based on the population size and the administration/governance hierarchy [82]. Currently, Dhaka is the only megacity in Bangladesh that includes an agglomeration of Dhaka City Corporation, along with several municipalities and a number of villages or union parishads (UP) [82,83]. The Statistical Metropolitan Areas, also known as City Corporations, have been defined as larger areas showing urban characteristics, where the definition of urban, as given by the Bangladesh Bureau of Census, includes a central place with infrastructure such as paved roads, electricity, water supply and sanitation, and where a majority of the population is working in the non-agricultural sector [84]. Smaller urban areas that also show urban characteristics have been defined as Pourashavas or Municipalities by the Ministry of Local Government, Rural Development and Cooperatives [82,85]. Finally, Other Urban Areas (OUA) are usually Upazila headquarters that have big market areas but have not been officially declared as Pourashava [82].
The administrative boundaries for city corporations and municipalities, as defined by the Bangladesh Bureau of Statistics for 1991, were used as the baseline for the spatial extent to identify built-up urban areas, where human activities and urban agglomerations occur, using Google Earth Engine (GEE) based land cover classification. The population data for the municipalities and city corporations was derived from the 10-year census from 1991, 2001 and 2011 for these revised urban areas. Although there was not a perfect overlap between the number of municipalities listed for population data and those used to examine the built-up urban areas, due to the revision of municipalities in the later years, a comparison between the two can provide a valuable understanding of the distribution of population and urban areas in Bangladesh.

Google Earth Engine
As a result of the limited availability of land cover maps for the whole of Bangladesh for multiple time periods, we employed GEE to produce land cover data. This GEE-based land cover map was specifically targeted to extract built-up urban areas in Bangladesh. GEE has gained much traction with regards to its application in the field of remote sensing. It stands out as a platform that provides simultaneous and instant access to a wide variety of public and private data as well as a diverse set of scalable computation and remote sensing applications, all into one online framework, thus creating a powerful tool for researchers, practitioners as well as other stakeholders interested in using remote Sustainability 2020, 12, 4643 9 of 26 sensing and geographic information methodologies [86]. More importantly, the integrated GEE system, making use of the cloud-based infrastructure, considerably reduces the need for local machines with high processing power and decreases the need for bandwidth intensive download of large remotely sensed data files [87]. These features within the GEE platform have provided the access and processing power needed to work with freely available, large and relatively high spatial and temporal resolution data such as Landsat [88]. Additionally, the ease of use of the platform for users with limited knowledge of programming languages and remote sensing techniques [89], as well as capabilities within the system for advanced programmers to incorporate sophisticated methodologies, make it a flexible and a powerful tool [90]. As Azzari & Lobell [86] state, the GEE platform, "by simplifying access and processing of large amount of satellite data, are changing the paradigm in land cover monitoring from a static, product-based approach into a more dynamic and application-specific one without any loss of accuracy" (p. 65).
The GEE platform has been previously implemented to examine a wide variety of topics, ranging from inspection of wetland degradation in Costa Rica [89] to the identification of mining areas in Brazil [91], snow mapping and hydrologic modeling [92] and automated mapping of crop land [93], as well as disaster management and earth sciences related research [94], to name a few. However, a significant portion of the literature associated with GEE still involves its application in land cover and land use mapping.

Data for Classification
Data used for land cover classification using the GEE platform comprised the Landsat 5 Top of Atmosphere (TOA) for 1990, 2000 and 2010 and the Landsat 8 TOA for 2019. Spectral bands corresponding to Blue, Green, Red, Near Infrared, Short Infrared 1, Short Infrared 2 and Thermal Infrared available within the Landsat images were used for classification. A four-image collection corresponding to the four years were created using 14 individual Landsat image tiles mosaicked together (except for 2019, where only 13 tiles were used). The Landsat TOA tiles used in the mosaic were limited to images collected during the months of January to May and with a cloud cover of less than 5 percent. Although not ideal, in cases when the cover was less than 5 percent, data from the following year for those months was used.

Training and Validation Points
The classification scheme used in this study was loosely based on the Fine Resolution Observation and Monitoring of Global Land Cover (FROM-GLC) level I scheme, one of the most commonly used global land cover products [95], that was originally introduced by Gong et al. [96]. For this study, the FROM-GLC level I scheme containing seven categories (cropland, forest, grassland, shrubland, water, bareland and snow/ice) was simplified to 5 categories. The simplification process involved combining forest, shrubland and grassland into a single category classified as vegetation, renaming bareland to barren/sand, renaming cropland to agriculture, eliminating the snow/ice class and adding the urban area category to the classification.
The use of higher-resolution images for the creation of training and validation points has been a well-used methodology. This methodology has been previously implemented in studies such as Hansen et al. [97] and Midekisa et al. [88], where Landsat images were used to create training and validation data points. Traganos et al. [98] used ArcGIS's based World Image base map to manually digitize the training data and [91] used high resolution Rapid-eye images for validation purposes. The training and validation data points used in this study for the purpose of classification were derived from a visual inspection of freely available high-resolution images ranging from 30 cm to 15 cm at its highest [99] to medium resolution images from Landsat, using Google Earth (GE). Through the inspection of high-resolution images, the areas identified as built-up urban areas included concrete buildings, metal and clay rooftops, asphalt roads as well as other built infrastructures; the barren/sand category included areas along the beach and riversides that had sandy deposits and rockfaces on the side of hills; and the agriculture category included areas that were being actively farmed at the time when the image was acquired as well as areas that were not being actively farmed at the time of the image acquisition. Finally, the snow/ice class was eliminated because it was largely irrelevant for Bangladesh. For the validation of the classified land cover, considering the entirety of Bangladesh, 1150 to 1200 sample validation points were manually digitized such that they were evenly distributed amongst the five land cover classes (see Figure 2). When only the cities and municipalities of Bangladesh were considered, 500 random points were generated within these cities and the municipalities using ArcGIS (Figure 3). These validation points were manually referenced, using high-resolution images within the GE, when available, or using Landsat images when high-resolution images were not available. Additionally, changes to the location and the number of the validation  (Figure 1). Additionally, changes to the shape, location and number of the training polygons were made for each of the four years in order to reflect the changes in land cover conditions over the years.
Sustainability 2020, 12, x FOR PEER REVIEW 12 of 28 points were made for each of the four years in order to reflect the changes in land cover conditions at the given point over the years.  For the validation of the classified land cover, considering the entirety of Bangladesh, 1150 to 1200 sample validation points were manually digitized such that they were evenly distributed amongst the five land cover classes (Figure 2). When only the cities and municipalities of Bangladesh were considered, 500 random points were generated within these cities and the municipalities using ArcGIS (Figure 3). These validation points were manually referenced, using high-resolution images within the GE, when available, or using Landsat images when high-resolution images were not available.
Additionally, changes to the location and the number of the validation points were made for each of the four years in order to reflect the changes in land cover conditions at the given point over the years.

Classification Using Classification and Regression Trees (CART)
In this study, the classification of land cover was preformed using Classification and Regression Trees (CART), a machine-learning-based supervised classifier [100], available within the GEE platform. According to Steinberg & Colla [101], "The CART decision tree is a binary recursive partitioning procedure capable of processing continuous and nominal attributes as targets and predictors" (p. 181). In other words, CART is a machine-learning-based method of creating prediction

Classification Using Classification and Regression Trees (CART)
In this study, the classification of land cover was preformed using Classification and Regression Trees (CART), a machine-learning-based supervised classifier [100], available within the GEE platform. According to Steinberg & Colla [101], "The CART decision tree is a binary recursive partitioning procedure capable of processing continuous and nominal attributes as targets and predictors" (p. 181). In other words, CART is a machine-learning-based method of creating prediction models by recursively partitioning the data and fitting a simple prediction model in each of these partitions [102]. De'ath & Fabricius [103] describe the methodology as being flexible, robust and easily interpretable, with the ability to handle missing values and the capability of handling a variety of data types such as numeric, categorical, ordered as well as survival data. Additionally, the well-known problem of imbalanced classes in machine learning has been explicitly addressed by the authors of CART. Regardless of the extent of imbalance in the training data, CART automatically adjusts to the imbalance without the need for external actions such as sampling or weighing [101]. The classification methodology associated with CART has been best described as a tree-like structure where data from the root node splits into two children, which further split into its grandchildren [101]. The best split of the node is based on the split that minimizes the sum of squares [104]. Once the maximum number of splits is reached, the tree is pruned based on the split that contributes the least to the overall performance of the tree, resulting in the smallest misclassification error [102]. This mechanism produces sequences of nested trees, each of which is a candidate for the optimum tree [101]. Finally, based on cross validation, the predictive performance of each of the pruned trees is evaluated and the final right-sized tree is identified [101,103].
Amongst other supervised classifiers available within GEE, based on the past studies, CART was determined to have the best performance in comparison. A previous study by Farda [105] revealed that when using supervised machine learning for classification in GEE, CART was shown to outperform the Support Vector Machines (SVM) classifier. Comparing the classification accuracy of CART, SVM and Random Forest (RF), Goldblatt et al. [106] indicated that among the three, SVM had the lowest classification accuracy, followed by CART and RF. These results were further verified by Shelestov et al. [107] in a review of the performance of classification algorithms available within GEE, where the authors indicated CART as having the best performance as compared to other algorithms.
The land cover classification results produced using CART were exported to the Google drive in the Tiff format for the 9 individual sectors. These exported sectors were combined as a single map using ArcMap. In order to get a smoother representation of the classified areas, the majority filter tool also available within the ArcMap was run. Finally, areas classified as built-up urban areas within the cities and municipalities were extracted from the results. The steps involved in the land cover classification process are highlighted in the flowchart (Figure 4).

Validation and Classification Results
The assessment of accuracy is an important component within the context of remote sensing, as it associates the remotely sensed and classified images with the observations in reality [98]. The accuracy of the classifier was examined by using the widely implemented confusion matrix or error matrix methodology [108]. The overall accuracy of the land cover classification for the whole of Bangladesh was calculated between 70.45 and 72.68 ( Table 2). The Kappa coefficient, which indicates the measure of agreement between the actual land cover and the predicted land cover [109], was calculated between 57.74 and 62.68 for their respective time frames (Table 2). Although the Kappa value is relatively low, based on the commonly cited scale for interpretation of Kappa values by Landis & Koch [110], the agreement between the predicted and the actual land cover was defined as being between moderate and substantial levels. Furthermore, when only the cities and municipalities were taken into consideration, the accuracy of the classification showed significant improvement. The overall accuracy of land cover classification within the cities and municipalities was calculated between 84.2 and 87.4 percent, and the kappa coefficient was between 73.06 and 77.81 for the 4 years, respectively ( Table 2). Based on Landis & Koch [110], this level of agreement between the predicted and actual land cover was shown to be substantial. Sustainability 2020, 12, x FOR PEER REVIEW 15 of 28

Validation and Classification Results
The assessment of accuracy is an important component within the context of remote sensing, as it associates the remotely sensed and classified images with the observations in reality [98]. The accuracy of the classifier was examined by using the widely implemented confusion matrix or error matrix methodology [108]. The overall accuracy of the land cover classification for the whole of Bangladesh was calculated between 70.45 and 72.68 (See Table 2). The Kappa coefficient, which indicates the measure of agreement between the actual land cover and the predicted land cover [109], was calculated between 57.74 and 62.68 for their respective time frames (See Table 2). Although the Kappa value is relatively low, based on the commonly cited scale for interpretation of Kappa values by Landis & Koch [110], the agreement between the predicted and the actual land cover was defined as being between moderate and substantial levels. Furthermore, when only the cities and municipalities were taken into consideration, the accuracy of the classification showed significant improvement. The overall accuracy of land cover classification within the cities and municipalities was calculated between 84.2 and 87.4 percent, and the kappa coefficient was between 73.06 and 77.81 for the 4 years, respectively (See Table 2). Based on Landis & Koch [110], this level of agreement between the predicted and actual land cover was shown to be substantial.   As the study is largely concerned with built-up urban areas, the producer accuracy and the user accuracy for the built-up urban areas are also computed. Producer accuracy provides information regarding how well a particular region is being classified and is derived by taking the ratio of validation samples that were correctly classified and the actual number of validation samples in the class [98,111]. User accuracy provides information regarding the reliability of the map and is derived by dividing the total number of validation samples in a given class by the total number of validation samples that are claimed to be in the class [98,111]. The producer accuracy of the built-up urban areas for the entirety of Bangladesh was calculated within the range of 40.4 and 62.93 percent, and the user accuracy was calculated in the range of 57.14 and 91.58 percent (Table 2) between 1990 and 2019. When considering only the cities and municipalities, the producer accuracy ranged from 73.53 to 85.96 percent, and the user accuracy ranged from 79.03 to 93.62 percent over the given years.
One of the major challenges faced during the classification process was due to the typical structure of built-up urban areas in the smaller municipalities in Bangladesh. In these locations, the built-up areas were generally comingled with, and surrounded by, vegetation such as trees and shrubs ( Figure 5). Considering the Landsat resolution of 30 m, this proximity of vegetation alongside the built-up areas contributed to a higher degree of error in the classification, particularly for urban areas.
by dividing the total number of validation samples in a given class by the total number of validation samples that are claimed to be in the class [98,111]. The producer accuracy of the built-up urban areas for the entirety of Bangladesh was calculated within the range of 40.4 and 62.93 percent, and the user accuracy was calculated in the range of 57.14 and 91.58 percent (Table 2) between 1990 and 2019. When considering only the cities and municipalities, the producer accuracy ranged from 73.53 to 85.96 percent, and the user accuracy ranged from 79.03 to 93.62 percent over the given years.
One of the major challenges faced during the classification process was due to the typical structure of built-up urban areas in the smaller municipalities in Bangladesh. In these locations, the built-up areas were generally comingled with, and surrounded by, vegetation such as trees and shrubs ( Figure 5). Considering the Landsat resolution of 30 m, this proximity of vegetation alongside the built-up areas contributed to a higher degree of error in the classification, particularly for urban areas. The land cover classification results for the cities and municipalities in Bangladesh revealed agriculture as being the most dominant land cover class, followed by vegetation (Figures 6 a-d). Agriculture and vegetation accounted for just over 90 percent of the land cover within the The land cover classification results for the cities and municipalities in Bangladesh revealed agriculture as being the most dominant land cover class, followed by vegetation (Figure 6a-d). Agriculture and vegetation accounted for just over 90 percent of the land cover within the municipalities, changing from 93.92 percent of the total land cover in 1990 to 90.11 percent in 2019. In the case of built-up urban areas, although it accounted for a relatively small percentage of the total land cover occupied within the cities and municipalities of Bangladesh, the actual amount of built-up urban areas within the cities and municipalities more than doubled in the last 30 years, increasing from 152.68 sq. km. to 333.64 sq. km. in 2019, which was an increase from 2.89 percent in 2010 to 6.31 percent in 2019 (Table 3). Conversely, the largest overall decrease in land cover was seen for agricultural land cover, with a decrease of just under 230 sq. km. during the same period (Table 3).

Conformity to Zipf's Law
The conformity to Zipf's law for the distribution of population and built-up urban areas within cities and municipalities in Bangladesh was tested using the OLS regression (see Eq. 4), and Wald's T-test was used to examine if the value of Zipf's coefficient (q) significantly differed from 1. To When specifically examining the location of the built-up urban areas, a majority of these areas were located within the two major cities, Dhaka and Chittagong. These built-up urban areas in Dhaka and Chittagong accounted for 50 to nearly 60 percent of the total built-up urban areas within the cities and municipalities in Bangladesh over the last 30 years (Table 4).

Conformity to Zipf's Law
The conformity to Zipf's law for the distribution of population and built-up urban areas within cities and municipalities in Bangladesh was tested using the OLS regression (Equation (4)), and Wald's T-test was used to examine if the value of Zipf's coefficient (q) significantly differed from 1. To examine the distribution of population and built-up urban areas in top end cities and municipalities, particularly in a heavy-tailed distribution, a grouping of these top end cities and municipalities was derived using head/tail breaks. Head/tail breaks provide a naturally determined method of categorizing the data based on inherent hierarchical levels within the data [112] and hence reduced the arbitrariness when truncating the samples. Based on the head/tail breaks depending on the year, for population, the upper 8.52 percent to 13.18 percent, and for built-up urban areas, the upper 10.36 percent to 13.83 percent of cities and municipalities were considered as top tier (Table 5). When examining the conformity of population among cities and municipalities to Zipf's law, results from the OLS regression for 1991 indicated the Zipf's coefficient to be 1.01, and the T-test showed that the value was not significantly different from 1. This signaled that the distribution of population within the cities and municipalities for 1990 conformed to Zipf's law. Over the next 20 years, the Zipf's coefficient in reference to population decreased to 0.87 and to 0.85 for 2001 and 2011, respectively, both of which were significantly different than 1 at 1 percent significance level (Table 6). From this result we can infer that there was an overall increase in evenness in the distribution of population amongst the cities and municipalities in Bangladesh, as seen in the more horizontal profile of the plot showing the relation between population and rank (Figure 7). These results can be attributed to an increase in population within the mid-sized and smaller sized cities as compared to that of larger cities from 2001 to 2011, hence leading to an evening out of the overall distribution of population amongst the cities and municipalities in Bangladesh.
When limiting the distribution to the top end cities and municipalities, the Zipf's coefficient of population was calculated to be 1.18, 1.32 and 1.27 for 1991, 2001 and 2011, respectively (Table 6). For all 3 years, the coefficients were seen to be statistically different than 1 at 1 percent significance level, and thus they did not conform to Zipf's law. From this we concluded that when only the top end cities and municipalities were considered, the population was still seen to be unevenly distributed, with most of the population being concentrated in the top portion of the top end cities and municipalities.  Likewise, when examining the conformity of built-up urban areas within the cities and municipalities to Zipf's law, the results from the OLS regression indicated the value of Zipf's coefficient for all four years to be significantly different from 1. The Zipf's coefficient for these top end built-up urban areas ranged from 1.68 in 1990 to 2.04 in 2019 (See Table 7). The increasing values of Zipf's coefficient indicated an increased unevenness in the distribution of built-up urban areas within the cities and municipalities over the years (See Figure 8). Here, the bigger cities were seen to host a larger share of the built-up urban areas as compared to the mid and smaller sized cities. These results can be attributed to a continued increase in urban buildup and the extension of the urban footprint occurring in the larger cities as compared to mid and smaller sized cities over the last three decades.
When the truncated sample of top end cities and municipalities was considered, the results also indicated an unevenness in the distribution of built-up urban areas in these cities. The Zipf's coefficient values in the truncated sample were, however, smaller than their corresponding full sample values and ranging from 1.21 in 1990 to 1.11 in 2019 (Table 7). Not all of the Zipf's coefficient values were statistically significant at 1 percent significance level, but they were statistically significant at 5 percent significance level. This indicates that while the distribution of built-up urban areas in the cities tended towards conforming to Zipf's law-a major portion of these built-up urban areas is still concentrated in the upper portion of the top tier cities.  Likewise, when examining the conformity of built-up urban areas within the cities and municipalities to Zipf's law, the results from the OLS regression indicated the value of Zipf's coefficient for all four years to be significantly different from 1. The Zipf's coefficient for these top end built-up urban areas ranged from 1.68 in 1990 to 2.04 in 2019 ( Table 7). The increasing values of Zipf's coefficient indicated an increased unevenness in the distribution of built-up urban areas within the cities and municipalities over the years (Figure 8). Here, the bigger cities were seen to host a larger share of the built-up urban areas as compared to the mid and smaller sized cities. These results can be attributed to a continued increase in urban buildup and the extension of the urban footprint occurring in the larger cities as compared to mid and smaller sized cities over the last three decades.  These results examining the conformity to Zipf's law highlight the inequalities and skewness of urban development while highlighting the increase in population and a subsequent expansion of the urban footprint within the cities and municipalities of Bangladesh. As indicated by the value of the Zipf's coefficient for population, while there has been a greater increase in the concentration of population in the small and mid-sized cities between 1991 and 2011 as compared to larger cities, in contrast, the growth of built-up urban areas has largely been concentrated in the biggest cities in Bangladesh, as indicated by the value of the Zipf's coefficient for built-up urban areas based on all the cities and municipalities and the truncated sample of top end cities and municipalities. The most recent data, from 2019, further indicates that this trend of greater expansion of the urban footprint in the bigger cities has continued during the recent period between 2010 and 2019. When examining the truncated sample of the population within the top end cities and municipalities for the same period of time, we see that the increase in the population of the top end cities has been concentrated in the biggest cities. This is further confirmed by inspecting the density of these cities and municipalities. Based on the 100 × 100 meter population density data from worldpop.org, as published by the University of Southampton [113], we find that while the footprint of the urban areas has increased, the population densities within these built-up urban areas have also increased at a similarly fast pace. For the top two cities in Bangladesh, Dhaka and Chittagong, between 2000 and 2019, the population density per 100 sq. meters within the built-up urban areas increased from 291.34 and 158.91 to 565.94 and 242.28, respectively (Table 8). This indicates that even though there has been an overall expansion in the urban footprint of the cities and municipalities in the top end cities, the increase in population in these cities has resulted in a subsequent increase in population density in these top end cities.  When the truncated sample of top end cities and municipalities was considered, the results also indicated an unevenness in the distribution of built-up urban areas in these cities. The Zipf's coefficient values in the truncated sample were, however, smaller than their corresponding full sample values and ranging from 1.21 in 1990 to 1.11 in 2019 (Table 7). Not all of the Zipf's coefficient values were statistically significant at 1 percent significance level, but they were statistically significant at 5 percent significance level. This indicates that while the distribution of built-up urban areas in the cities tended towards conforming to Zipf's law-a major portion of these built-up urban areas is still concentrated in the upper portion of the top tier cities.
These results examining the conformity to Zipf's law highlight the inequalities and skewness of urban development while highlighting the increase in population and a subsequent expansion of the urban footprint within the cities and municipalities of Bangladesh. As indicated by the value of the Zipf's coefficient for population, while there has been a greater increase in the concentration of population in the small and mid-sized cities between 1991 and 2011 as compared to larger cities, in contrast, the growth of built-up urban areas has largely been concentrated in the biggest cities in Bangladesh, as indicated by the value of the Zipf's coefficient for built-up urban areas based on all the cities and municipalities and the truncated sample of top end cities and municipalities. The most recent data, from 2019, further indicates that this trend of greater expansion of the urban footprint in the bigger cities has continued during the recent period between 2010 and 2019. When examining the truncated sample of the population within the top end cities and municipalities for the same period of time, we see that the increase in the population of the top end cities has been concentrated in the biggest cities. This is further confirmed by inspecting the density of these cities and municipalities. Based on the 100 × 100 meter population density data from worldpop.org, as published by the University of Southampton [113], we find that while the footprint of the urban areas has increased, the population densities within these built-up urban areas have also increased at a similarly fast pace. For the top two cities in Bangladesh, Dhaka and Chittagong, between 2000 and 2019, the population density per 100 sq. meters within the built-up urban areas increased from 291.34 and 158.91 to 565.94 and 242.28, respectively (Table 8). This indicates that even though there has been an overall expansion in the urban footprint of the cities and municipalities in the top end cities, the increase in population in these cities has resulted in a subsequent increase in population density in these top end cities.

Discussion
The distribution of population and built-up urban areas conforming to Zipf's law provides much information on the expected size of these urban areas with regards to their respective ranks. However, based on previous research, we see that a perfect conformity of city rank-size to Zipf's law does not always occur [19,36,41,72,114]. These deviations away from Zipf's law offer distinct insights into the specifics of the distribution of population or urban areas within these cities. In particular, when looking at Bangladesh, the departure from Zipf's law for the rank-size of built-up urban areas and the population within the cities and municipalities sheds light on the changes in the distribution of the urban footprint of these areas relative to the change in distribution of population within these urban areas over the years, since 1991. Furthermore, it draws attention to the rapid pace of development of built-up urban areas and the expansion of the urban footprint occurring in the largest cities, such as Dhaka and Chittagong, compared to the relatively low levels of increase in built-up urban areas in small and mid-sized cities and municipalities, where most of the population growth has been seen.
These spatial inequalities in urban development have led to the unsustainable expansion of urban areas within Bangladesh and have also contributed towards inducing social conflicts within the country [84]. The literature suggests that the strategic distribution of public expenditure can significantly aid in decreasing the level of development inequalities between the cities [115], but in Bangladesh these public investments, along with the private capital investment and administrative decisions, disproportionately favor a select few major cities [1,84,116]. The majority of programs, such as those aimed at building roads and bridges, along with ventures in technological innovations, have been targeted for areas that are already urbanized and well-off. Meanwhile, there is a lack of adequate investment flowing into mid and smaller sized cities that are experiencing an increase in population; a decision that has largely been driven by political influences [116]. These disparities in investment tilted towards the bigger cities have led to a further rise of these cities as the primary economic drivers of the country and have contributed to the pull factor of these cities, resulting in an increase in population as well as in the urban footprint of these areas. Conversely, economies in smaller and mid-sized cities have mostly remained agriculture-based, with diminished levels of revenue potential and an overall lower income per capita [84]. This has contributed to the push factor of these regions. While there has been an expansion of urban built-up areas in the bigger cities of Bangladesh, these push and pull factors have led to the migration of poorer and vulnerable population from rural areas into these cities, which has contributed to a similarly large increase in the population density in these cities [117]. A considerable portion of the increase in the urban footprint within these cities has been unplanned and unguided, with disregard to environmental and human considerations, leading to major concerns regarding the sustainability of these urban areas. The increased pressure on the local capabilities and natural resources of these cities due to the increases in population and urban footprint of the cities have not been addressed with adequate infrastructural upgrades and facilities to support the expansion, leading to issues with pollution, traffic, crime and environmental degradation [118]. Despite higher influxes of investments in these bigger cities, basic amenities such as water, roads, electricity and sewage and the overall planning necessary to support the increasing population and urban footprint of these cities have been insufficient. Furthermore, the lack of infrastructure and services is particularly prominent in slum settlements within these bigger cities [119]. This uncontrolled and unrestricted expansion can be attributed to the lack of-or, at times, selective-enforcement of policies, inefficient allocation and distribution of resources, along with an inadequate urban planning strategy, which has largely stemmed from a lack of good governance [1]. There is a distinct need for equitable provisions for physical infrastructure [120] and a targeted pro-poor strategy within the context of urbanization and urban development [121], to promote a current and future sustainable urban development that is in balance with the environment [120].

Conclusion
The issues of urban development and sustainable urban growth are extremely complex, particularly for a developing country such as Bangladesh, with significant constraints on financial resources. With regards to the increasing urban footprint of the larger cities, as indicated by the nonconformity to the Zipf's coefficient for built-up urban areas, policies that are aimed at curtailing the uncontrolled haphazard urban expansion, particularly in the largest cities, would need to be implemented. Moreover, for a successful realization of these urban growth restriction strategies, supporting policies targeting smart development, the equitable distribution of resources and a pro-poor urban development would also need to be prioritized and implemented simultaneously. These policies would need to focus on the planned and guided development of built-up urban areas within the city limits, in a manner that minimizes the impact on the natural environment and surrounding area while addressing the socio-economic issues of informal settlements, slums and the reduction of non-income poverty, such as the availability of urban services [121]. With regards to the increasing concentration of population in the mid and smaller sized cities, as indicated by the nonconformity to Zipf's coefficient for population, policies primarily targeted at the reduction of economic stagnation, increased opportunities and an overall economic development of the region, with an emphasis on guided urban growth that ensures environmental consciousness in these cities, would need to be implemented. In both cases, the policy tools put into place would have to be multifaceted and aimed at environmental conservation and natural resources management with sustainable urban expansion, while keeping in mind the underlying factors of economic growth, poverty alleviation strategies and an overall improvement in living standards that would impact the push and pull factors of the region. These policies would need to be employed at different administrative levels, from the central to the local government, with a common policy goal. Most importantly, the success of these agendas, targeted at sustainable development and growth, as well as improved livability, will primarily depend on the good governance, efficient implementation and enforcement of these policies, without its misuse [115,122].
Traditionally, developmental planning and policy making in Bangladesh has been based on a sectoral format, as opposed to a spatial approach [84]. The integration of a spatial approach within the overall planning and policy framework would entail the examination of anthropogenic characteristics and natural features, along with socioeconomic attributes such as economic conditions and disparity in access and living conditions [123]. Using a spatial framework would allow for identifying geographically explicit development issues while defining clear outcomes through the incorporation of different economic, environmental, cultural and social agendas specific to the particular locality, and thus allowing for the implementation of targeted policies for the region [123,124]. Additionally, the use of spatial planning aimed at territorial development could act as the strong base needed to improve intra-governmental partnerships and coordination, and to incorporate other stakeholders under its umbrella. However, the implementation of a spatially explicit plan has generally been seen as one of the weakest links within the overall planning framework, and has been frequently carried out without sufficient theoretical exploration or methodological rigor [125]. The incorporation of robust theoretical underpinning, such as Zipf's law, to substantiate empirical observations would augment the spatial component of the overall urban developmental planning framework of Bangladesh, and would no doubt pave the way to usher in targeted policies such as the Urban Growth Boundaries, aimed at restricting and guiding the direction of urban expansion for cities and municipalities in Bangladesh.