Tracking Urban Expansion Through Satellite Imagery

December 12, 2023

Tracking Urban Expansion Through Satellite Imagery

African cities will be the epicenter of urbanization by the middle of this century. The best estimates project a doubling of the population by 2050 and two-thirds of that increase will be concentrated in cities. Because of this, there is some urgency for countries and cities to engage in urban planning in a way that helps foster economic development. However, for low- and middle-income countries even as large as Nigeria, there is a dire lack of reliable and frequent population counts or censuses and limited direct management of urban development.

Satellite imagery can help us in this regard by revealing how development has occurred over time, and identifying likely areas for future development over longer periods of time than are covered by more traditional data sources. Satellite imagery has shown its usefulness for development economics, as it can help us study key economic questions even in settings where data is not easily available by relying on information captured from space.

In this post, I will show how I have used remote sensing along with supervised classification to look at the evolution of Nigerian cities between 2000 and 2020. The random forest classification allows me to identify land cover types, more specifically the built environment from basic satellite images and using tools from by Ujaval Gandhi’s Spatial Thoughts, we can use this classification to compute what areas and land use have been lost to growing urbanization in the time period studied.

Finding the right satellite images

There are many, many satellite images available - one easy way to access and get to know them is through the GEE Code Editor’s data catalog. For each of the rasters, the GEE data catalog identifies the specific bands contained in the data and also points to potential issues and needed corrections in the images selected.

Here, I used LandSat 7, which goes back to 2000. In order to have the best image possible for a given year of data, I correct for any of the issues mentioned above (such as correcting surface reflectance issues or removing clouds) and then aggregate images from a given year by taking the median of these images. This allows me to obtain as clear an image as possible.

Running the supervised classification

The goal of this exercise is to go from this satellite image above to an image that identifies the specific land use for each of the 30m pixels in the image. For this, I use a supervised classification, meaning that I tell the machine learning algorithm the different categories that I want identified in the new images based on an existing categorization. Because I am focused here on the built environment only, I only have 4 categories: bare (untouched land), vegetation, urban, and water.

Existing land cover data usually have more specific categories, such as the National Land Cover Database for the US that has 16 categories including 4 for urban land, and/or are collected on a finer spatial and temporal resolution like the daily, 10m resolution Dynamic World database. Unfortunately, such sources as Dynamic World only start from more recent years (2015 onwards) and richer land cover maps are more likely to cover higher-income countries.

Creating the categories

Instead of relying on existing training data here, I manually create my own classification by sampling observations for each of the categories of interest (bare, urban, water, and vegetation) in the original image from Landsat 7 in 2020. To do this, I use pins and/or small polygons in the GEE map window; for each category, I create a new geometry to which I assign a specific value but give a similar class name for all the geometries of the 4 classes. To make sure my sampling is accurate, I also compare the same object (here the city of Lagos) during different seasons.

After building my 4 geometries for bare, urban, water, and vegetation, they can be merged into one variable that will be used by the classifier.

Improving on the image

The classifier is a random forest classifier with 50 trees using the merged geometries as a training sample. Here, instead of using the LandSat 7 median image constructed above along with my collection of observations for each category, I also include more information to the image in order to improve on the classification. 

In this case, since I am interested in distinguishing between water, vegetation, and built environment, I include additional bands composed of indexes that help highlight these different features, such as NDVI, NDBI, mNDWI and BSI. These indices are normalized differences between the values of the different bands contained in the image: for instance, NDVI (normalized difference vegetation index) is the normalized difference between the values of the NIR (near-infrared) and red bands. It is most commonly used in remote sensing as an indicator for greenness. This relies on the fact that healthy vegetation will absorb the red band and reflect off the near-infrared band, so high values of the index (close to 1) indicate vegetation. The other normalized differences use the same principle to identify water (modified NDWI), built environment (NDBI and built-up surface index, which is also a normalized difference but one that combines multiple bands: near-infrared and blue on one side, and red and short-wave infrared 1 on the other). 

Additionally, bands that track other parameters such as slope can also be included to the image that will be classified. All these are added on the original bands of the median image and the resulting image is denoted below as composite

var addIndices = function(image){

  var ndvi = imagenormalizedDifference(["SR_B4","SR_B3"]).rename(["ndvi"]);

  var ndbi = imagenormalizedDifference(["SR_B5","SR_B4"]).rename(["ndbi"]);

  var mndwi = imagenormalizedDifference(["SR_B2","SR_B5"]).rename(["mndwi"]);

  var bsi = imageexpression(

"(( X + Y ) - (A + B)) /(( X + Y ) + (A + B))",{

'X': imageselect("SR_B5"),//swir1

'Y': imageselect("SR_B3"),//red

'A': imageselect("SR_B4"),// nir

'B': imageselect("SR_B1"),// blue

}).rename("bsi");

return imageaddBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi);

};


var composite = addIndices(landsat7_2020);

Running the classifier

This is better explained in the code below, but we use machine learning (a random first classifier below) to detect land cover. Using the training sample we constructed above, and based on the starting image augmented with the additional bands, the random forest classifier will try to detect (with some error) in a new image the same categories that were identified in the training sample. For the researcher, the only choice variables here are the number of trees, where each tree can be thought of as a data sample drawn from the image to apply the algorithm.

// Create a random forest classifier with 50 trees.

var classifier = eeClassifiersmileRandomForest({

numberOfTrees:50

}).train({// Train the classifier.

// Use the sampled pixels.

features: training,

// This is the class that we want to predict.

classProperty:"class",

// The bands the classifier will use for training and classification.

inputProperties: compositebandNames()

});


// Apply the classifier on the 2020 image.

var classified20 = compositeclassify(classifier);

After running the classifier, we can create a new map showing the classified image. The end result is shown below, where red denotes the built environment, green denotes vegetation and blue denotes water from the 2020 composite image. Since we used the pixels from Lagos in 2020 to train the classifier, this classification is not as interesting as looking at what happened in 2010 and in 2000, and where changes in built environment are happening. We can then apply this same classifier to classify images from 2010 and 2000.

Following images show Lagos in 2010 (top) and in 2000 (bottom), based on the random forest classification and the 2020 training sample.

To more systematically look at changes, we can map the difference between the two maps. Recall that each pixel corresponds to a specific value (bare, urban, water, vegetation) thanks to the classification. There are multiple ways to do this: the simplest way is to tabulate the difference between 2010 and 2020 and map it. We can obtain a map,showing only the areas that have changed their classification between 2010 and 2020.

Since this does not tell us much about the changes in the built environment or even what types of land use has been moved to urban land, we can try to pin down exactly which land has become urban over time (see code below). We can further use the function ee.Reducer.frequencyHistogram() to create a transition matrix, with the caveat that the matrix will provide us with a count (the number of pixels that changed from one category to the other), which might not be useful information.

## since we had previously assigned the value of 0 to one of the categories, we reassign new values here

var classified10_new = classified10remap([0,1,2,3],[1,2,3,4]);

var classified20_new = classified20remap([0,1,2,3],[1,2,3,4]);


## this allows us to keep track of class in 2010 and in 2020, by multiply the class value in 2010 by 100

## the first digit gives us class in 2010 and the last one class in 2020

var envt2020 = classified10_newmultiply(100).add(classified20_new).rename("transitions20");


## we use this new variable to identify land that was not urban in 2010 (first number is not 2) and has become urban by 2020 (units digit is equal to 2)

var urbandiff20_10 = envt2020eq(102).or(envt2020eq(302)).or(envt2020eq(402));

## this second variable tracks the land that remained urban between 2010 and 2020

var urban20102020 = envt2020eq(202);



## this adds the map to visualize the changes, as seen below

MapaddLayer(urbandiff20_10clip(lagos),{min:0,max:1,palette:["white","#7a0177"]},"newly urban - in 2020");

Wrapping up the classification

In order to properly finish the classification, we want to check the accuracy, that is whether the classifier has correctly identified the different categories in our classified images.

Usually this can be done by taking the training dataset, splitting it into two parts: one part is used to train the classifier and the second part is used as a validation sample. In this case, since the training data was created directly, and there is a dearth of landcover maps going back to 2000 for Nigeria, the solution is to find another source of information that will act as validation.

References

  1. https://spatialthoughts.com/

  2. https://www.linkedin.com/pulse/ndvi-ndbi-ndwi-calculation-using-landsat-7-8-tek-bahadur-kshetri

  3. Lall, Somik Vinay, J. Vernon Henderson, and Anthony J. Venables. 2017. “Africa’s Cities: Opening Doors to the World.” World Bank, Washington, DC. License: Creative Commons Attribution CC BY 3.0