Introduction
Built-up places can be found almost everywhere nowadays. The majority of built-up places are found on land. Thus, it is a no-brainer why a majority of remote sensing applications are focused on land surfaces. Built-up refers to areas, dominated (more than 50% in cover) by non-vegetated, human-constructed elements, such as roads, buildings, runways and industrial facilities.
Going by the definition, built-up area assessment is a key factor in urban studies, where there is a high prevalence of built-up surfaces compared to rural areas. Just as there are many vegetation indices such as Normalized Difference Vegetation Index (NDVI), Normalized Difference Moisture Index (NDMI) and Enhanced Vegetation Index (EVI) for vegetation studies, urban studies employ just as many in assessing built-up cover. These include the Normalized Difference Built-up Index (NDBI), Index-Based Built-up Index (IBBI), Urban Index (UI) and Built-up Index (BUI). The last is what shall be our major focus but not before passing a stealing glance at NDBI.
In this exercise, we shall explore remote sensing applications. We shall use the BUI to extract the built-up area in Sao Paulo, Brazil. BUI is a relatively more accurate index compared to NDBI. It is also relatively easy to calculate. In the BUI, positive values indicate built-up areas. In all built-up indices, a major problem that still exists is that it is hard to delineate both built-up and bare areas since the reflectance characteristics are quite similar. However, we cross our fingers bare areas will not influence our results too much since Sao Paulo is a relatively green city! At least from a bird’s eye view of course!
Let’s load an administrative shapefile for Sao Paulo to work with. The shapefile is accessible from here.
Let’s first begin from R. Then we shall go to GEE for Google Earth Engine.
As always, the usual spatial R packages have to be loaded first.
library(sp)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/gachuhi/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/gachuhi/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(rgeos)
## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-0
## Polygon checking: TRUE
library(raster)
Let’s read the downloaded Sao Paulo Municipalities shapefile into R Studio.
# Read the Sao Paulo municipalities shapefile into R Studio
sao <- readOGR(dsn="E:/documents/gis800_articles/building_detection/sao_paulo/sao_paulo.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "E:\documents\gis800_articles\building_detection\sao_paulo\sao_paulo.shp", layer: "sao_paulo"
## with 650 features
## It has 16 fields
With that done, we have saved the shapefile into the object sao. Time to view it.
# View the sao object
plot(sao)
Wow. Sao Paulo isn’t just a small city, as presumed. It’s a whole region! However, we don’t want to work with this conglomeration of municipalities. We just want a single unitary shapefile, you know, for clipping a raster into. The above sao object is only good if were are planning to clip a raster, say a satellite image or DEM into various municipalities.
To create a single unitary shapefile where all municipalities have been merged into a single region, we have to aggregate them. In our case, the municipalities will be aggregated into a single attribute that is common to all of them. To do this, let’s first view the attributes of sao and choose a field value that is common to the hundreds of municipalities in this region.
# View the first few attributes of sao
sao@data[0:15, ]
## CODMUNIC NOMEMUNICP NOMEMUNIC CODUF NOMEUF_1
## 0 <NA> <NA> <NA> <NA> <NA>
## 1 <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA>
## 5 3500105 ADAMANTINA Adamantina 35 SAO PAULO
## 6 3500204 ADOLFO Adolfo 35 SAO PAULO
## 7 3500303 AGUAI Agua\xa1 35 SAO PAULO
## 8 3500402 AGUAS DA PRATA \xb5guas da Prata 35 SAO PAULO
## 9 3500501 AGUAS DE LINDOIA \xb5guas de Lind\xa2ia 35 SAO PAULO
## 10 3500550 AGUAS DE SANTA BARBARA \xb5guas de Santa B\xa0rbara 35 SAO PAULO
## 11 3500600 AGUAS DE SAO PEDRO \xb5guas de S\xc6o Pedro 35 SAO PAULO
## 12 3500709 AGUDOS Agudos 35 SAO PAULO
## 13 3500758 ALAMBARI Alambari 35 SAO PAULO
## 14 3500808 ALFREDO MARCONDES Alfredo Marcondes 35 SAO PAULO
## NOMEUF_2 REGIAO_1 REGIAO_2 CODMESO NOMEMESO CODMICRO
## 0 <NA> <NA> <NA> <NA> <NA> <NA>
## 1 <NA> <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA> <NA>
## 5 SP SUDESTE SE 08 PRESIDENTE PRUDENTE 35
## 6 SP SUDESTE SE 01 SAO JOSE DO RIO PRETO 04
## 7 SP SUDESTE SE 07 CAMPINAS 29
## 8 SP SUDESTE SE 07 CAMPINAS 30
## 9 SP SUDESTE SE 07 CAMPINAS 33
## 10 SP SUDESTE SE 04 BAURU 22
## 11 SP SUDESTE SE 06 PIRACICABA 28
## 12 SP SUDESTE SE 04 BAURU 20
## 13 SP SUDESTE SE 11 ITAPETININGA 42
## 14 SP SUDESTE SE 08 PRESIDENTE PRUDENTE 36
## NOMEMICRO CODUFMESO CODUFMESOM AREA_97 ID
## 0 <NA> <NA> <NA> <NA> 1
## 1 <NA> <NA> <NA> <NA> 2
## 2 <NA> <NA> <NA> <NA> 3
## 3 <NA> <NA> <NA> <NA> 4
## 4 <NA> <NA> <NA> <NA> 5
## 5 ADAMANTINA 3508 350835 411,7792 6
## 6 SAO JOSE DO RIO PRETO 3501 350104 210,8163 7
## 7 PIRASSUNUNGA 3507 350729 473,3618 8
## 8 SAO JOAO DA BOA VISTA 3507 350730 142,5886 9
## 9 AMPARO 3507 350733 59,99524 10
## 10 AVARE 3504 350422 408,5039 11
## 11 PIRACICABA 3506 350628 3,64473 12
## 12 BAURU 3504 350420 967,5698 13
## 13 ITAPETININGA 3511 351142 159,2092 14
## 14 PRESIDENTE PRUDENTE 3508 350836 119,5107 15
Looks like the value SUDESTE under the REGIAO_1 field is common to all, if not most, municipalities. It’s like the attribute of a Province in a Kenya Counties shapefile. If one would like to merge all counties, say to their former respective provinces, the Province attribute would do the trick.
Let’s use the aggregate function from the sp package to aggregate the municipalities into one shapefile. The REGIAO_1 field.
# Aggregate municipalities using the REGIAO_1 field
sao_region <- aggregate(sao, by = "REGIAO_1", do_union = T)
# Plot the aggregated shapefile
plot(sao_region)
Done! We know that a few municipalities did not have any value under the REGIAO_1 and these may have still retained their boundaries. However, they are so small they are hardly visible. Methinks they could have been islands.
# Check attributes of the aggregated shapefile
sao_region@data[0:5, ]
## [1] NA "SUDESTE" NA NA NA
Let’s save our aggregated shapefile from R into our computer’s directory. Thereafter, we shall transfer it to Google Earth Engine (GEE) where the bulk of our work shall reside.
# Save the aggregated shapefile into your directory
sao_region_shp <- writeOGR(sao_region, dsn = ".", layer = "sao_region",
driver = "ESRI Shapefile", overwrite_layer = T)
After saving your shapefile, head over to your GEE account.
GEE
First things first. In your GEE account, head over to the assets tab.
Click on the New drop-down appearing with the yellow halo above. Select Shapefiles from the Table upload section. You will get a dialog box like the one below.
Select the sao_region_shp shapefile from your directory ensuring you have selected the dependencies such as .shp, .shx, .dbf and .prj as explained in the dialog box.
Once done, click Upload. The shapefile will appear under your legacy assets. In order to import it into your GEE script, just click the arrow with the hover text ‘Import into script’. The shapefile should appear in your GEE code editor. Give it the name sao_region.
For assessing the built-up area in Sao Paulo, we created the following variables.
- city
- point1
- point2
Below is the code that will load all the above variables, including that of sao_region into your GEE code editor.
var sao_region = ee.FeatureCollection("users/sammigachuhi-test/sao_region"),
city =
/* color: #d63000 */
/* shown: false */
/* displayProperties: [
{
"type": "rectangle"
}
] */
ee.Geometry.Polygon(
[[[-46.949648518267026, -23.182191189569593],
[-46.949648518267026, -23.842046005948035],
[-46.103701252642026, -23.842046005948035],
[-46.103701252642026, -23.182191189569593]]], null, false),
point =
/* color: #98ff00 */
/* shown: false */
ee.Geometry.Point([-46.598294953027256, -23.677920629371847]),
point2 = /* color: #0b4a8b */ee.Geometry.Point([-46.598086018267026, -23.33359429234038]);
Built-up area assessment
To undertake our built-up assessment, let’s first zoom in on Sao Paulo.
// zoom to shapefile of Sao Paulo region
Map.centerObject(sao_region, 6);
Let’s add our sao_region shapefile into the code editor and view it.
// Add the Sao Paulo shapefile to google earth engine map canvas
Map.addLayer(sao_region, {}, "sao_paulo");
Remember the two points that we had created, point1 and point2. They are about to point us somewhere. The two points are crucial in downloading Sentinel-2 images to those specific locations. Let’s go ahead and load our Sentinel-2 images into GEE.
// Load Sentinel 2 image for clearest month of the year which is August
var sao_clear = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterBounds(point)
.filterDate('2021-08-01', '2021-08-30')
.sort('CLOUDY_PIXEL_PERCENTAGE')
.first();
var sao_clear2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterBounds(point2)
.filterDate('2021-08-01', '2021-08-30')
.sort('CLOUDY_PIXEL_PERCENTAGE')
.first();
print("Images over Sao Paulo July - September 2021: ", sao_clear);
print("Images over Sao Paulo July - September 2021: ", sao_clear2);
Let’s add our Sentinel 2 images to the GEE canvas.
Map.addLayer(sao_clear, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021");
Map.addLayer(sao_clear2, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021vs2");
Below are our two images. They in fact line up to each other so well that there is absolutely no space between them. However, they are separate. Unless we look for a way to join them, it will be double work to assess the built-up area in both of them.
What shall we do? Mosaic! Like, more psyche!!
Mosaicing. Our two images may not overlap, but the technique will be the same. In GEE, mosaicking is the technique of combining several images with overlapping parts (the images may be obtained at different times, different viewing angles or different sensors) into a largescale seamless high resolution image. Our two images may not be overlapping each other, but the technique will be the same. In GEE, mosaicking is done using the mosaic() method.
Before mosaicking the two images, however, we have to first combine both to form a single collection using the ee.ImageCollection() method. Thereafter, the mosaicking is done on the image collection. This two-step procedure is done using the code script below.
// Mosaic sao_clear and sao_clear2
var sao_images = ee.ImageCollection([sao_clear, sao_clear2]);
var sao_images2 = sao_images.mosaic();
Map.addLayer(sao_images2, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021 Mosaicked");
Remember the city bounding box we had created earlier? Its purpose was to serve as a bounding box to clip the larger mosaicked image to the city core. We will clip the larger mosaicked image, which is too large an area for built-up assessment anyway, to the city core.
Let’s go ahead and clip the mosaicked image – sao_images2 to the city core.
var sao_images_clip = sao_images2.clip(city);
Map.addLayer(sao_images_clip, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021 Clipped");
Calculating the Normalized Difference Built-up Index (NDBI)
Now that we have our image reduced to a sufficient scale for our analysis, we shall first perform NDBI on it followed by other calculations that will build up to BUI. We calculate NDBI using the formula shown below:
NDBI = (SWIR – NIR)/(SWIR + NIR)
GEE’s ee.Image expression method simplifies mathematical expressions involving images. We define the mathematical formula using the expression parameter and select the images to use in the expression using the map parameter. We use the above method to calculate NDBI like so:
var ndbi_expression = '(SWIR - NIR)/(SWIR + NIR)';
var ndbi2021 = sao_images_clip.expression({
expression: ndbi_expression,
map: {
SWIR: sao_images_clip.select("B11"),
NIR: sao_images_clip.select("B8")
}
}).rename("ndbi");
Map.addLayer(ndbi2021, {palette:['blue', 'yellow', 'red']}, "NDBI 2021");
We have also defined the bandnames of the output images to read “ndbi” using the rename method.
Calculating the Normalized Difference Vegetation Index
NDVI may be confusing you when we are primarily dealing with built-up areas. But remember, we previously mentioned that we calculate BUI by subtracting NDVI from NDBI. That is, BUI = NDBI – NDVI.
Let’s calculate the NDVI speedily.
The formula for NDVI is (NIR – Red)/(NIR + Red). We shall use the image.expression() function to calculate the NDVI just like we did for NDBI.
// Calculate the Normalized Difference Vegetation Index (NDVI) for sao_images_clip
var ndvi_expression = '(NIR - Red)/(NIR + Red)';
var ndvi2021 = sao_images_clip.expression({
expression: ndvi_expression,
map: {
NIR: sao_images_clip.select('B8'),
Red: sao_images_clip.select('B4')
}
}).rename('ndvi');
Map.addLayer(ndvi2021, {palette: ['red', 'yellow', 'green']}, "NDVI 2021");
The NDVI image looks lovely!
Now without further ado, let’s calculate the BUI for Sao Paulo. Sao Paulo looks to be located in a lush green environment. Therefore, bare areas should not constitute much of a large percentage of our end results. Therefore, if the area’s greenery favours us, the BUI should mostly collect the impervious surfaces.
Let’s do this thing. BUI formula is NDBI – NDVI.
// Calculate the Built Up Index (BUI)
// BUI = NDBI - NDVI
var bui2021 = ndbi2021.subtract(ndvi2021).rename('bui');
Map.addLayer(bui2021, {palette:['green', 'yellow', 'red']}, 'BUI 2021');
Notice that towards the southeast there seems to be a high reading of built-up areas (more positive values indicate built-up areas) in yellow. However, on closer inspection, this yellow patch is in fact a lake! It seems that even though bare areas are rare in Sao Paulo, the BUI also captures water values as built-up too! The explanation could be that built-up and water register low NIR values, and thus are difficult to distinguish from each other. The same issue is also, albeit to a lesser extent, prevalent in the NDBI image. Oh well, we shall work with what we have.
In BUI and NDBI, higher positive values indicate built-up and barren areas. In our case, barren areas are almost non-existent. As such, built-up areas carry the most part. Since we want to only extract the most probable built-up area(s) from our bui2021 image, we shall work with a threshold of 0. Extracting all pixels with a value above this threshold should comfortably give us an image of built-up only, with the additional problem of water bodies also captured therein.
First, let’s get the minimum and maximum values for the bui2021 image.
// Get min max of bui2021
var bui2021_minmax = bui2021.reduceRegion({
reducer: ee.Reducer.minMax(),
geometry: city,
scale: 10,
bestEffort: true
});
print("Min max values for bui2021: ", bui2021_minmax);
print("BUI for 2021", bui2021);
The values of both extremes are:
bui_max: 1.867632850241546
bui_min: -1.741782898901804
Let’s select all bui2021 pixels above a threshold of 0.
// Select all areas above a threshold of 0
var datamask2021 = bui2021.select("bui");
var mask2021 = datamask2021.gte(0);
var built2021 = bui2021.updateMask(mask2021);
Map.addLayer(built2021, {palette:['red']}, "Built 2021");
Time to calculate the total built-up area within Sao Paulo’s city bounds. The pixelArea() function calculates the area of each pixel in square metres. Thus, calculating the total sum of the areas of all pixels–which are in square metres–will give us the area of that image in…what? …square metres.
// Calculate area of the Built-up area (BUI) in 2021
var builtArea2021 = built2021.multiply(ee.Image.pixelArea());
var stats2021 = builtArea2021.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: city,
scale: 10,
bestEffort: true
});
Let’s print the total built-up surface area within the city bounds we had set for Sao Paulo.
print("Total impervious surface area is: ", stats2021);We get a value of 107454854.88727954 or roughly 107 square kilometres. However, these values are false since we also captured waterbodies. For example, the water feature encircled in yellow in the figure below. In order to get the most authentic built-up area values, one will need to mask out the waterbodies.
The collection of water pixels is a disappointment, but the BUI calculation did a fair job of extracting built-up areas and leaving out vegetated surfaces.
Conclusion
Several urban studies have used built-up indices, assessing the built-up cover and studying the contribution of buildings to the urban heat island effect. However, the confusion about bare soil, or barren areas with built-up cover has led to the development of new indices that take both into account. Examples are Dry Built-up Index (DBI) and Dry Bare-Soil Index (DBSI). In our case, we didn’t find barren areas prevalent in our region of interest, but we counted water features as built-up.
Similar reflectance values mostly cause confusion about land covers, therefore we must conduct additional processes to get the true result. Of late, machine learning has been promising. For example, instead of collecting all built-up areas indiscriminately, machine learning looks for specified patterns and collates them to the end result. For example, if one only wanted to collect building features from an image, machine learning would look for rectangular and other defined polygonal patterns. We extract rectangular and clear-cut polygon shapes with matching reflectance characteristics as features. This process eliminates roads from the end result, thus giving the researcher a true estimate of buildings in an area.
Here is the full code script for this exercise.
var sao_region = ee.FeatureCollection("users/sammigachuhi-test/sao_region"),
city =
/* color: #d63000 */
/* shown: false */
/* displayProperties: [
{
"type": "rectangle"
}
] */
ee.Geometry.Polygon(
[[[-46.949648518267026, -23.182191189569593],
[-46.949648518267026, -23.842046005948035],
[-46.103701252642026, -23.842046005948035],
[-46.103701252642026, -23.182191189569593]]], null, false),
point = /* color: #98ff00 */ee.Geometry.Point([-46.598294953027256, -23.677920629371847]),
point2 = /* color: #0b4a8b */ee.Geometry.Point([-46.598086018267026, -23.33359429234038]);
//////////////////////////////
// zoom to shapefile of Sao Paulo region
Map.centerObject(sao_region, 6);
// Add the Sao Paulo shapefile to google earth engine map canvas
Map.addLayer(sao_region, {}, "sao_paulo");
// Load Sentinel 2 image for clearest month of the year which is August
var sao_clear = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterBounds(point)
.filterDate('2021-08-01', '2021-08-30')
.sort('CLOUDY_PIXEL_PERCENTAGE')
.first();
var sao_clear2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterBounds(point2)
.filterDate('2021-08-01', '2021-08-30')
.sort('CLOUDY_PIXEL_PERCENTAGE')
.first();
print("Images over Sao Paulo July - September 2021: ", sao_clear);
print("Images over Sao Paulo July - September 2021: ", sao_clear2);
Map.addLayer(sao_clear, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021");
Map.addLayer(sao_clear2, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021vs2");
// Mosaic sao_clear and sao_clear2
var sao_images = ee.ImageCollection([sao_clear, sao_clear2]);
var sao_images2 = sao_images.mosaic();
Map.addLayer(sao_images2, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021 Mosaicked");
// Clip the Mosaicked image to the variable city bounds
var sao_images_clip = sao_images2.clip(city);
Map.addLayer(sao_images_clip, {min: -67, max: 1700, bands:['B4', 'B3', 'B2']}, "Sao Paulo 2021 Clipped");
// Calculate Normalized Difference Built-up Index (NDBI) for sao_images_clip
var ndbi_expression = '(SWIR - NIR)/(SWIR + NIR)';
var ndbi2021 = sao_images_clip.expression({
expression: ndbi_expression,
map: {
SWIR: sao_images_clip.select("B11"),
NIR: sao_images_clip.select("B8")
}
}).rename("ndbi");
Map.addLayer(ndbi2021, {palette:['blue', 'yellow', 'red']}, "NDBI 2021");
// Calculate the Normalized Difference Vegetation Index (NDVI) for sao_images_clip
var ndvi_expression = '(NIR - Red)/(NIR + Red)';
var ndvi2021 = sao_images_clip.expression({
expression: ndvi_expression,
map: {
NIR: sao_images_clip.select('B8'),
Red: sao_images_clip.select('B4')
}
}).rename('ndvi');
Map.addLayer(ndvi2021, {palette: ['red', 'yellow', 'green']}, "NDVI 2021");
// Calculate the Built Up Index (BUI)
// BUI = NDBI - NDVI
var bui2021 = ndbi2021.subtract(ndvi2021).rename('bui');
Map.addLayer(bui2021, {palette:['green', 'yellow', 'red']}, 'BUI 2021');
// Get min max of bui2021
var bui2021_minmax = bui2021.reduceRegion({
reducer: ee.Reducer.minMax(),
geometry: city,
scale: 10,
bestEffort: true
});
print("Min max values for bui2021: ", bui2021_minmax);
print("BUI for 2021", bui2021);
// Select all areas above a threshold of 0
var datamask2021 = bui2021.select("bui");
var mask2021 = datamask2021.gte(0);
var built2021 = bui2021.updateMask(mask2021);
Map.addLayer(built2021, {palette:['red']}, "Built 2021");
// Calculate area of the Built-up area (BUI) in 2021
var builtArea2021 = built2021.multiply(ee.Image.pixelArea());
var stats2021 = builtArea2021.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: city,
scale: 10,
bestEffort: true
});
print("Total impervious surface area is: ", stats2021);
Remote Sensing Applications: Built-up Assessment