Detecting, extracting, and mapping of inland surface water using Landsat 8 Operational Land Imager: A case study of Pune district, India

Introduction

By 2030, the world is anticipated to face a 40 percent water crisis as a result of urbanization, strain on food production, population increase, and increasing demand from the industrial, domestic, and energy sectors^1,2. The current water crises are a result of poor water management and not because of adequacy of water. Sustainability of water can be maintained with proper governance and appropriate water management practices³. Inland water bodies are critical for the survival of all forms of life and need to be monitored frequently⁴. The regular monitoring and mapping of surface water is of paramount importance for better water resource management⁵.

Surface water bodies are dynamic in nature as they shrink and expand with the course of time⁶. Remote sensing is an effective and low cost alternative to conventional in-situ water quality monitoring systems^7,8. Recent advancements in optical remote sensing have catapulted the field of surface water sensing into a new era. It has been used in water resource assessment and management in the last four decades^9–19. These applications involve delineation of surface water and assessing water quality using thematic information extraction techniques²⁰.

Both active and passive remote sensing have been used to detect and extract surface water bodies⁹. Optical remote sensing of water is based on the difference between spectral reflectance of land and water. Water absorbs most of the energy beyond near-infrared wavelength whereas other land features have higher reflectance in those wavelengths. Thus, water appears as a dark spot in near-infrared and higher wavelength bands of multispectral images which makes it easy to differentiate water from non-water topographic features. The accuracy depends on the weather conditions during the image captured and topography of the study area. Microwave is very effective during the presence of clouds, rain fields, the presence of water vapor and aerosols. Microwave remote sensing of surface water is based on the backscattering properties of water in the microwaves region. The surface of water behaves as a specular reflector for incident microwaves, and responds with no or very low signal returns. Thus, water is easily distinguished from other topographic features in and around the area.

Applications of remote sensing for various purposes of water resource management were proposed for the first time for surface water bodies in India²¹. Moderate Resolution Imaging Spectroradiometer (MODIS), a coarse resolution but high temporal resolution sensor have been successfully applied to near-real time surface water monitoring^22–25. The Visible Infrared Imaging Radiometer Suite onboard Suomi National Polar-orbiting Partnership (Suomi NPP-VIIRS) has been tested for detecting surface water²⁶. The low resolution satellites have major limitations in observing small water bodies. Cloud cover is another issue often associated with low resolution satellites. Landsat is the most successful series of satellites in monitoring surface water bodies in the history of remote sensing⁹. Landsat 8 Operational Land Imager (OLI) is the most recent satellite in the Landsat series and is widely used in detecting surface water bodies^11–19. Systeme Probatoire d’Observation dela Tarre (SPOT), which has an improved spatial resolution of 10 m but that is not available free of cost, is used in a few studies for applications of surface water^27,28. The MultiSpectral Instrument (MSI) onboard Sentinel-2 with spatial resolutions of 10m and a revisit time of 10 days at equator and five days with two satellites under cloud free conditions is used to detect surface water^29–31. Many works specifically reviewed recent accomplishments, current status, and limitations in this domain, including geographic scale, integration with in-situ hydrological data, cloud obscuration, and other land feature classes.^{7–9,32–34}.

Water extraction methods using remote sensing have been pursued by advancing sensors and through the development of various methods like thematic classifications, Linear Unmixing, single band thresholding and multispectral water indices^5,9,20,35. The method selection from the list above is based on the scale of the water map and the proportion of water bodies to non-water area in the image.

The method of extracting features of surface water using water indices remains popular because it is simple to implement, has a relatively fast processing time, and produces fairly accurate results^36,37. These methods have produced excellent results when used with Landsat imagery^11,38,39. Surface water features can be distinguished using a variety of water indexes, including Normalized Difference Water Index (NDWI), Modified Normalized Difference Water Index (MNDWI), and Automated Water Extraction Index (AWEI). NDWI uses the green and near infrared (NIR) band to delineate water features in the image⁴⁰. Nevertheless, The NDWI is unable to suppress the built-up land signal which make it difficult to differentiate water and built up area effectively²⁰. MNDWI is a modified version of NDWI that uses green and short-wave infrared (SWIR) bands to overcome limitations encountered by NDWI for extracting water using space observation. AWEInsh, an automated water extraction index with no shadow, and AWEIsh, an automated water extraction index with a shadow, are two more additional water indices that have been proposed for water delineation⁴¹.

Selecting an optimum threshold to discriminate water accurately from land is a challenging and complex task⁴¹. Ideally, spectral index images assign higher values to water and negative values to the non-water pixels. Thus, zero is ideally a reference value to separate water and non-water pixels. However, in practice most of the times the threshold values are non-zero. In this study automated Otsu’s threshold method and trial-and error methods are used to estimate optimum threshold values^42,43.

Pekel (2016) presented a comprehensive validation protocol for accuracy analysis of the Global Surface Water Explorer (GSWE) dataset⁴⁴. This dataset has been validated using the total 40,124 control points distributed geographically and temporarily all across the globe. GSWE is produced from Landsat imagery at the global scale over the past 32 years.

In this study, we demonstrated a case study of Pune district for detection and mapping of surface water using the MNDWI spectral index method. The objective of the study is evaluation of the MNDWI water indices method for extraction of surface water from Landsat 8 OLI images having multiple water bodies located spatially apart. The time series surface water map is prepared for the entire year using the cloud free images. The surface water maps so developed are assessed for accuracy by comparing them with a reference surface water map prepared by The European Commission’s Joint Research Centre in the framework of the Copernicus Programme.

Methods

The processes followed for achieving the objectives of study are represented in Figure 2. In summary, the process starts with downloading the cloud free Landsat 8 OLI of the study area (see Underlying data) and carrying out initially radiometric and geometric corrections, and cloud detection. In the next step MNDWI spectral indices were calculated from the Landsat images in order to enhance all water pixels in the image. Water bodies were extracted using the threshold slicing method using Otsu and trial and error method. Finally, the images were assessed for accuracy by comparing with a reference image. All processes followed are explained in the separate subsections below.

Study area

The study area for development of surface water maps is between 18.01 N and 19.36 N and 73.34 E and 74.76 E which is a portion of the Pune district administrative boundary in India. Pune lies on the west of the Deccan plateau (Figure 1), at an altitude of 560 m (1,840 ft) above mean sea level. It is located in western Maharashtra along the foothills of the Sahyadri Mountains⁴⁵.

Figure 1. Study area location and extent.

Landsat 8 image courtesy of the U.S. Geological Survey.

Figure 2. Flow chart of the methodology used.

OLI, Operational Land Imager; ToA, Top of Atmosphere; MNDWI, Modified Normalized Difference Water Index; JRC, Joint Research Centre; GSWE, Global Surface Water Explorer. Landsat 8 images courtesy of the U.S. Geological Survey.

This area is chosen for two reasons. First, it has heterogeneous topography comprising hills on the western side and plain terrain on the eastern side of the district. It covers all the types of surface water bodies, such as rivers, lakes, and reservoirs. Secondly, the authors are from this region, which provides prior knowledge and understanding of the selected area.

According to the 2011 Census of India, Pune district had a total population of 9.429 million, meaning it lies within India’s top ten most populated cities⁴⁵. The monsoon season in Pune is observed from June to September, and the cloud presence makes it challenging for water extraction from remote sensing. Therefore, surface water extraction analysis were performed for the non- rainy season.

Image pre-processing

The first step in the pre-processing is to transform raw digital numbers (DN) to Top of Atmosphere (ToA) reflectance (Figure 3). The pre-processing is performed using Semi-Automatic Classification Plugin of QGIS (3.22.3) open source software (please note that all pre-processing of the raw images, calculation of MNDWI, cropping of images and thresholding using trial and error presented in the methods is performed using QGIS open source software). Landsat Collections Level-1 data can be re-scaled to top of atmosphere (ToA) reflectance using radiometric re-scaling coefficients (see Underlying data⁴⁶). DN is a mathematical value of spectral data and not physical measured data so it is to be converted to ToA. DN can be converted to ToA reflectance using the re-scaling coefficients in the MTL using (Equation 1):

$${\rho }^{\prime }=M_p\ast Q_{cal}+A_p(1)$$

Figure 3. Raw Landsat 8 image of 6th February 2016 representing conversion of Digital Numbers (DN) to Top of Atmosphere (ToA) Reflectance.

Landsat 8 images courtesy of the U.S. Geological Survey.

Where M_p and A_p are the metadata defined multiplicative and additive band specific re-scaling factors, respectively. The quantized and calibrated standard product pixel values are denoted by Q_cal. It is not corrected for the sun angle, which results in the ToA reflectance as (ρλ),

$$\rho \lambda =\frac{\rho {\lambda }^{\prime }}{cos{\theta }_{SZ}}=\frac{\rho {\lambda }^{\prime }}{sin{\theta }_{SE}}(2)$$

where θ_SE is the local sun elevation angle as specified in the metadata and θ_SZ denotes the local solar zenith angle as calculated using (90⁰ -θ_SE).

MNDWI spectral index

MNDWI were initially calculated to enhance the difference between water and non-water areas. Surface water bodies tend to have positive values in MNDWI images, whereas, soil, vegetation and built up areas are expected to represent negative values²⁰. Xu (2006) proposed that the Short-wave Infrared (SWIR) band can reflect subtle features of water more effectively compared to Near Infrared (NIR) band and substituted SWIR for NIR band in NDWI²⁰. The MNDWI method has been commonly used and is a powerful index that can extract water bodies. It is expressed by an Equation (3).

$$mNDWI=\frac{\rho 3-\rho 7}{\rho 3+\rho 7}(3)$$

Where the subscript ρ denotes the surface reflectance value computed from bands 3 (green), and 7 (visible) (SWIR-2).

This study restricts the surface water bodies mapping to the Pune district administrative area. The MNDWI images are subsetted using the shapefile to remain focused on the surface water bodies in the study area (Figure 4). All the image processing algorithms mentioned in this section are performed using QGIS (3.22.3) open source software.

Figure 4. Time series Modified Normalized Difference Water Index (MNDWI) derived from Landsat 8 Operational Land Imager (OLI).

Landsat 8 image courtesy of the U.S. Geological Survey.

Trial and error method of threshold

The optimal threshold estimated using the Otsu threshold method did not work satisfactorily and the image obtained encountered ‘salt-pepper error’ (Figure 5). The alternative, iterative trial and error method based on visual interpretation was used in this study to find the optimal threshold to discriminate water and non-water from the MNDWI image. An initial threshold of 0 was applied to all the water index images to judge classification status by visual interpretation. Based on the initial screening, the threshold is adjusted to get higher kappa values. The optimal threshold so selected was used to generate surface water map. This water map was analyzed for classification accuracy compared with the JRC monthly map, available in Underlying data⁴⁸.

Figure 5. Result of Otsu’s threshold for Modified Normalized Difference Water Index (MNDWI) image on 6th February 2016.

Landsat 8 image courtesy of the U.S. Geological Survey.

Results and discussion

After the radiometric corrections and the creation of the subset to study area extent, all images are processed (using QGIS 3.22.3) for derivation of water indices (see Figure 4). The figures use the YYYYMMDD format to express dates corresponding to the Landsat image. The MNDWI images clearly visualizes enhancement in open surface water features. Surface water features appear bright compared to other non-water features. The western side of the study region is mountainous, which is where the majority of surface water bodies are located. Among all the derived images the MNDWI index is in the range of -1 to 0.86. The brightness variation is observed because of the topographic corrections characterized by undulating and complex terrain. MNDWI images derived on 22nd February 2016 and 9th March 2016 indicate cloud cover of 23.09 % and 9.83 %, respectively.

The water indices images are processed for separation of water and non-water regions using a threshold slicing method. The Otsu’s threshold method is evaluated for threshold slicing for all available images. The derived binary classified image, consisting of water and non-water pixels, is compared for accuracy to the JRC GSWE water map for the corresponding month using a confusion matrix. The threshold value and confusion matrix parameters for each image are presented in Table 2. The result obtained using the Otsu’s threshold method were unsatisfactory and the model failed to produce decent results (see Figure 5). The sample calculation for clear and cloud free image captured on 6th February 2016 have been presented for discussion. The results for the rest of the images are similar. The input MNDWI image, corresponding grey scale histogram and output surface water image from Otsu’s threshold are presented in the Figure 5. MNDWI images have the range of values(-1 to 0.86), resulting in the formation of multi-modal histograms. The output image is a result of ’salt-and-pepper’ noise. The threshold value obtained by Otsu method for this image is -0.04. The threshold divides images into two classes: pixel values less than -0.04 and pixel values greater than -0.04.

Forward class pixels with values less than the cutoff have a majority of non-water characteristics. This portion of the binary image appears as ’salt’. The pixels with values more than threshold includes maximum pixels with water and other features as well. Background classes have a large number of pixels compared to forward classes. However, it has been demonstrated that this method produces unstable results when the proportions of water and non-water pixels are not equal and do not form a bimodal histogram³⁴. Data transmission error while calculating the threshold in Otsu’s algorithm could be another reason for the unsatisfactory result.

The iterative trial and error method was used as an alternative to the Otsu method (Figure 7). The optimum threshold value was determined for each image. The result classified binary image is compared with the reference image. The confusion matrix parameters are performed for all images of MNDWI derived by both Otsu’s Method and trial and error method (Table 2). The threshold by the Otsu’s method is less than zero, whereas, the threshold by trial and error method are all positive values. Based on the classification point of view, the producer’s accuracy (PA) represents error of omission. Otsu’s method analyses the gray scale histogram and classifies more pixels as water. PA of surface water map derived using the Otsu’s method and trial and error method, indicate the proportion of the GSWE image that is classified as water. More pixels of the water class reduces the percentage of omission in comparison with the GSWE water map. Both methods result in high PA for classification. UA of MNDWI surface water map derived using the Otsu’s method and trial and error method, represent the probability that pixels classified as water actually represent water in the GSWE data. UA measures error of commission, which represents the number of pixels assigned to the incorrect class. The Otsu’s threshold method classified all pixels with a gray scale level more than the respective threshold as water. This class has other topographic features along with water. Thus, the error of commission obtained is high for the MNDWI surface water maps derived by the Otsu’s method. However, iterative trial-and-error method performed better in terms of UA for surface water mapping. The minimum threshold value 0.12 is applied to the image captured in May and maximum threshold value 0.46 is applied to images acquired in February. The threshold values depend on the MNDWI range of the image. The MNDWI index observed is higher in February and lower in May among all the images of the study area. The summer season in the study area peaks around May during which the majority of surface water bodies shrink and have less water availability (see Table 3). The reduced amount of water causes turbidity and exposes silt affecting the satellite reflectance, resulting in a reduction in the number of classified water bodies. This is the reason for less MNDWI index during the summer season, between April and May. The monsoon season is observed during June to September in the study area during which the reservoirs are in full water condition (i.e. the reservoir is filled with maximum capacity of water). However the presence of cloud limits the availability of clear images. The images of immediate post monsoon season are observed with high (>0.7) MNDWI values. The non-water area has a wetness factor, which is why the lower MNDWI index has been raised from -1 to -0.7. The reservoirs with full water conditions can be easily visualized and delineated by the MNDWI method.

Mountain shadow

Mountain shadows have sometimes been misclassified as water. The DEM have been used to detect the misclassified portion of the surface water by mountain shadow noise (Figure 6 (a)). The mountain shadow depth in relation to sun elevation angle and DEM have been derived, and the mountain shadow depth can be seen as the white spots in the Figure 6 (b). The surface water derived by trial and error method to MNDWI spectral image is the overlay with the mountain shadow depth in Figure 6 (c). This misclassification is removed by applying logical operation (Figure 6 (d)). Surface water bodies located towards the western hilly regions of the study area are extracted successfully (Figure 8). The MNDWI result is presented in Table 2, and all are followed by shadow depth analysis for removal of misclassification of water bodies.

Figure 6.

(a) Digital Elevation Model (DEM) of Pune. (b) Mountain shadow depth derived from DEM. (c) Misclassified mountain shadows as water (red area is water and white area is mountain shadow). (d) Integrating shadow depth with binary Modified Normalized Difference Water Index (MNDWI) image to remove misclassified water. Landsat 8 image courtesy of the U.S. Geological Survey.

Figure 7. Sample image of trial and error threshold images on the date 21st January 2016 and 6th February 2016.

Landsat 8 image courtesy of the U.S. Geological Survey.

Annual surface water occurrence map

The intra-annual distribution of surface water is prepared by mapping maximum water extents during eight months of the year in 2016. The comprehensive picture of the dynamics of the surface water on yearly basis is prepared. The percentage of occurrence at which water was present on the surface during all months except monsoon season, June to September, was captured in a single image (Figure 8). Performing the accuracy analysis of the map is not possible and can be compared with a similar database. The JRC GSWE map has weekly, monthly and yearly data of surface water. The annual surface water map, JRC GSWE annual image of 2016, is delineated for the study area extent (Figure 9). The JRC GSWE yearly image is a single image representing all months of water occurrence in that year, in terms of permanent and seasonal water.

The overall water occurrence computed by MNDWI during the dry period only, would be biased because it would give more weight to the dry period than the wet period.

This could be a limitation of the study and proposed for a future aspect of the research.

Figure 8. Annual water occurrence by Modified Normalized Difference Water Index (MNDWI) method and Joint Research Centre Global Surface Water Explorer.

Figure 9. Identified reservoirs location and shape.

Reservoir surface area comparison

The surface area of seven clearly visible reservoirs in the MNDWI maps and with surface area more than 500 hectors, were estimated for each of the eight months and compared to the surface area of those reservoirs calculated using GSWE data for the same months (Table 3). The seven reservoirs are Bhama Askhed, Bhatghar, Chaskaman, Mulashi, Panshet, Shivrata and Varasgaon. The information on water occurrence and recurrence is presented in Figure 9 and Figure 10. The pattern in changing surface area by both the methods is the same. The surface area of all reservoirs is higher by GSWE methods during all the eight months except October and November, which is immediately after the monsoon months. The maximum difference in the values of surface area is observed in May and October. During May, the water storage is less in reservoirs with deteriorated water quality. In both the cases, the MNDWI method classified water for less surface area values. After the immediate monsoon, the surrounding dry land areas of the water bodies are generally misclassified as water bodies class areas due to wetness during the monsoon the in spectral indices method. The thematic maps and temporal profile were used to characterize the transition between the months of the years. MNDWI performed better in delineating large water bodies, as demonstrated by comparison to the reference data. The method also monitored the changes in quantity of water on a yearly basis.

Figure 10. Comparative representation of surface area of distinct reservoirs in study area calculated by Modified Normalized Difference Water Index (MNDWI) and Global Surface Water Explorer (GSWE) methods.

Conclusions

Inland surface water bodies in Pune district were detected and mapped using MNDWI spectral index derived from Landsat 8 OLI. The surface water was detected and mapped using cloud free images from eight months of the year 2016. Long term monitoring of surface water using spectral indices derived from Landsat is possible, however it is susceptible to cloud presence and some topographic restrictions.

Automatic Otsu’s threshold method is tested to separate water and non-water classes from the MNDWI image. The results obtained with Otsu’s method are not satisfactory. The MNDWI images input to Otsu’s method represent a threshold in the range -1 to 0.86. The image is multi-modal and has more than two classes. Otsu’s method requires that the number of pixels of both the classes should be proportionate. The MNDWI image has large discrepancy between the number of water and non-water pixels. The trial and error method of threshold slicing is tedious and time consuming. The trial and error method was used in iterative manner to achieve a better result. The threshold is different for all MNDWI images acquired on different dates. The Kappa coefficient represent better agreement between MNDWI and GSWE reference maps. The result obtained by trial and error methods are satisfactory.

The shadow depth derived using DEM is integrated with the MNDWI image to remove the shadow effect. The shadow effect is successfully removed and better results are achieved.

An annual surface water occurrence map is prepared from monthly surface water maps to study the dynamics of surface water. Permanent water occurrence at seven reservoirs with more than 100 Ha surface area were analyzed and compared for surface water areas on a yearly basis with reference images.

Accuracy analysis of the result is challenging because of the unavailability of physically measured data. The physically measured data for such cases is rarely available. The result obtained in this study are compared with the reference map. Hence, the comprehensive accuracy analysis could be a limitation of this study. The performance of the method may differ from the result presented in this study, by implementation of improved automatic threshold method, consideration of vegetation cover, and water quality as an aspect of further research.

Data availability