Focusing on the responsible application of geographic technologies across a range of human rights and humanitarian issues and how these technologies can be used to improve research and documentation.
Investigating Land Use and Land Cover Change in Bahrain: 1987-2013
Large-scale land cover and land use change is apparent between 1973-2013. Images NASA | Analysis AAAS
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.
Figure 1: Bahrain study area
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.1  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%.2  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 m3 in 1980 to 269 million m3 in 2000.3  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.4 
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.5 
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.
Table 1: Landsat imagery
|1 August 1987||Landsat 5 - TM||LT51630421987213AAA01|
|17 August 1990||Landsat 5 - TM||LT41630421990229XXX03|
|28 June 1998||Landsat 5 - TM||LT51630421998243XXX01|
|22 July 2004*||Landsat 7 - ETM +||LE71630422004204ASN01|
|6 July 2004**||Landsat 7 - ETM +||LE71630422004188ASN01|
|20 June 2004***||Landsat 7 - ETM +||LE71630422004172ASN01|
|16 August 2013*||Landsat 7 - ETM +||LE71630422013228PFS00|
|13 June 2013**||Landsat 7 - ETM +||LE71630422013164ASN00|
|1 September 2013***||Landsat 7 - ETM +||LE71630422013244ASN00|
|23 July 2013||Landsat 8 - OLI||LC81630422013204LGN00|
**Fill image- used to correct for SLC-off
***Secondary fill image- used to correct for SLC-off
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.”6  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 * Grescale + Brescale
where Grescale and Brescale are the band-specific rescaling gain and bias factors provided in each image’s metadata.7 
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).”8  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 (Lp) by identifying the darkest pixels in a band and using the following equation:9 
Lp = DNmin * Grescale + Brescale - [(.001(E0cos(θ0)Tz + Edown)Tz)/π]
Here, θ0 is the solar zenith angle provided in the metadata of each image while E0 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.10 
Tz is the atmospheric transmittance along the path from the sun to the ground surface; Tv is the atmospheric transmittance along the path from the ground surface to the sensor; and Edown is downwelling spectral irradiance.11  It has been shown that Tz and Tv can be assumed to equal 1.0 and Edown can be assumed to be 0.0 while still adequately correcting for atmospheric interference.12  Therefore, the equation can be rewritten as:
Lp = DNmin * Grescale + Brescale - [(.001(ESUNλcos(θ0))/πd2]
A histogram of each band was generated to identify the DNmin 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:13 
ρ = [πd(Lλ - Lp)]/[(E0cos(θ0)Tz+Edown)Tv]
which simplifies to:
ρ = [πd(Lλ - Lp)d2]/[(ESUNλcos(θ0)]
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.14  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,15  filling each image with additional images to fill in the gaps created by SLC-off.
Each image was classified using ERDAS Imagine’s ISODATA unsupervised classification algorithm.16  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.17  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.
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 21st century, the two years selected for validation were 2004 and 2013 (Table 2).
Table 2: High-resolution validation images
|22 July 2013||Worldview-2||103001002584C100|
|29 July 2004||Quickbird-2||1010010003229F00|
Imagery of two 25km2 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.
Figure 2: Validation image locations
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.
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).
Table 3: Error matrix for 2013 Landsat-8 classification
|Validation||Water||Built Up||Soil||Vegetation||Sum||User's Accuracy|
|Producer's Accuracy||98%||65%||84%||90%||Overall Accuracy||82%|
Table 3: Error matrix for 2013 Landsat-7 classification
|Validation||Water||Built Up||Soil||Vegetation||Sum||User's Accuracy|
|Producer's Accuracy||97%||57%||82%||88%||Overall Accuracy||80%|
Table 3: Error matrix for 2004 Landsat-7 classification
|Validation||Water||Built Up||Soil||Vegetation||Sum||User's Accuracy|
|Producer's Accuracy||98%||65%||84%||90%||Overall Accuracy||82%|
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.
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 km2 to 730km2, or 12.5% (Figure 3).
Figure 3: Total land area (km2)
Blue represent all images except for 2013 where green is the 2013 Landsat-7 image, and red is the 2013 Landsat-8 image.
Figure 5 shows that the total land area increased between 1987 and 1990, from 649.5 km2 to 660.8 km2. 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).
Figure 4: Water misclassified as land in Tubli Bay
A false color infrared display of the 1990 image (A) indicates that organic material is present in Tubli Bay (green arrows) and was subsequently classified as land during the density slice (B). Images NASA | Analysis AAAS. Coordinates 26.18 N, 50.6 E.
Between 1990 and 1998 the analysis revealed a slight decrease in the total land area, from 660.8 km2 to 659.1 km2. 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).
Figure 5: Tubli Bay water line 1990-1998
In 1990 (A) areas of Tubli Bay (green arrows) were misclassified as land due to organic material. In 1998 (B) these areas were classified as water. Images NASA | Analysis AAAS. Coordinates 26.18 N, 50.6 E.
Figure 6: Land added to Muharraq Island
Between 1990 (A) and 1998 (B), substantial amounts of land were added to Muharraq Island (green arrows). Images NASA | Analysis AAAS. Coordinates 26.25 N, 50.64 E.
Figure 7: Island expansion between 1987 and 2013
By 2005, the amount of land in the study area had increased to 680.8 km2. Land reclamation continued for the next seven years and by 2013, approximately 50 km2 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 km2 for Landsat-8 and 735.0 km2 for Landsat-7.
Land additions were concentrated in the northern section of Bahrain (Figure 7). The approximately 10 km2 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.
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.
Table 6: Area classified as vegetation
|Year||Area Vegetation (km2)||Percent Vegetation*|
*Percent vegetation based on total land area calculated for each image.
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 km2 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.
Figure 8: Vegetation loss around Tubli Bay
Between 1987 (A) and 2013 (B), much of the vegetation around Tubli Bay was lost. Images NASA | Analysis AAAS. Coordinates 26.18 N, 50.6 E.
Figure 9: Vegetation loss west of Manama
Between 1987 (A) and 2013 (B), much of the vegetation west of Manama was lost. Images NASA | Analysis AAAS. Coordinates 26.24 N, 50.53 E.
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).
Figure 10: Expansion of agriculture 1987-1998
Between 1987 (A) and 1998 (B) several rectangular fields (blue arrows) and circular fields (green arrow) were created in the center of the island. By 1998, more rectangular fields (blue arrow) and one large circular field (green arrow) were also present. Images NASA | Analysis AAAS. Coordinates 26.14 N, 50.51 E.
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).
Figure 11: Landscaped vegetation added to center of the island
By 2004 (A) a golf course had been constructed (green arrow) as well as a landscaped compound (blue arrow). By 2013 (B) the golf course had been expanded (green arrow) and vegetation had been planted along the road accessing the compound (blue arrow). Images NASA | Analysis AAAS. Coordinates 26.14 N, 50.51 E.
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 km2 to 155 km2.
Table 7: Area classified as built up
|Year||Area Built Up (km2)||Percent Built Up*|
*Percent built up based on total land area calculated for each image.
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.
Figure 12: Urbanization west of Manama 1990-2013
Large increases in the urban extent were observed between 1990 (A) and 1998 (B). Major increases were again observed between 2004 (C) and 2013 (D). Images NASA | Analysis AAAS. Coordinates 26.19 N, 50.50 E.
Figure 13: Built up area 1987-2013
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.
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 km2 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 km2 to 155 km2. 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 
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 
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.