Land-cover change analysis in 50 global cities by using a combination of Landsat data and analysis of grid cells

Global urban expansion has created incentives to convert green spaces to urban/built-up area. Therefore, understanding the distribution and dynamics of the land-cover changes in cities is essential for better understanding of the cities’ fundamental characteristics and processes, and of the impact of changing land-cover on potential carbon storage. We present a grid square approach using multi-temporal Landsat data from around 1985–2010 to monitor the spatio-temporal land-cover dynamics of 50 global cities. The maximum-likelihood classification method is applied to Landsat data to define the cities’ urbanized areas at different points in time. Subsequently, 1 km2 grid squares with unique cell IDs are designed to link among land-cover maps for spatio-temporal land-cover change analysis. Then, we calculate land-cover category proportions for each map in 1 km2 grid cells. Statistical comparison of the land-cover changes in grid square cells shows that urban area expansion in 50 global cities was strongly negatively correlated with forest, cropland and grassland changes. The generated land-cover proportions in 1 km2 grid cells and the spatial relationships between the changes of land-cover classes are critical for understanding past patterns and the consequences of urban development so as to inform future urban planning, risk management and conservation strategies.


Introduction
Urban areas occupy a relatively small fraction of the Earth's land area, but at present more than half of the global population lives in urban areas, and this proportion is expected to increase in the coming decades (http://esa.un.org/unup/). Already the intensive burning of carbon fuels in the world's urban areas accounts for about 70% of global greenhouse gas emissions (Solecki et al 2013). In addition, previous research suggests that a 10% increase in urban land cover in a country is associated with an increase of more than 11% in the country's total CO 2 emissions (Angel et al 2011).
Urban area expansion and associated land-cover change also contribute to the loss of terrestrial carbon stored in vegetation biomass (Seto et al 2012); threaten biodiversity (McKinney 2002); and result in waste, pollution, and environmental degradation (Millington et al 1999). Therefore understanding how urbanization affects vegetation dynamics will provide information useful for determining the amount and timing of ecosystem services availability for urban inhabitants (Jenerette et al 2013).
The monitoring and mapping of urban growth and the development of effective urban planning strategies require knowledge of the spatio-temporal extent and expansion trends of cities. Urban expansion can take place through substantially different forms: the redevelopment of built-up areas at higher densities, infill of remaining open spaces in already built-up areas, or the development of greenbelt land around cities (Angel et al 2011).
Although urban prosperity takes different forms and characterizations, urban growth patterns diverge amongst regions and countries. For example, in the United States people tend to live in low-density single-family homes and commute by car to work. By contrast, in Japan high-rise residential buildings dominate and workers commute by public transportation (e.g., metro or rail line). Thus urban land-cover change may vary from region to region because of the diversity of interurban variability as a result of socioeconomic, cultural, historical and environmental differences among different cities across the world (Small et al 2005). Therefore, it is important to compare cities that have different urban growth patterns for their effect on sustainable urban planning and development.
Various available global data sets are used for measuring, analyzing and, hence, understanding the complex processes of urbanization . Examples are global urban extent maps based on, e.g., NOAA Air Force Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS) sensor night-time lights imagery (Elvidge et al 2007, Doll 2008 or MODIS data . The distribution of cities around the world broadly corresponds to the brightness distribution of DMSP night-time lights (Zhang and Seto 2011, Parés-Ramos et al 2013. However, no single DMSP brightness threshold is valid for extracting the urban extent of all cities because small settlements that are not frequently lit are likely to be excluded (Small et al 2005). Although coarse spatial resolution (from 250 m to 2 km) monitoring provides global and national estimates of urban growth, coarse data may be less reliable for correctly estimating the urban area of cities and often results in either overestimation or underestimation of urban areas (Potere and Schneider 2007). Townshend and Justice (2002) argued that 'a substantial proportion of the variability of land cover change has been shown to occur at resolutions below 250 m', and Giri et al (2013) reported that land parcels managed at a local scale are often smaller than the resolution of coarse spatial resolution satellite data. Thus, coarse geometric resolution is a clear restriction for the tracing of smallscale urban outlines, extents and patterns (Taubenböck et al 2012).
Today, higher resolution sensor systems are available (e.g., IKONOS, QuickBird, WorldView and GeoEye1-2) for monitoring the spatial effects of urbanization that provide spatial information content that is hundreds of times better than coarse spatial resolution data sets. Small (2003) compares 14 cities at a very high geometric level using QuickBird data with a sub-metre geometric resolution and derives parameters such as vegetation fraction ;Berger et al (2013) extract urban land cover information from high spatial resolution multi-spectral and light detection and ranging (LiDAR) data; Taubenböck et al (2012) analyze the spatial effects of urbanization over a span of almost 40 years using Landsat data. However, little attention has been paid to the quantitative analysis of correlations among land-cover category changes.
Monitoring urban land-cover dynamics with 30 m resolution is possible by using Landsat TM/ETM+ satellite data that are freely available from 1982. Landsat satellite images are ideal for identifying and characterizing both natural and anthropogenic changes over large areas of land because of the system's acquisition, processing and distribution strategies (Hansen and Loveland 2012).
Recently, spatio-temporal analyzes of land-cover changes using 1 km 2 grid cells have demonstrated that grid cells provide a new way to obtain spatio-temporal information about areas that are smaller than the municipal scale and uniform in size (Bagan and Yamagata 2012) and to further develop the change dynamics analysis to better characterize the phenomena using limited available data. The relationships among changes in urban land-cover patterns can be better understood if these data are mapped onto a grid composed of square grid cells. Furthermore, it allows the commonly used multivariate linear regression model to be used to analyze urbanization trends.
The primary goal of this research is to improve the understanding of urban growth during the period around 1985-2010 in 50 cities across the globe. A secondary objective is to quantify the spatio-temporal dynamics of urban growth. To characterize the nature of the changes and to compare and contrast trends across the 50 cities, we used the maximum-likelihood classification (MLC) method to produce land-cover maps from Landsat images recorded between 1985 and 2010. These land-cover maps are then combined with 1 km 2 grid cells to analyze the spatio-temporal land-cover changes and investigate their statistical properties. These grid cell data allow us to examine the different urban expansion forms in a global sample of 50 cities, help us to develop sustainable urban neighbourhoods and determine the impact of urbanization on the environment over a period of more than two decades.

Distribution of the 50 cities
We selected 50 global cities from five continents, including ones from developed (Europe, North America and Japan) and developing regions (Asia, Africa and South America). Converting land to urban land use leads to increases of greenhouse gas emissions (Angel et al 2011). For further investigate the relationship between urban land cover changes and greenhouse gas emissions, out of the 50 cities, 44 (not including Amsterdam, Berlin, Bern, Rome, Ottawa and Lisbon) are studied by Kennedy et al (2011) to estimate urban greenhouse gas emissions. Hence, if we have detailed urban growth patterns of these cities, it may allow us to explore relationships between the greenhouse gas emissions and urban growth patterns, which are critical to evaluate land usedriven greenhouse gas emissions from urban regions. This Environ. Res. Lett. 9 (2014) 064015 H Bagan and Y Yamagata was not the focus of our work here. However, we will investigate this point in our future research. Urbanization is characterized by notable gains in urban/ built-up areas at the expense of green or open spaces. In this study, the urban/built-up land-cover class includes all nonvegetative, human-constructed elements such as buildings, asphalt and concrete, including residential, commercial, industrial, office use and transportation space. Other types of urban land use, such as golf courses, urban green parks and nature areas, are not included in the urban/built-up class.
The MODIS Land Cover Type product (MCD12Q1: Land Cover Type Yearly L3 Global 500 m SIN Grid V005) for 2005 was used to delineate the urban boundaries of all 50 cities to further refine the study landscape instead of using the administrative boundary of the city or metropolitan area (figure 1). First, we extracted the urban/built-up land-cover class and then buffered urban/built-up pixels up to a distance of 10 km to define each urban boundary. Table 1 lists the 50 cities and the corresponding Landsat images used for land-cover classification. We acquired Landsat TM and ETM+ images to interpret land-cover change on two separate dates (nominally around 1985 and 2010). Landsat data from these approximate time-points covered nearly all 50 cities with 73 and 72 Landsat scenes per time step (145 Landsat scenes overall). However, owing to cloud coverage over Cape Town (South Africa), Glasgow (United Kingdom), Delhi and Kolkata (India) around these times, we used Landsat data from 1995 and 2000 for Glasgow, 1985and 2000for Cape Town, 1999for Delhi and 2000 for Kolkata (table 1).

Landsat data
All Landsat standard data products are processed using the Level 1 product generation system with the following parameters applied: Cubic Convolution resampling method, 30 m (TM, ETM+) pixel size, universal Transverse Mercator map projection, World Geodetic System 84 datum and Northup image orientation. The standard terrain correction, which provides systematic radiometric and geometric accuracy by incorporating ground control points while employing a Digital Elevation Model for topographic accuracy, is used (http://landsat.usgs.gov).

Land-cover classification
In both developing and developed countries, cities from different regions are characterized by different land-cover types and urban sprawl patterns. They vary in shape and size (e.g., overstretched blocks and low-density development) and are surrounded by a variety of landscapes (e.g., agriculture, forest, grassland and desert) (Schneider and Woodcock 2008). Therefore, a key requirement for this study was to create landcover classes to standardize the classifications and therefore facilitate comparisons. We adopted the six broad land-use categories described by the Intergovernmental Panel on Climate Change (IPCC) (IPCC 2003) for land-cover classes (table 2).
Very high resolution satellite data (from, e.g., QuickBird, IKONOS and WorldView-2), acquired from web-based resources such as Google Earth and visual interpretation of the remote sensing data were used to acquire the ground truth dataset. Land-cover classification accuracy at different times of the year is not constant due to seasonal variations in spectral characteristics of different land cover classes. Hence, knowledge about suitable seasons and timing for land-cover mapping is important for generating accurate land cover maps (Sinha et al 2012). In consideration of the weather conditions, to maximize the vegetation information content for each monitoring date, only images from the green vegetation season and with low cloud cover were used. All analyzes were based on the optical and thermal infrared bands of the TM and ETM+ data, and panchromatic bands were excluded. If the images contained cloud cover and associated shadows, we isolated the clouds and their associated shadow portions by  visual image interpretation, and then substituted the selected portions with clear pixels from Landsat imagery of the same area from the nearest dates. The MLC method was used to produce land-cover maps from the Landsat images. MLC is a classification technique widely used in the remote sensing community. With the MLC method, probability density functions are built for each class based on the training samples' spectral values. During classification all pixels are assigned class membership based on the relative likelihood (probability) of that pixel occurring within each class' probability density function (Lillesand et al 2008). MLC is easy to apply and widely available in image processing and statistical software packages. The newly developed classifier may be improve classification accuracy when compared with the MLC method. However, gathering specific knowledge is an enormous task and tuning parameters is very costly. Therefore, the newly developed classifier are not universally applicable (Xie et al 2008).

Grid cell-based analysis
We created 1 km 2 grid square cells for the 50 cities. The grid square cells enable us to represent the land-cover maps in 1 × 1 km grid square cells. To evaluate the spatio-temporal changes of land-cover categories and compare them, the classified images were intersected with empty 1 km 2 grid square cells and the percentages of the six land-cover types within each cell were computed. The results were stored in a new attribute table. To calculate the percentage of a landcover type within a cell, we divided the sum of the landcover-type area by the area of the cell. For a pixel located on a grid line (or spanning multiple grid lines), we assigned the pixel to the cell that covered the largest part of the pixel. We then calculated the correlation coefficients of the land-cover categories based on the 1 km 2 grid square cells to investigate the relationships among changes in the land-cover categories. Thus, the grid cells enable us to aggregate the categories for each map and to calculate their proportions. Furthermore, they enable us to evaluate the spatio-temporal changes in land-cover categories to allow a much easier statistical comparison of the land-cover changes.
Classification by the MLC method was performed for the 50 cities. Overall accuracy for the resulting land-cover maps for each year ranged between 77 and 95%. The lowest classification accuracy was for New York City in 2009 (77.05%) and the highest was for Seoul in 2006 (94.92%). However, Cropland classes tended to be more strongly affected by changes in observation time than other classes. Therefore, the Cropland class was the most difficult among all categories to classify and was confused mainly with Grassland. Using these maps, we assessed the land-cover changes in the cities caused by urban expansion over the study period, and then examined the relationship between the expansion of urban area and other changes in the land-cover categories.  Land-cover class Class description 1. Forest All land with woody vegetation consistent with thresholds used to define forest land in the national greenhouse gas inventory.

Cropland
Arable and tillage land, and agro-forestry systems where vegetation falls below the thresholds used for the forest land category, consistent with the selection of national definitions.

Grassland
Rangelands and pasture land that is not considered as cropland.

Wetlands
Land that is covered or saturated by water for all or part of the year (e.g., peatland) and that does not fall into the forest land, cropland, grassland or settlements categories.

Settlements
Developed land, including transportation infrastructure and human settlements of any size, unless they are already included under other five categories. 6. Other land Bare soil, rock, ice, and all unmanaged land areas that do not fall into any of the other five categories.

Urban land-cover change analysis
There are several major trends evident in the changes of land cover over the approximate 1985-2010 period. The suburban areas surrounding these core cities experienced enormous increases in Settlements area, and there was a marked decrease in green space. In Beijing (figure 2) Figure 3 shows the rates of increase in Settlements area during the study period. More than 32 of the cities grew rapidly by more than 30% between around 1985 and 2010. The fastest growing city was Bangkok, Thailand, (259% increase) followed by Seoul, Republic of Korea, (173% increase) and Tianjin, China (147%). Of the cities included in this study, the fastest growing were found predominantly in Asia, whereas the slower growing cities were almost all located in North America and Western Europe.
As shown in figure 3, by 2006 the Settlements area of Bangkok was already more than 2.5 times what it was in early 1988. In fact, over the past several decades, rural-urban migration has been a major contributor to urban growth in Bangkok, where the population increased from 2.1 million in 1960 to 6.3 million in 2000 (Choiejit and Teungfung 2005).

Grid-cell-based urban structure analysis
To investigate the relationships among the land-cover changes caused by urban expansion, we calculated the correlation coefficients between Settlements and the other land-cover categories (i.e., Forest, Cropland, Grassland and Wetlands) based on the 1 km 2 grid square cells. Table 3 presents a summary of the linear correlation coefficient matrix for the changes of land-cover categories in the 50 cities during the study period. Table 3 shows that the Settlements changes had a strong negative correlation with changes of greenfields (i.e., Forest, Cropland and Grassland) and Wetlands, which implies that the enormous expansion of Settlements area resulted in rapid decreases in green spaces and wetland. Note the following in particular: As shown in table 3, the differences in correlation coefficients of Settlements changes between Forest, Cropland and Grassland suggest that the availability of land for urban expansion varied amongst cities. However, some of the unrealistically high or low correlation coefficient values were possibly due to the fact that the land cover classification results affected by changes in observation time. Figure 4(a) shows the grid-cell-based spatial change of Settlements area from 1984 to 2010 in Beijing. The value of each grid square cell was calculated by subtracting the Settlements area of 1984 from that of 2010 in each grid cell and then dividing the changed area by the cell area. As figure 4(a) illustrates, the Settlements area rapidly expanded mainly in three directions: north, east and south. Figure 4(b) shows the grid-cell-based spatial change of Cropland area from 1984 to 2010 in Beijing, which was calculated in the same way as for figure 4(a). The Cropland area decreased rapidly around the core city, typically in three directions: north, east and south. As shown in figure 2, and figures 4(a) and (b), the Cropland area has been mainly converted to the urban/built-up area and Grassland. Results from our correlation analyzes also show that Cropland changes was negatively correlated with Grassland changes (r = −0.59) in Beijing. The conversion of Cropland to Grassland was partly due to several government initiatives such as the Desertification Combating Program around Beijing and Tianjin, and the Conversion of Cropland to Forest and Grassland Program (also known as the 'Grain for Green' Program) (Liu et al 2010). Upon completion of these programs, part of the Cropland was converted to urban land use, grassland and shrubs (Hu 2007). Another reason is that Cropland is temporarily converted to bareland when the land is cleared to prepare the site for construction and this bareland would usually be covered with weed grass during the transition period.
As figures 4(a) and (b) illustrate, in Beijing, the Settlements area rapidly expanded to the surrounding suburban area (north, east and south parts) where it was mainly flat; in contrast, the urban/built-up area decreased in the centre of the city. In fact, many high-rise (e.g., office or commercial) Brussels1987 Brussels2010  buildings replaced the dense low-rise buildings in the city centre concurrently with the development of housing estates in the suburbs. Meanwhile, the greatest decreases in agricultural land area took place in areas suitable for urban land use where the demand for land for urban purposes is high. This suburbanization process expanded urban areas enormously over a short time. This finding is consistent with previous studies of Tokyo, Japan, where urban areas expanded at the expense of agricultural land (Bagan and Yamagata 2012). Figures 4(c) and (d) show the spatio-temporal dynamics of Settlements and Grassland changes in Los Angeles, and figures 4(e) and (f) shows them for Settlements and Wetlands changes in Seoul. Figures 4(c)-(f) are also calculated in the same way as for figure 4(a). In Los Angeles, urban area rapidly expanded to suburban area causing either fragmentation or a complete loss of Grassland (figures 4(c) and (d)). Conversions to the developed area from grasslands/shrublands in the Central California Valley are also reported by Sleeter et al (2011). In Seoul from 1986 to 2011, the Wetland area rapidly decreased around the human development area (figures 4(e) and (f)). Figure 5(a) shows the relationship between Settlements and Cropland changes from 1984 to 2010 in Beijing. We found a strong, negative linear relationship (r = −0.73), suggesting that a vast area of Cropland had been converted to urban/built-up area during the last two decades in Beijing. Figures 5(e) and (f) show the relationships between Settlements and Wetlands changes in Seoul and Bangkok, respectively. The correlations are very weak in both cities, implying that Wetlands have been converted to other landcover types. In fact, in Seoul, the correlation between Table 3. Summary of the linear correlation coefficients between changes of settlements and changes of other land-cover categories (Forest, Cropland, Grassland and Wetlands).
Correlation between settlements and other land cover types. The correlation coefficient is denoted by r.  Wetlands and Grassland changes is −0.45, which implies that part of Wetlands has been converted to Grassland. As figures 3-5 illustrate, the loss of green spaces was spatially concentrated in human development areas. These results confirm that the tendency to continuously increase the amount of land with residential and industrial facilities and rapid urban development directly cause loss of farmland and other vegetation on the urban periphery.
The rapid growth in global urban land cover is likely to continue as long as urban populations continue to grow. The future expansion of cities into arable lands causes concern that the displacement of agricultural land by urban land use will require, where possible, bringing new land into cultivation as well as increasing land productivity (Angel et al 2011).

Conclusion
The combination of remote sensing data and grid cells analysis enabled us to quantify and compare the spatial and temporal forms, patterns and structures of urban growth for 50 cities across the globe. The urban growth rate of almost half of the cities was greater than 50%, with much of this increase happening in suburbs. The rapid expansion of urban land use area around cities accompanied the massive decreases in green spaces that occurred during the study period of around 1985-2010. This resulted in a reduction of natural habitat and a sharp decrease in the biodiversity of the cities.
The spatio-temporal analysis based on 1 km 2 grid cells revealed that there are considerable differences among the current urban growth trends in many cities in the world. The use of grid cells makes it possible to examine the extent to which different cities fall into different urban expansion types, so that the commonalities and differences among them can be better understood. The results from urban growth analyzes based on grid cells should be of great importance to city planners, environmental managers and policy makers. Urban expansion processes differ between areas with different levels of human development. For example, in Tokyo there has been large population loss in the city centre as many residents have migrated to the outlying suburban areas (Bagan and Yamagata 2012), whereas the population has tended to increase in urban centres in India (Taubenböck et al 2009) and China (Seto and Fragkias 2005). Thus, further research is also needed to investigate the relationships between population density and settlement density to inform future land development and conservation strategies. In addition, we have plans to combine urban land cover changes with the greenhouse gas emissions analysis and to extend the full analysis-greenhouse gas emissions and urban land cover changes-to other regions. Taubenböck H, Wegmann M, Roth A, Mehl H and Dech S 2009 Urbanization in India-spatiotemporal analysis using remote sensing data Comput. Environ. Urban Syst. 33 179-88 Townshend J R and Justice C O 2002 Towards operational monitoring of terrestrial systems by moderate-resolution remote sensing Remote Sens. Environ. 83 351-9 Xie Y, Sha Z and Yu M 2008 Remote sensing imagery in vegetation mapping: a review J. Plant Ecol. 1 9-23 Zhang Q and Seto K C 2011 Mapping urbanization dynamics at regional and global scales using multi-temporal DMSP/OLS nighttime light data Remote Sens. Environ. 115 2320-9 Environ. Res. Lett. 9 (2014) 064015 H Bagan and Y Yamagata