Subpixel Mapping of Surface Water in the Tibetan Plateau with MODIS Data

: This article presents a comprehensive subpixel water mapping algorithm to automatically produce routinely open water fraction maps in the Tibetan Plateau (TP) with the Moderate Resolution Imaging Spectroradiometer (MODIS). A multi-index threshold endmember extraction method was applied to select the endmembers from MODIS images. To incorporate endmember variability, an endmember selection strategy, called the combined use of typical and neighboring endmembers, was adopted in multiple endmember spectral mixture analysis (MESMA), which can assure a robust subpixel water fractions estimation. The accuracy of the algorithm was assessed at both the local scale and regional scale. At the local scale, a comparison using the eight pairs of MODIS / Landsat 8 Operational Land Imager (OLI) water maps demonstrated that subpixels water fractions were well retrieved with a root mean square error (RMSE) of 7.86% and determination coe ﬃ cient (R 2 ) of 0.98. At the regional scale, the MODIS water fraction map in October 2014 matches well with the TP lake data set and the Global Lake and Wetland Database (GLWD) in both latitudinal and longitudinal distribution. The lake area estimation is more consistent with the reference TP lake data set (di ﬀ erence of − 3.15%) than the MODIS Land Water Mask (MOD44W) (di ﬀ erence of − 6.39%).


Introduction
The Tibetan Plateau (TP) is known as the earth's third pole [1,2] and the water tower of Asia, and is an important factor greatly affecting the earth's climate system [3,4]. As the headwater area of ten major rivers in Asia, the TP supplies water resources to the billions of inhabitants in the plateau and surrounding regions [5]. Surface water such as lakes, rivers, and artificial reservoirs in the TP, which have been greatly affected by climate change [6] and human disturbance [7,8], are playing an essential role in climate balance, hydrological cycle, and ecosystem balance [9,10]. Therefore, long-term knowledge of the state of water bodies on the TP can improve understanding of the impact of climate change on the cryosphere and hydrosphere.
The algorithms and techniques proposed for mapping water bodies can be categorized into two major types: (1) manual compilation method based on remote sensing images or geographic information system (GIS) database, and (2) automated mapping method based on remote sensing images.
There are several global and regional surface water data sets developed with the compilation method. The Global Lake and Wetland Database (GLWD) was compiled [11] based on seven GIS databases gathered from the 1990s. Wan et al. [12] published a data set of lakes over the TP for three periods: 1960s, 2005, and 2014. The manual interpretation and digitization approaches were applied to the ground survey data and the high spatial resolution satellite data mainly from the China-Brazil Earth Resources Satellite (CBERS) and the GaoFen-1 satellite (GF-1). For the compilation method, digitization through visualization with prior knowledge can ensure relatively high accuracy of water boundary extraction and provide databases with strict quality control [13,14]. However, the use of these data sets is constrained by their static characteristics, which require higher temporal resolution water maps to reveal long-term changes of water bodies.
Many water mapping algorithms have been developed to extract and delineate open water features by utilizing various satellite images with different spatial, temporal, and spectral resolutions, among which high spatial resolution sensors such as Landsat satellites are the most widely used. The techniques mainly include the thresholding method [15][16][17][18][19], spectral water index method [20,21], linear spectral unmixing model [22], and data exploration techniques [23]. However, due to the long revisit intervals, the main limitation of large-scale use of Landsat data is the availability of cloudless ground observations [24]. For example, multi-year time series Landsat images are required to compile global cloudless images [15].
Moderate Resolution Imaging Spectroradiometer (MODIS) data with daily repeat time and medium spatial resolution has the advantage of monitoring detailed water area variation and has been applied to detect seasonal and interannual changes of water bodies regionally and globally [25][26][27]. Using the MODIS 8-day reflectance product, the support vector machine classification [28], and a threshold segmentation algorithm [29] were applied to map the spatial distribution of large lakes over the TP. However, the low spatial resolution of MODIS data will affect the purity of pixel and the clarity of water boundary, thus limiting the accuracy of the extracted feature [30].
Measuring long-term TP surface water changes with both high temporal resolution (~10 days) and high accuracy remains a challenge. Multiple endmember spectral mixture analysis (MESMA) is commonly used to map subpixel land cover distribution from coarse resolution images [31][32][33][34][35]. However, because environmental factors on the TP (such as clouds, snow, glaciers, and shadows) complicate the classification process, there are few or no examples of applying MESMA to generate water fraction maps in an automated and operational manner.
In this paper, we present a comprehensive method that is well-suited for mountain environments to map water bodies over the TP. A multi-index threshold endmember extraction method was applied to automatically select the endmembers from MODIS images. To incorporate endmember variability, a combination called typical and adjacent endmembers is used in MESMA endmember selection strategy. Section 4 details the accuracy assessment of the algorithm at both the local scale and whole regional scale. The cross-comparison was performed with the Landsat 8 Operational Land Imager (OLI) images and four public data sets, namely, GLWD, MODIS Land Water Mask (MOD44W), the TP lake data set, and the newly published European Space Agency (ESA) surface water data. Section 5 discuss the advantages, limitations, and implementation of the water mapping algorithm.

Area of Interest
Using the National Aeronautics and Space Administration (NASA) Space Shuttle Radar Terrain Mission (SRTM) digital elevation model (DEM) database v4.1, the boundary of the TP (Figure 1) is defined as 2500 m a.s.l (above sea level) with a total area of 3.08 × 10 6 km 2 . There are more than 1200 lakes with an area greater than 1 km 2 in the Tibetan Plateau, accounting for half of the total number and area of lakes in China [36]. Lakes on the TP are rarely affected by human activities and provide important indicators for studying climate change. Many major Asian rivers such as the Indus, Ganges, Brahmaputra, Yangtze, and Yellow rivers, the sources of which are all located on the TP, provide water for more than 1.4 billion people [37]. For comparison and analysis, following the definition of basin defined by Wan et al. [12], the TP is divided into 12 basins (Figure 1). These basins include nine exorheic drainage basins (i.e., AmuDarya, Brahmaputra, Ganges, Hexi, Indus, Mekong, Salween, Yangtze, and Yellow) and three endorheic drainage basins (i.e., Inner TP, Tarim, and Qaidam).

Elevation Data
The SRTM DEM data [38] is widely used to derive the boundary of watersheds and basins, and to delineate the shorelines of lakes, rivers, and oceans [39,40]. The gap-filled seamless SRTM DEM Version 4.1 (90 m, accessed at [41]) was used to compute shaded relief maps in ENVI v5.3. The threshold that shaded relief less than 0.1 was applied to produce the shadow map.

MODIS MOD09A1 Product
The MODIS MOD09A1 Version 6 product was downloaded from the NASA's Earth Observing System Data and Information System (EOSDIS) [42]. MOD09A1 provides an 8-day composite 7-band reflectance product at 500 m resolution that is atmospherically-corrected along with the associated quality assurance (QA) information. By selecting the high observation range, low viewing angle, absence of cloud or cloud shadow coverage and low aerosol loading, each MOD09A1 pixel contains the best possible observation during an 8-day period [43].

Validation Data
The accuracy of the algorithm was assessed at both the local scale and regional scale. At small scale, we used the Landsat 8 OLI images as reference data to validate the MODIS water unmixing method. Over the basin scale, the MODIS water fraction maps in 2014 (8 October 2014) are compared with four public data sets: the TP lake data set, MOD44W, ESA surface water data, and GLWD.

Landsat 8 OLI
The Landsat-8 OLI product provides reflectance at 9 bands with a spatial resolution of 30 m, including a 15-m panchromatic band with an imaging width of 185 km × 185 km. The C1 Level-1 OLI data were downloaded from the US Geological Survey (USGS) website [44], and are in the World Geodetic System (WGS84) datum and the Universal Transverse Mercator (UTM) projection. The calibration process converts the digital numbers value to radiance, then uses the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) algorithm for atmospheric correction to obtain OLI surface spectral reflectance. To assess the algorithm accuracy at the local scale, eight pairs of OLI/MOIDS images acquired within eight days (see Table 1) were used (Section 3.1). Since the co-registration between OLI and MODIS data was found to be unsatisfactory through overlaying and visual comparison, the OLI water map was converted to Sinusoid projection and carefully registered to the MODIS water map with a root mean square error (RMSE) of less than 10 m. The cloud and associated shadow areas in the OLI image were detected using the OLI band quality assessment (BQA) ancillary data [45] and excluded from comparison. The data set of lakes over the TP [12] was downloaded from [46] https://figshare.com/articles/ Data_TPLakes/3145369. The data set includes three time series: 1960s, 2005, and 2014. Through manual interpretation and digitization of the high spatial resolution satellite images (such as CBERS and GaoFen-1), the boundaries of lakes (area greater than 1 km 2 ) were determined using the median lines of the wet (August-October) and dry (April and May) season boundaries.

MOD44W
The MODIS Land Water Mask (MOD44W) is a global water mask at 250 m resolution produced annually from 2000 to 2015 [47]. For collection 6, a decision tree classifier is applied to MODIS 16-day composite image to determine the presence of water, which is accumulated at each pixel to produce an annual probability summary layer. Pixels with a probability greater than 50% are classified as water in MOD44W. The Version 6 data have a water mask layer and a quality assurance (QA) layer stored in a Hierarchical Data Format 4 (HDF4) file.

GLWD
The GLWD [11] was downloaded from [48]. The data set was compiled based on various forms of digitized archives, including GIS database and satellite images acquired during the 1990s. The GLWD comprises large lakes and reservoirs (level 1), smaller lakes, reservoirs and rivers (level 2), and other wetland types (included in level 3) such as marsh, peatland, etc. In this study, level 1 and level 2 data were combined for comparison.

The ESA Surface Water Data
The ESA surface water data [23] were released by the European Commission, Joint Research Centre. Three million Landsat satellite images (Landsat 5, 7, and 8) were used and many techniques for big data exploration and information extraction were exploited to quantify the changes in global open water from 1984 to 2018 at 30 m resolution. The data are available both in Google Earth Engine (GEE) and the website [49]. The image collection "MonthlyHistory" (October 2014) in GEE was used in this study. Since it is hard to separate lakes from other water bodies such as rivers automatically, buffer zones of the TP lake data set (buffer distance 2 km) were created to extract the lakes in ESA data. With the support of cloud-based GEE computing platform, the water areas divided by basins were calculated.

Method
The following sub sections detail the entire process from the MODIS MOD09A1 product to the maps of subpixel water fraction, which are illustrated in Figure 2. Firstly, a multi-index method is applied to the MODIS-derived indexes to extract endmembers, including water, vegetation, barren, and snow endmembers. Based on an assumption that the mixed pixels which contain water bodies only exist in the neighborhood of water endmembers (pure water pixels), the water mixed pixels are identified using a morphological dilation method. Then MESMA is applied to the possible mixed water pixels to generate subpixel water fraction map. To further clean the effect of mountain shadow, the shaded relief algorithm is implemented to the SRTM DEM to generate a mountain shadow map, which is finally masked out from the original water fraction maps.

Endmember Selection
Appropriate endmember selection method is the key to successful spectral mixture analysis [50]. There have been a number of researches on selecting the optimal number and type of endmembers for a certain satellite image [51][52][53][54]. In this study, we apply a multi-index threshold method [55][56][57] to select the endmembers, which are classified as vegetation, snow, barren, and water. Based on three satellite-derived indexes listed in Table 2, the rules of identifying endmembers are shown in Table 3.
The water endmembers must satisfy two constraints simultaneously: (1) Normalized Difference Water Index (NDWI) values are greater than 0.1 and (2) reflectance in MODIS band 2 is less than 0.2. These criteria must be carefully examined to make sure the selection of only pure water pixels and minimum omission errors. Since the criteria can ensure the endmembers within each endmember class have a wide range of spectral signatures, our endmember selection method can address the spectral variability which can be caused by numerous factors including illumination conditions, atmospheric effects, vegetation phenology, etc. Table 2. Satellite-derived indexes used for endmember selection (Bi is the ith band of MOD09A1).

Index Equation Remark
Reference Water has positive value [20] Normalized Difference Snow has a greater value [58] Normalized Difference Vegetation has a greater value [59]

Identification of the Mixed Pixels Contain Water Bodies
Since the water and shadow share a similar spectral signature of MODIS band, it is difficult to distinguish between them, especially for mixed pixels (e.g., the pixels which contain water and bare land and the pixels which contain shadow and land are indistinguishable in MODIS spectral). Therefore, it is a better choice that the unmixing process is applied to only the mixed pixels which contain water bodies rather than to all the mixed pixels. In this study, we assume that the possible water mixed pixels only exist surrounding the pure water body (e.g., the lake shoreline). This assumption can avoid many spurious pixels which contain shadow rather than water bodies especially in mountainous regions. Based on the water endmembers (pure water pixels) extracted in Section 3.1, a basic operation of mathematical morphology was implemented to identify the mixed water pixels, which were extracted through the morphological dilation of the water endmember pixels by a 3 × 3 kernel. Due to the coarse resolution of the MODIS, our algorithm only accounts for water bodies with a width of more than 500 m; therefore, small isolated water bodies with a width less than 500 m are not considered in this study.

Spectral Mixture Analysis
MESMA is a technique based on the assumption that the reflectance of each pixel can be expressed as a linear combination of reflectance of discrete endmembers, as shown in Equation (1). To achieve physically meaningful endmember fraction estimates, this study utilizes the sum-to-one and non-negative constraints (Equations (2) and (3)). MESMA analyzes the individual linear spectral mixture for each permutation of two or more endmembers, so that selected endmembers can change in number and type for the target mixed pixel. Since there are numerous endmember candidates selected from the MOD09A1 images, to reduce the computing time, it is assumed that only two endmember types mainly account for the contribution of the reflectance ρ i in Equation (1). All possible combinations of water endmembers and nonwater endmembers are tried by an iterative procedure to find the best fit model (lowest RMSE in Equation (4)) for each mixed pixel.
where ρ i is the reflectance of the mixed pixel in band i, r ij is the reflectance of endmember j in band i, f j is the fractional cover of endmember j, and ε i is the residue error for band i.

Approaches of Using Endmembers in Subpixel Water Mapping
Optimally, the ideal choice for spectral mixture analysis would be the minimum number of endmembers describing the spectral variability of the image scene [60]. Too many endmembers that are spectrally similar will increase the computation cost and may cause the endmember fractional cover to be physically inaccurate based on field assessments [53,61]. However, our endmember extraction method selects a great number of endmembers for a certain class type (e.g., water and vegetation), since each MODIS tile covers an area of about 1200 km × 1200 km and contains 2400 × 2400 pixels. In addition, the theory of spatial dependence believes that by sharing the same light and environmental conditions spatially-adjacent pure pixels of the same land cover class are more likely to have the same spectral reflectance than distant ones [22,62].
To reduce this within-class variation, we adopted an endmember selection method [55][56][57] which selects the most representative endmember from both local scale and global image. The selection strategy is called the combined use of typical and neighboring endmembers. The typical endmembers refer to the global endmembers in an image which share the most common spectral signature. The spectrum of the typical endmember is calculated as the mean of that of all the endmembers within a land cover class. The neighboring endmembers are the local endmembers within a 9 × 9 grid (window size 4.5 km × 4.5 km) of which the mixed pixel is located at the center. The spectra of the neighboring endmembers were used to consider the spectral similarity among spatially-adjacent pure pixels.
The procedure of the typical and neighbor endmember selection method is as follows. Firstly, all the endmembers are extracted using the multi-index method as stated in Section 3.1. Then, the spectra of typical endmembers are calculated as the means of all the endmembers for each class respectively. For a mixed water pixel, the neighboring endmembers are searched within a 9 × 9 grid. Every possible combination of the spectra of the neighboring endmembers and those of the typical endmembers are all tried in MESMA. Finally, the best-fit model (lowest RMSE) together with the corresponding subpixel fraction are assigned to the targeted pixel.

Cleaning of the Mountain Shadow Area
To remove mountain shadow pixels, cleaning of the water maps was processed. Shaded relief maps were computed using the SRTM DEM data according to the sun elevation and the azimuth angle at acquisition time of MODIS images in the ENVI v5.3. The shaded areas were extracted by the criteria that the relief value is less than 0.1. Then a shadow map was produced by applying the threshold and was finally masked out from the water map.

Results
In this study, the accuracy assessment of the water fraction unmixing algorithm was conducted at both the local scale and regional scale. At local scale, Landsat 8 OLI data were used as accurate reference maps of the water cover to evaluate the subpixel water fraction estimates. The Landsat satellites images are the most widely used optical satellite images to extract and analyze surface water information. Two tests were conducted: one was with the MOD09A1 data; the other was with the simulated MODIS reflectance data from the OLI image. At regional scale, the water map of the Tibetan Plateau was compared with four public data sets: the TP lake data set, the GLWD, the ESA surface water data, and MOD44W.

Validation of the Algorithm with Pairs of OLI/MODIS Images
To evaluate the accuracy of the water maps, eight pairs of OLI/MODIS image subsets (see Table 1) were selected, which cover the three largest lakes in the TP, namely Qinghai Lake, Selin Co, and Nam Co, as well as five other large lakes. The acquisition dates of the OLI and MODIS data are within eight days. The OLI water maps were produced from the maximum likelihood classification with the Landsat 8 OLI data, followed by manual removal of the shaded area caused by mountains. The cloud and the associated shadow areas were detected using the OLI BQA data and excluded from the comparison. The OLI binary water maps with a resolution of 30 m were resampled to 25 m resolution and aggregated to 500 m, 1000, m and 2000 m resolution using an averaging method to generate the reference maps. Figure 3 shows examples of the false color OLI map and comparison of the water fraction map from MODIS (raster maps) and OLI data (vector maps) for Qinghai Lake, Selin Co, and Nam Co in UTM projection. The OLI water boundary vector features were generated by an automatic vectorization tool (ArcScan) in ArcGIS 10.2.
The MODIS and the reference OLI water maps were compared using: (1) mapped water area; (2) pixel-based correlation analyses (RMSE, bias, coefficient of determination (R 2 ), linear regression). Table 4 shows the comparison results between Landsat and MODIS derived subpixel water fraction maps for each OLI/MODIS image pair listed in Table 1. Table 4 reveals that the algorithm estimates the water area well: the error is within 5% for all eight OLI/MODIS pairs, with a total error of −1.4%. The underestimation is mainly attributed to: (1) the small water bodies with widths less than 500 m are not mapped in the MODIS Water Map; (2) although MOD09A1 provides the best possible 8-day composite Level 2 Gridded (L2G) observation, even with the absence of cloud and shadow, the cloud effect still exists in some pixels. By comparing with the QA data of MOD09A1, we found that in some areas (e.g., the red rectangle in Figure 3d), the surface water cannot be mapped due to the clouds in the MODIS image.
Pixel-based correlation analyses were used to analyze the difference between the reference Landsat water fractions (x) and the MODIS water fractions (y). At 500 m resolution, the overall RMSE is 7.86%. The RMSE of image pair no. 2 is the biggest (10.1%) partially because there are many small puddles which are not mapped in the MODIS image. As the resolution became coarser (from M500 to G2000), pixel-based correlation analysis shows the R 2 increased by between 0.02 and 0.05. The biases at three resolutions (500 m, 1000 m, and 2000 m) are the same due to the linear aggregation resampling of the pairs of maps. The total bias is −0.2%, which reveals a small underestimation of the water fraction in the spectral unmixing algorithm, which is mainly because the small water bodies were not considered in the algorithm.

Performance of the Water Mapping Algorithm at the Boundary of Water Bodies
High precision co-registration of different resolution satellite images is difficult to accomplish and may cause inaccuracy in the test. In order to focus on the accuracy of the algorithm at the boundary of water bodies, the algorithm was tested using the simulated MODIS image from Landsat 8 OLI image, which can provide perfect co-registration accuracy. The OLI sensor has bandwidths corresponding to the MODIS sensor, but its bandwidths are wider than that of MODIS (see Table 5). The consistency between the OLI data and the MODIS data has been demonstrated [63,64]. The OLI image over Ngoring Lake [the worldwide reference system (WRS): 143/036] obtained on 16/08/2014 (Figure 4a) was used as the reference data. By aggregation of the OLI data to 480 m resolution using the average resampling method, the MODIS reflectance (Figure 4b) was simulated. The 480 m coarse resolution data can achieve a perfect co-registration accuracy with the OLI data which the MODIS 500 m product cannot. The size of the simulated MODIS image is 100 × 100 pixels (48 km × 48 km). Since band 5 of MODIS data does not correspond to any band of OLI data, only six other bands were used in the test. By comparing the OLI water fraction reference map (Figure 4c) with the MODIS water fraction map (Figure 4d), it can be seen that the algorithm can delineate the lake shoreline well. The total area of water body estimated from Landsat is 687.54 km 2 , while the total water area from the simulated MODIS data is 685.82 km 2 . The main reason for the underestimations (0.25%) is because the water mapping algorithm ignores small water bodies, for example, the southwest bottomland of Ngoring Lake (the red rectangle in Figure 4c).
To evaluate the performance of the algorithm at marginal water pixels, the scatter plot of the fractional water body retrieved from simulated MODIS data versus the Landsat 8 OLI water maps is shown in Figure 5. The small isolated water bodies were removed from the comparison. The linear regression is good with RMSE of 14.7% and R 2 of 0.78. Since the comparison is performed only at the mixed pixels of the lake boundary, when all pixels of the scene participated in the comparison, the RMSE and R 2 become 5.81% and 0.98. Therefore, the result shows that the algorithm performs well on the boundary of the lakes and better for the whole scene.

Validation with the Public Data Sets
Seven MOD09A1 tiles (h23v05, h24v05, h24v06, h25v05, h25v06, h26v05, h26v06) on 2014/10/08 were used to generate the water fraction map of the Tibetan Plateau, which is shown in Figure 6. To assess the performance of the subpixel water mapping algorithm at regional scale, firstly the TP water fraction map in this work was compared with two public data sets: the TP lake data set and the GLWD. The three water maps were projected into Albers Equivalent Conical Projection, from which the latitudinal and longitudinal water area distributions were aggregated at steps of 1 • and compared with each other (see Figure 6). The MODIS water fraction map in this study matches very well with the TP lake data sets by Wan et al. [12] in both latitudinal and longitudinal distributions. However, a small difference was found between these mentioned two and the GLWD, because they reflect lake areas at different time periods and some small lakes are typically omitted in GLWD maps [12,15]. In terms of the longitudinal distribution between 93.5 • and 95.5 • East, the lake area in this study is greater than that found by Wan et al. [12], although it shows a similar distribution as the GLWD. This is mainly because human-exploited salt lakes (especially in Qaidam Basin) are not included in the TP lake data set. Figure 6. MODIS water fraction map for the Tibetan Plateau region from MOD09A1 (2014/10/08 to 2014/10/15), with comparisons of the latitudinal and longitudinal distributions of water areas derived from the three data sets: MODIS water fraction map produced in this study, Tibetan Plateau (TP) lake data set by Wan et al. [12], and Global Lake and Wetland Database (GLWD) by Lehner et al. [11].
To further evaluate the accuracy of the algorithm at the basin scale, the lake area of this study, the MOD44W in 2014, and ESA surface water data in Oct 2014 were aggregated by basins and compared with the TP lake data set, which was used as reference data since the manual compilation method can ensure relatively high accuracy of water boundary extraction. Table 6 summarizes the detailed results of the comparison. The human-exploited salt lakes, manmade reservoirs, and big rivers were masked out from both the MODIS water fraction map and the MOD44W like the TP lake data set. It should be noted that the satellite images used in the four water maps are different in both temporal and spatial resolution. The MODIS water fraction map was generated using the best MODIS observations during an 8-day period (October 8-15) at 500 m resolution; the MOD44W provides water maps using all MODIS 16-day composite images for the reported year at 250 m resolution; the ESA surface water map used Landsat Images to produce monthly open water occurrence (Oct 2014) at 30 m resolution; the TP data set utilized the multi-temporal (wet and dry season) images (mainly from the GF-1 wide-field-of-view (WFV) sensor at 16 m resolution) to generate the lake water boundaries. As shown in Table 6, the total lake areas for the four data sets are 45,125.04 km 2 (this study), 43,616.72 km 2 (MOD44W), 46,551.18 km 2 (ESA), and 46,595.08 km 2 (the TP lake data set), respectively. ESA surface water maps have the highest accuracy (−0.09) due to highest resolution (30 m) and advanced big data exploration and information extraction technology. At the MODIS resolution (500 m and 250 m), the difference of the total area (−1470.04 km 2 , −3.15%) between the MODIS water fraction map and the reference TP lake data set is less than that (−2978.36 km 2 , −6.39%) between the MOD44W and the TP lake data set. This is mainly because the water area in the Inner F and Inner E basin for the MOD44W water mask is greatly underestimated. There is also an underestimation of the lake area in this study, (e.g., Inner E, Brahmaputra and Yellow Basin), which is attributed to: (1) saline lakes (especially in the Inner E basin) which sometimes have a salt layer above the water surface, which greatly increases the reflectance of the water pixels and causes the pixels to be classified as 'nonwater'; (2) sometimes ice and snow cover the lakes, which cause the pixels to show reflectance of the snow; (3) although the MOD09A1 is processed using 8-day composite data, cloud effect still exists at some water pixels which causes underestimation; (4) as mentioned in Section 4.2, the small lakes with widths less than 500 m are not mapped by the algorithm. Using the MOD09A1 image from 17 May 2014, we also compared the MODIS water mapping results (water area 38,533.048 km 2 ) with the ESA surface water data (water area 39,894.66 km 2 ) in May 2014. The data are not listed in this paper. The area difference is −3.41%, which means that these two data (this study and ESA surface water data) can consistently reveal the seasonal variation of TP surface water to a certain extent. Area difference* refers to the area difference between the data and the TP lake data set by Wan et al. [12]. ESA: European Space Agency. Figure 7 illustrates the scatter plots of the basin lake area derived from this study versus the reference TP lake data set and the MOD44W versus the TP lake data set. From Figure 7, we can see the linear relationship between the MODIS water fraction map in this study and the TP lake data set is slightly stronger than the relationship between the MOD44W and the TP lake data set, which reveals a better water mapping accuracy.

The Performance of the Endmember Selection Strategy
To further illustrate the process of the endmember selection strategy, MOD09A1 h25v05 tile Day of Year (DOY) 233-240 in 2014 was used as an example. With the multi-index threshold method, 142,734 water endmembers, 64,373 vegetation endmembers, 137,635 barren endmembers, and 3167 snow endmembers were selected. The maps of the selected endmembers and the spectra of the typical endmembers are shown in Figure 8. It can be seen that the spectrum of each class endmember covers a certain range (within the error bar in Figure 8c), while the intra-class variability (within an endmember class) is less than the inter-class variability (among endmember classes). For water endmember class, the intra-class spectra variability at the green band (SD = 0.067) and red band (SD = 0.067) is larger than that at the other band (SD ≈ 0.03), which means that the spectral variations of the water body should be accounted for in the unmixing method. In the unmixing process, 27,502 (73.44%) water mixed pixels used neighboring endmembers and 9948 (26.56%) ones used typical endmembers.

Limitations of the Algorithm
The limitations of the method are worth noting, and include the following: (1) Due to the coarse resolution of the MODIS images and the mixed water pixel extraction strategy, the algorithm only considers water bodies with a width greater than 500 m. (2) Although MOD09A1 provides the best possible 8-day composite L2G observation in the absence of cloud and shadow, the cloud effect still exists at some pixels and causes undetected water pixels, which are the main reason for the underestimation of the water area. (3) Over some periods, the salt layer above the salt lake and the snow and ice cover on the lake (especially winter and spring) greatly increase the reflectivity of all segments and also lead to an underestimation of the water area; (4) The chlorophyll concentration, total suspended solids, and many other factors that affect water spectrum characteristics [65][66][67] can also affect the accuracy of water body mapping algorithms. The eutrophication and increase of sediment load can cause an increase of the near infrared reflectance and decrease of the NDWI, which will lead to the ineffective extraction of water members with the thresholding method. Although this effect may not be severe for relatively clear, high altitude lakes in the TP.

The Potential Implementation of the Proposed Surface Water Unmixing Algorithm
The potential applications of the proposed surface water unmixing algorithm could be: (1) Overcoming the coarse spatial resolution limit of the MODIS data and monitoring lakes area changes with higher temporal resolution (~8 days); (2) Dynamic changes of water surface can be used to evaluate climate change impacts on the TP and guide human activity in long-term trends [19]; (3) The generalization of this method can be spatially extended to the global scale.  Figure 8b) using neighboring endmembers (red) and typical endmembers (green); (d) the typical endmember spectra. Error bar is the standard deviation of the reflectance for each endmember class.

Conclusions
Prior works focus on mapping the TP surface water using high-resolution satellite data, which cannot monitor long-term changes with both high temporal resolution (e.g., 10 day) and high accuracy. This study presents a method that is very suitable for mountain areas to map open water in the TP. A multi-index threshold method was implemented to select endmembers. To overcome the confusion between the water and hill shade spectrum, two steps were implemented: (1) only the neighboring pixels around the pure water pixels participate in the unmixing process; (2) a shadow map is generated based on a SRTM DEM and is masked out from the water fraction map. The accuracy of the algorithm was assessed using the Landsat 8 OLI data (local scale) and four data sets (regional scale): the TP lake data set, the MOD44W, and the GLWD. At local scale, the comparison with the Landsat water maps revealed good accuracy. To overcome the co-registration error due to resolution differences between MODIS and OLI data, the performance of the algorithm at the boundary of the lakes was evaluated using a simulated MODIS image generated from the OLI image. The comparison illustrates the unmixed water fraction map from the simulated MODIS image is highly correlated with the water fraction map retrieved from the original OLI image (R 2 = 0.78), which reveals that the algorithm performs well at the boundary of lakes. At the regional scale, the water fraction map in 2014 matched well with the TP lake data set and the GLWD in both latitudinal and longitudinal distribution. The lake area estimation is more consistent with the reference TP lake data set (difference of −3.15%) than the MOD44W (difference of −6.39%).
The limitations of the method are worth noting. Due to the coarse resolution of the MODIS images and the mixed water pixel extraction strategy, the algorithm only considers water bodies with a width greater than 500 m. In addition, although MOD09A1 provides the best possible 8-day composite L2G observation in the absence of cloud and shadow, the cloud effect still affects some pixels, which might cause undetected water pixels. Therefore, future work will consider time-series information on MODIS images to mitigate the cloud effect and provide more temporally consistent water maps.
Author Contributions: C.L. and J.S. conceived the main idea. C.L., J.S., and X.L. performed the experiments. The manuscript was written by C.L. and improved by J.S., X.L., Z.S., and J.Z. All authors have read and agreed to the published version of the manuscript.