Investigating Land Use and Land Cover Change in Bahrain: 1987-2013

I. Introduction

At the request of the Bahrain Center for Human Rights (BCHR), the Geospatial Technologies and Human Rights Project of the American Association for the Advancement of Science (AAAS) used a time series of satellite imagery to examine land use and land cover changes that have taken place in The Kingdom of Bahrain since 1987. BCHR had expressed concern regarding environmental change occurring in the country, particularly the impact of such change on non-affluent citizens, including with regard to access to water and the coastline; massive urbanization, including land reclamation projects; and changes to vegetation.

The Kingdom of Bahrain is located in the Persian Gulf between Saudi Arabia and Qatar and is comprised of 36 islands (Figure 1). The country was historically known for its palm trees and natural springs.¹ The largest of the islands, al-Awal Island, contains the capital city, Manama, and is home to the majority of Bahrain’s population and economic activity. This study focuses on al-Awal Island and those islands immediately adjacent to it.

Over the past 40 years, the population of Bahrain has increased from 213,000 in 1970 to 1,252,000 in 2010- an increase of 487%.² This rapid population growth has strained the ecological resources of the island. For example, the total amount of water used annually in Bahrain grew from 138 million m³ in 1980 to 269 million m³ in 2000.³ In addition, multiple land reclamation projects have substantially increased the size of the country while negatively impacting the environment. Reports also indicate that al-Awal Island and the islands immediately adjacent have lost a substantial portion of their native palm trees.⁴

II. Data and Methods

Time series analyses of land use and land cover change require extensive data preparation to ensure the comparability of the imagery. Images must be geometrically corrected so that the images are properly aligned at the pixel level. If this is not done, detected changes may reflect alignment errors rather than actual changes in land use and land cover. In addition, atmospheric effects must be accounted for. There can be substantial variation in atmospheric absorption depending on local meteorological factors at the time of image acquisition. As a result, radiometric correction techniques that adjust for atmospheric effects have been developed. The following section describes the images acquired and corrections applied.

For this study, AAAS acquired images from NASA’s Landsat series of satellites, which have been collecting imagery for the past 40 years. Though over its lifetime there have been some problems with data retention, acquisition, and transfer, the Landsat program provides the most comprehensive and long-term data set of land cover available today. As a result, it has been used extensively to monitor land use and land cover changes around the world.

The Landsat program began in 1973 with the launch of the Earth Resources Technology Satellite 1, which was later renamed Landsat-1. Landsat-1 carried a 4-channel multispectral scanner (MSS), which acquired imagery at a 60m (meter) spatial resolution. Landsat-1 was followed by Landsat-2 in 1975 and Landsat-3 in 1978. Both satellites carried the MSS sensor. In 1982, Landsat-4 was launched. In addition to the MSS sensor, Landsat-4 carried the new 7-channel Thematic Mapper (TM) sensor, which acquired imagery at a 30m spatial resolution. Next, in 1984, Landsat-5 was launched, and carried both the MSS and TM sensors.

Following the failure of Landsat-6 to reach orbit in 1993, the Landsat-7 satellite was launched in 1999. Landsat-7 carried the Enhanced Thematic Mapper Plus (ETM +) sensor, which was designed to mimic the previous TM sensor while adding a panchromatic band with a 15m spatial resolution and an improved thermal detector. In April 2003, the scan line corrector on Landsat-7 malfunctioned and all Landsat-7 scenes captured since include stripes of missing data comprising approximately 22% of the total scene, a malfunction referred to as SLC-off. As a result, additional pre-processing steps, which will be discussed in a later section, are necessary when using Landsat-7 imagery. Finally, in February 2013, Landsat-8 was launched, carrying the Operational Land Imager (OLI) sensor, which mimicked the prior ETM + sensor and added two bands designed to detect coastal aerosols and cirrus clouds.⁵

The AAAS analysis focused on three types of changes: the extent and location of land reclamation; vegetation loss, particularly the loss of date palm trees; and quantification of the extent of urbanization across the island.

Landsat imagery of the study area from five different years was acquired for this analysis (Table 1). All images were captured in the summer months to ensure that detected changes were the result of actual changes in land use and land cover rather than seasonal variation. Anniversary date images were not available, but every effort was made to use images captured on as close to the same date as possible. For dates between April 2003 and the present, multiple Landsat-7 images were used to create a single composite image for the summer to correct for the scan line corrector malfunction of Landsat-7. In 2013, imagery from both Landsat-7 and Landsat-8 was acquired in order to compare the classification results of the two sensors.

Landsat Pre-Processing

Each image was visually inspected to determine if geometric correction was necessary. All of the images had been processed by the United States Geological Survey (USGS) to the Standard Terrain Correction (Level 1T). This process “(p)rovides systematic radiometric and geometric accuracy by incorporating ground control points while employing a Digital Elevation Model (DEM) for topographic accuracy.”⁶ Additional visual inspection revealed that the images were properly aligned at the pixel level. Therefore, no further geometric corrections were required.

Each image was then radiometrically corrected using a simple dark object subtraction technique. For data management purposes, Landsat data are scaled to quantized calibrated pixel values (DNs) before being distributed to the public. These DNs must then be converted to at-satellite radiance values (L_λ) using the following equation:

L_λ = DN * G_rescale + B_rescale

where G_rescale and B_rescale are the band-specific rescaling gain and bias factors provided in each image’s metadata.⁷

The at-satellite radiance values are next converted to surface reflectance values. This process adjusts for any atmospheric effects so that the reflectance values of different images can be directly compared to one another. There are many different techniques for calculating the level of atmospheric interference. The best rely on in-situ observations of the atmosphere at the time of image capture. In situations where such observations are unavailable, atmospheric effects can be estimated from dark objects. According to Chavez (1996), this process assumes that, “(w)ithin the image some pixels are in complete shadow and their radiances received at the satellite are due to atmospheric scattering (path radiance).”⁸ Chavez goes on to note that since few objects on earth are absolutely black, assuming a one-percent minimum reflectance is better than zero percent. Therefore, it is possible to calculate the path radiance (L_p) by identifying the darkest pixels in a band and using the following equation:⁹

L_p = DN_min * G_rescale + B_rescale - [(.001(E₀cos(θ₀)T_z + E_down)T_z)/π]

Here, θ₀ is the solar zenith angle provided in the metadata of each image while E₀ is the solar spectral irradiance. This can be calculated by dividing the mean exoatmospheric solar irradiance (ESUN_λ), which is provided in Landsat’s calibration parameter files (CPF), by the earth-sun distance in astronomical units (d) squared. The earth-sun distance is a function of the day of the year and the table of earth-sun distances provided by Chandler et al. (2009) was used for this analysis.¹⁰

T_z is the atmospheric transmittance along the path from the sun to the ground surface; T_v is the atmospheric transmittance along the path from the ground surface to the sensor; and E_down is downwelling spectral irradiance.¹¹ It has been shown that T_z and T_v can be assumed to equal 1.0 and E_down can be assumed to be 0.0 while still adequately correcting for atmospheric interference.¹² Therefore, the equation can be rewritten as:

L_p = DN_min * G_rescale + B_rescale - [(.001(ESUN_λcos(θ₀))/πd²]

A histogram of each band was generated to identify the DN_min value, which was set as the lowest value representing at least 1,000 pixels. The surface reflectance (ρ) can then be calculated using the following equation:¹³

ρ = [πd(L_λ - L_p)]/[(E₀cos(θ₀)T_z+E_down)T_v]

which simplifies to:

ρ = [πd(L_λ - L_p)d²]/[(ESUN_λcos(θ₀)]

The bands for all images were then stacked and clipped to include only the area of al-Awal Island and adjacent other islands as shown in Figure 1. A gapfill algorithm was then applied to the Landsat 7 images. This process, developed by the USGS, takes advantage of the fact that the gaps in images of a single area do not perfectly align. It is, therefore, possible to use pixels from one image to fill the gaps in another.¹⁴ The pixels of the fill image are calibrated to the primary image by linearly matching each image’s histogram. When the entire scene’s histogram is used, it is referred to as global linear histogram matching. However, a moving window can be used to calculate local histograms for each pixel. Local linear histogram matching is more computationally intensive, but performs better in areas with multiple types of land use and land cover. Therefore, in this study, local linear histogram matching was used, as Bahrain includes multiple land cover types. This process was performed using a free extension to Exelis ENVI 5,¹⁵ filling each image with additional images to fill in the gaps created by SLC-off.

Image Classification

Each image was classified using ERDAS Imagine’s ISODATA unsupervised classification algorithm.¹⁶ The algorithm uses multiple iterations to sort similar pixels into a user-defined number of clusters. Each cluster is then assigned by an analyst to a land cover class. If a cluster appears to contain multiple land cover types, the ISODATA algorithm can then be performed again on that cluster to further separate the pixels in a process known as “cluster busting.”

Prior to performing the classification, however, a simple density slice of the mid-infrared (band 5) band was used to separate land from water. This method relies on water’s extreme lack of reflectance in mid-infrared wavelengths and has proven to be an effective method of delineating water boundaries.¹⁷ The analyst examines the imagery to identify the threshold reflectance value that separates land pixels from water pixels. For this study, reflectance values of 0.1 were used as the threshold value. All pixels below 0.1 were assigned to the “water” class and ignored during the classification.

For the classifications, six spectral bands of each Landsat image were used as inputs for the ISODATA algorithm. Forty-eight clusters were used during the initial classification. In each classification, between four and five clusters required cluster busting. Only 12 clusters were used for cluster busting to ensure that each cluster would contain an adequate number of pixels. The clusters were then assigned to one of three land cover types: built up, soil, and vegetation. The built up class includes human built features such as roads and structures. It should be noted that human built features created from soil (e.g., new islands or dirt roads) were assigned to the “soil” class for the purpose of this study.

Validation

High-resolution satellite images from two dates were acquired to validate the classification results. As commercial high-resolution imagery was not available prior to the early 21^st century, the two years selected for validation were 2004 and 2013 (Table 2).

Imagery of two 25km² areas of the island was acquired on each date (Figure 2). The areas were selected by examining imagery from Landsat and Google Earth to identify locations that included all potential land cover types. The high-resolution images covered different geographic areas, so only one scene, in the north of the island, was common to each date.

For each scene, 500 random points were created, resulting in a total of 1,000 points for each image date. Each one of these points was manually assigned to one of four land cover classes: water, built up, soil, or vegetation. These points were then used to validate the 2004 and 2013 classifications and calculate overall accuracies, as well as user’s and producer’s accuracies.

Producer’s accuracy is the probability that a pixel of a given category was correctly classified. It is derived by dividing the number of correctly classified pixels of a class by the total number of validation points of that same class. User’s accuracy is the probability that a pixel that has been classified as being in a given class is actually of that class. It is derived by dividing the number of correctly classified pixels of a class by the total number of pixels assigned to that same class.

Including both types of accuracy is necessary, as it is possible to obtain a high accuracy for one metric at the expense of the other. For example, high producer’s accuracies can be obtained by over-classifying, i.e., including a large number of pixels not of the target class in the classification, and thereby including a high percentage of pixels of the target class. Similarly, under-classifying can yield a high user’s accuracy though it will miss a number of pixels of the target class. As a result, balancing user’s and producer’s accuracies provides the most accurate classification and offers a more robust accounting of the accuracy.

III. Results

Classification Accuracies

The classifications of the 2004 and 2013 images resulted in overall accuracies of 82% for both the 2013 Landsat-8 and 2013 Landsat-7 images. The overall accuracy of the 2013 Landsat-7 image was 80% (Tables 3-5).

In the three classifications, the user’s and producer’s accuracies for the water class are between 97% and 99%. These results demonstrate that the density slice method (described in the previous classification discussion), adequately distinguishes between land and water. However, the classifier struggled to distinguish between built up and soil pixels. All three error matrices show that a significant number of pixels of each class were misclassified. User’s and producer’s accuracies for these classes range from 57% producer’s accuracy for the built up classification of the 2013 Landsat-7 image to 88% for the producer’s accuracy of the soil classification in the same image. This classification demonstrates the problems with over- and under-classifying, as the user’s accuracies for built up and soil are 79% and 75%, respectively. Clearly, soil has been over-classified while built up has been under-classified. The classifications of the 2013 Landsat-8 and 2004 Landsat-7 images reveal similar, though less extreme, issues in accurately distinguishing between these two classes.

The classifications were more successful when identifying vegetation pixels. The 2013 Landsat-8 classification resulted in a significantly higher producer’s accuracy, implying that vegetation was over-classified. The 2013 Landsat-7 classification, however, resulted in balanced user’s and producer’s accuracies of 80% and 82% respectively. The 2004 Landsat-7 classification also resulted in acceptable user’s and producer’s accuracies.

Examining these results reveals that the classification methodology used in this study is capable of identifying water with a high degree of accuracy. Vegetation can also be adequately distinguished, though the 2013 Landsat 7 classification performed better than the 2013 Landsat-8 classification. However, there were significant issues in distinguishing soil from built up land cover types.

Island Expansion

The total extent of al-Awal Island and the islands immediately adjacent to it was calculated using the land-water mask created during the density slice of the mid-infrared band. This revealed that the total land mass grew from approximately 650 km² to 730km², or 12.5% (Figure 3).

Figure 5 shows that the total land area increased between 1987 and 1990, from 649.5 km² to 660.8 km². However, much of that increase appears to be the result of either low water levels or organic matter in the waters of Tubli Bay, which caused the water to be classified as vegetation (Figure 4).

Between 1990 and 1998 the analysis revealed a slight decrease in the total land area, from 660.8 km² to 659.1 km². However, areas classified as land in 1990 were classified as water in 1998 but new land has been added to Muharraq Island (Figures 5 and 6).

By 2005, the amount of land in the study area had increased to 680.8 km². Land reclamation continued for the next seven years and by 2013, approximately 50 km² were added. It should be noted that the 2013 Landsat-8 and 2013 Landsat-7 classifiers provided slightly different results. The total land mass was calculated as 730.8 km² for Landsat-8 and 735.0 km² for Landsat-7.

Land additions were concentrated in the northern section of Bahrain (Figure 7). The approximately 10 km² of land reclaimed between 1987 and 1998 were primarily added to al-Muharraq Island in the North East. Between 2004 and 2013, however, multiple large-scale land reclamation projects had occurred. Substantial land was added to Muharraq Island as well as the northern coast of al-Awal Island. In addition, a set of fish and crescent shaped islands was added in the south.

Vegetation Analysis

Analyzing vegetation patterns presents several challenges, as a number of factors can influence the amount of vegetation detected. The dates of image acquisition ranged from June through September, therefore, seasonal variation may have influenced the results. In addition, arid vegetation can be particularly sensitive to climatic conditions so year-to-year climate variation may also be a factor. Both seasonal and year-to-year climate variation can impact classification results and make it difficult to determine if changes in detected vegetation were the result of actual changes in vegetation rather than errors in classification.

The total amount of vegetation in the study area changed throughout the study period but did not exhibit a clear trend. Table 6 shows the total area classified as vegetation in each year. As the land area available for vegetation increased over time, the percent vegetation has also been included.

One notable feature of this analysis is the disagreement in total vegetation area between the classifications of the two 2013 images. The Landsat-8 analysis assigned 10.1 km² more land to the vegetation category than the Landsat-7 classification. This is likely due to two factors. First, it was noted above that the Landsat-8 classification had a higher producer’s accuracy than user’s accuracy. This implies that it was more likely for a non-vegetation pixel to be incorrectly placed in the vegetation class than it was for a vegetation pixel to be erroneously categorized as not vegetation. In addition, the Landsat-7 image is a composite including imagery from August, June, and September, while the Landsat-8 image was captured in July. Both of these factors may have contributed to the disagreement between the classifications.

A second potential reason for the absence of a trend in total vegetation area is that this metric does not differentiate between old and new vegetation. It is possible that vegetation was lost in one area while other vegetation was added to a different area. Figures 8 and 9 show a dramatic decrease in vegetation between 1987 and 2013 in the areas around Tubli Bay and west of Manama.

In the center and on the west side of al-Awal Island, the amount of vegetation increased between 1987 and 2013. This new vegetation is largely either agricultural or landscaped. Beginning in 1990, it is possible to observe land being converted from soil into vegetation in rectangular and circular units. The circular units of new vegetation are consistent with center pivot irrigation and the rectangular units also imply agricultural fields (Figure 10).

By 2004, several large areas of landscaped vegetation were present. Construction had begun on a golf course and at least one new compound with landscaping had been constructed. By 2013, the golf course had expanded and the road accessing the new compound had been lined with vegetation (Figure 11).

Built Up Analysis

A major challenge in remote sensing is accurately differentiating between soil and built up areas. Both types of features include a wide variety of materials with various spectral signatures. As a result, there is significant overlap between the two classes. For example, bright soil is spectrally similar to white rooftops, while dark soil is similar to roads. These difficulties are reflected in the error matrices of the 2004 and 2013 classifications, which show substantial confusion between these two classes. However, in all three classifications, the built up user’s accuracy was greater than the producer’s accuracy. This indicates that the classifications were consistently underestimating the total built up area.

Despite these limitations, the classification of built up pixels revealed a clear trend. Table 7 shows the total area classified as built up as well as percent built up. Over the 26-year study period, the area classified as built up more than doubled, from approximately 76 km² to 155 km².

The classifications show that the amount of built up area increased considerably between 1990 and 1998 and remained stable between 1998 and 2004 before again increasing dramatically between 2004 and 2013. The increase in built up area was not evenly distributed throughout the island. Much of it occurred in the area to the west of Manama (Figure 12). In addition, increases in the built up area were observed on the west side of the island. Much of this increase can also be observed in Figures 7-10. Examining the classifications (Figure 13) corroborates these observations as the density of pixels classified as built up increases in the northern section of the island throughout the study.

Discussion

Using Landsat imagery for analysis of historic land use and land cover change presents several challenges, particularly if the study area includes urban features or arid vegetation. Urban areas often include multiple land cover types in a small area. With a spatial resolution of 30m, it is almost guaranteed that a Landsat pixel in an urban area will contain several types of land cover. Similarly, vegetation in arid climates is often sparse and dispersed. As a result, background soil can dilute the spectral signiture of the vegetation. These types of mixed pixels were difficult to classify and may have reduced the accuracy of the classifications.

The type of classifier used for this study is known as a per-pixel classifier, meaning each pixel is assigned to a discrete class. Other methods, such as linear spectral mixture analysis, are sub-pixel classifiers that attempt to determine the ratio of various land cover types within a given pixel. These methods can provide a more realistic assessment of the land cover contained within each pixel but rely on the presence of adequate training data to identify pure-pixels, i.e., pixels that contain a single land cover class. These data can be derived from high resolution imagery, land cover maps, or ground surveys. For this study, high resolution imagery was only available for two of the years of interest. Historic land cover maps were also unavailable. In addition, traveling to Bahrain to identify pure pixels based on specific features that could be used as training data was not a feasible option.

Given constraints on training data availability, the ISODATA unsupervised classification algorithm was an appropriate choice dispite its limitations. However, these limitations did result in confusion in certain cases. Analysis of coastal pixels was problematic, as they often contained some water in addition to other land cover types. As a result, coastal pixels were often clustered with built up pixels. Cluster busting was used to address the issue, but the maps of built up area (Figure 13) do show areas of the coast erroneously classified as built up.

The lack of validation data for historic classifications also impeded analysis. By validating the three classifications from 2004 and 2013, it was possible to determine a range of likely accuracies. However, it is unknown exactly how accurate each individual classification was. This may have contributed to confusion regarding the classification of vegetation.

Conclusions

An unsupervised classification algorithm was used to assess land use and land cover changes in Bahrain between 1983 and 2013. Six Landsat images were acquired and classified to assess three phenomena: 1) the extent of land reclamation; 2) changes in vegetation; and 3) urbanization. The classifier was able to separate land and water with a high degree of accuracy. Approximately 80 km² of land was added between 1987 and 2013, representing a 12.5% increase in total land area. Analysis of trends in the total amount of vegetation was inconclusive. However, substantial changes in vegetation patterns were observed including the loss of vegetation in the area around Tubli Bay and the conversion of land in the west and center of the island into agriculture and landscaped vegetation. Finally, the amount of built up area on the island more than doubled from 76.2 km² to 155 km². These observations represent major changes in the land use and land cover of Bahrain that will impact the environment in as yet unknown ways.

For a PDF version of this case study, click here.

---

1. Noor Al-Nabi, Mohammad. 2012. The History of Land use and Development in Bahrain. Information Affairs Authority, Directorate of Government Printing Press. Kingdom of Bahrain.

2. United Nations Department of Economic and Social Affairs/Population Division World Population Prospects: The 2012 Revision, Volume II: Demographic Profiles. http://esa.un.org/unpd/wpp/Demographic-Profiles/pdfs/48.pdf

3. The World Bank. 31 March 2005.A Water Sector Assessment Report on the Countries of the Cooperation Council of the Arab States of the Gulf. http://siteresources.worldbank.org/INTMNAREGTOPWATRES/Overview/20577193/GCCWaterSectorReport--Englishversion.pdf

4. http://bahraniat.blogspot.com/2004/11/island-of-million-palm-trees.html

5. http://bahraniat.blogspot.com/2004/11/island-of-million-palm-trees.html

6. https://landsat.usgs.gov/Landsat_Processing_Details.php

7. Chander, G., B. Markah, and D. Helder. 2009. Summary of Current Radiometric Calibration Coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI Sensors. Remote Sensing of Environment. 113:5 pp. 893-903.

8. Chavez, P. 1996. Image-Based Atmospheric Corrections – Revisited and Improved. Photogrammetric Engineering & Remote Sensing. 62:9 pp. 1025-1026.

9. Song, C., C. Woodcock, K. C. Seto, M. P. Lenney, and S. A. Macomber. 2001. Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects? Remote Sensing of Environment. 75:2 pp. 230-244.

10. See note 7.

11. See note 8.

12. See note 9.

13. See note 9.

14. Scaramuzza, P., E. Micijevic, G. Chander. 2004. SLC Gap-Filled Products Phase One Methodology. United States Geologic Survey. http://landsat.usgs.gov/documents/SLC_Gap_Fill_Methodology.pdf

15. Excelis Code Library. http://www.exelisvis.com/Learn/CodeLibrary.aspx Gapfill.sav

16. “In an unsupervised classification, the identity of the land-cover types to be specified as classes within a scene are not generally known a priori (i.e., before the fact). Consequently, the analyst instructs the computer to use a series of statistically-based heuristic rules to group (identify) pixels with similar spectral characteristics.” This is in contrast to other classification methods where pixels of known land cover types are inputted as training data and used by the algorithm to determine the land cover of all other pixels. Jensen, J.R., J. Im, P. Hardy, and R. R. Jensen. 2009. Image Classification. In Warner, T, M. D. Nellis, and G. M. Foody eds. The Sage Handbook of Remote Sensing. Sage Publications. Chapter 19 pp. 269-281.

17. Frazier, P and K. Page. 2000. Water Body Detection and Delineation with Landsat TM Data. Photogrammetric Engineering and Remote Sensing. 66:12 pp. 1461-1467.