Spatio-Temporal Change of Lake Water Extent in Wuhan Urban Agglomeration Based on Landsat Images from 1987 to 2015

Urban lakes play an important role in urban development and environmental protection for the Wuhan urban agglomeration. Under the impacts of urbanization and climate change, understanding urban lake-water extent dynamics is significant. However, few studies on the lake-water extent changes for the Wuhan urban agglomeration exist. This research employed 1375 seasonally continuous Landsat TM/ETM+/OLI data scenes to evaluate the lake-water extent changes from 1987 to 2015. The random forest model was used to extract water bodies based on eleven feature variables, including six remote-sensing spectral bands and five spectral indices. An accuracy assessment yielded a mean classification accuracy of 93.11%, with a standard deviation of 2.26%. The calculated results revealed the following: (1) The average maximum lake-water area of the Wuhan urban agglomeration was 2262.17 km2 from 1987 to 2002, and it decreased to 2020.78 km2 from 2005 to 2015, with a loss of 241.39 km2 (10.67%). (2) The lake-water areas of loss of Wuhan, Huanggang, Xianning, and Xiaogan cities, were 114.83 km2, 44.40 km2, 45.39 km2, and 31.18 km2, respectively, with percentages of loss of 14.30%, 11.83%, 13.16%, and 23.05%, respectively. (3) The lake-water areas in the Wuhan urban agglomeration were 226.29 km2, 322.71 km2, 460.35 km2, 400.79 km2, 535.51 km2, and 635.42 km2 under water inundation frequencies of 5%–10%, 10%–20%, 20%–40%, 40%–60%, 60%–80%, and 80%–100%, respectively. The Wuhan urban agglomeration was approved as the pilot area for national comprehensive reform, for promoting resource-saving and environmentally friendly developments. This study could be used as guidance for lake protection and water resource management.


Introduction
Urban lakes provide numerous ecosystem services that are closely related to human well-being in both the present and future.These ecosystem services include agricultural production, fishery resources, flood mitigation, water storage, and entertainment locations [1,2].All of these ecosystem services are affected by fluctuations in the lake-water extent over time and space [3].Because the urban lake-water extent has been dramatically impacted due to urbanization, understanding the spatiotemporal patterns of lake-water extent dynamics is necessary for sustainable urban development.
The Wuhan urban agglomeration-consisting of Wuhan, Huangshi, Ezhou, Xiaogan, Huanggang, Xianning, Xiantao, Qianjiang, and Tianmen-is one of five national urban agglomerations and has extensive lake-water extent, due to its wet climate.For example, Wuhan is nicknamed "Water City" (Jiang Cheng).Over the past three decades, the Wuhan urban agglomeration has experienced a rapid expansion and has faced increasing challenges in protecting the lake water, both qualitatively and quantitatively, especially for Wuhan city [4][5][6][7][8][9].However, few studies have focused on the whole urban agglomeration.Moreover, these studies adopted imaging data for several specific years.Obviously, mapping lake-water extent dynamics is unreasonable without considering the ephemerality in the lake-water extent.A clear understanding of the spatio-temporal change of the lake-water extent in the Wuhan urban agglomeration is thus important.
Two major methods exist for extracting the water extent, based on remote-sensing data.One approach is the use of water indices, such as the Modified Normalized Difference Water Index (MNDWI) [18], the Normalized Difference Water Index (NDWI) [19], and the Automated Water Extraction Index (AWEI) [20].However, an ideal single threshold for water indices, in order to distinguish between water bodies and non-water bodies, is difficult to determine because the spectral signature of water varies in space and time [13].Automatic threshold methods have been adopted for determining the optimal value, like Otsu [21,22].Some comparative analyses between these water indices have been conducted [23,24], but the best water indices cannot be determined.The other method is classification based on a combination of spectral bands and other variables (e.g., water indices, shape features) [25].These classification methods include Support Vector Machines (SVM), Maximum Likelihood (ML) and Decision Tree (DT), et al. [26,27].It is important to note that a best classification method doesn't exist for all cases.For long-term seasonally continuous water mapping in a large area, DTs are more popular, primarily for their ease of application [26].However, the DT method performs less well than algorithms such as SVM and neural networks.But, the performance of Random Forest (RF), an improved implementation of DT, is superior to a traditional DT [26,27].Therefore, some studies employed the RF to classify the continuous remote images [13,25,28].For example, Tulbure et al. adopted RF to produce continuous surface water time series at a subcontinental scale in an Australian semi-arid region, the Murray-Darling Basin [13].
With an increasing amount of available remote sensing data, the mapping of seasonally continuous lake water is possible.This study employed a method that utilizes 1375 seasonally continuous Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper Plus (ETM+)/ Operational Land Imager (OLI) data scenes from 1987 to 2015, to map the long-term lake-water extent in the Wuhan urban agglomeration.This study is the first to evaluate the lake-water extent dynamics in space and time in the Wuhan urban agglomeration, based on seasonally continuous Landsat images.The specific objectives are: (1) to obtain the long-term seasonally continuous lake water extent; (2) to map the water inundation frequency to explore the spatial distribution of the lake-water extent; (3) to analyze the temporal variation of the lake-water extent; and (4) to explore the impacts of urbanization and climate change on the lake-water extent.This study could be used as guidance for lake protection and water resource management for the Wuhan urban agglomeration.

Study Area
The Wuhan urban agglomeration, located in the eastern part of the Hubei province, China, covers an area of 580,151.9 km 2 (Figure 1).It consists of Wuhan, Huangshi, Ezhou, Xiaogan, Huanggang, Xianning, Xiantao, Qianjiang, and Tianmen and was approved by the Chinese Government in 2008.As one of five national urban agglomerations, the Wuhan urban agglomeration plays a critical role in central China.On the one hand, the Wuhan urban agglomeration has witnessed rapid urban expansion.For example, the built-up area in Wuhan city increased by 10.9% from 1990-2000 and increased by 154.5% from 2000-2013 [7].On the other hand, rapid urban expansion has encroached upon lake water bodies.Lake wetlands and marsh wetlands decreased by 18.71% and 50.3%, respectively, from 1987-2005.

Study Area
The Wuhan urban agglomeration, located in the eastern part of the Hubei province, China, covers an area of 580,151.9 km 2 (Figure 1).It consists of Wuhan, Huangshi, Ezhou, Xiaogan, Huanggang, Xianning, Xiantao, Qianjiang, and Tianmen and was approved by the Chinese Government in 2008.As one of five national urban agglomerations, the Wuhan urban agglomeration plays a critical role in central China.On the one hand, the Wuhan urban agglomeration has witnessed rapid urban expansion.For example, the built-up area in Wuhan city increased by 10.9% from 1990-2000 and increased by 154.5% from 2000-2013 [7].On the other hand, rapid urban expansion has encroached upon lake water bodies.Lake wetlands and marsh wetlands decreased by 18.71% and 50.3%, respectively, from 1987-2005.
The Wuhan urban agglomeration is located within eight scenes (Path/Row: 121/39, 122/38, 122/39, 122/40, 123/38, 123/39, 123/40, and 124/39), and Table 1 lists the number of the Landsat TM, ETM+, and OLI scenes processed per path/row.A total of 1375 scenes were used, including 1019 scenes from the TM from 1987-2011, 180 scenes from the ETM+ from 1999-2003, and 176 scenes from the OLI from 2013-2015.After mid-late 2003, we did not use the scenes from the ETM+ with the Scan-Line-Corrector-Off (SLC-Off) problem, because this problem creates data gaps in the ETM+ imagery, thereby leading to missing data.Since 2013, the Landsat 8 OLI datasets were utilized because Landsat 5 was terminated in 2011.The temporal distribution of the adopted Landsat images for the study area is shown in Figure 2.Although these images are unevenly distributed in time, they can generally be used for the long-term mapping of surface water.Overall, the available time series were relatively adequate for our study.
The Wuhan urban agglomeration is located within eight scenes (Path/Row: 121/39, 122/38, 122/39, 122/40, 123/38, 123/39, 123/40, and 124/39), and Table 1 lists the number of the Landsat TM, ETM+, and OLI scenes processed per path/row.A total of 1375 scenes were used, including 1019 scenes from the TM from 1987-2011, 180 scenes from the ETM+ from 1999-2003, and 176 scenes from the OLI from 2013-2015.After mid-late 2003, we did not use the scenes from the ETM+ with the Scan-Line-Corrector-Off (SLC-Off) problem, because this problem creates data gaps in the ETM+ imagery, thereby leading to missing data.Since 2013, the Landsat 8 OLI datasets were utilized because Landsat 5 was terminated in 2011.The temporal distribution of the adopted Landsat images for the study area is shown in Figure 2.Although these images are unevenly distributed in time, they can generally be used for the long-term mapping of surface water.Overall, the available time series were relatively adequate for our study.

Landuse, Precipitation and DEM data
The landuse data acquired in 1995 and 2009 for the Wuhan urban agglomeration was collected, to analyze the impacts of urbanization on the lake area.The data were collected from the LUCC database of China.Moreover, to explore the correlation between the lake-water extent and precipitation, daily precipitation levels from eight stations around the Wuhan urban agglomeration from 1987-2015 were collected from the China Meteorological Data Sharing System.Then, the average precipitations of the eight stations were calculated.Moreover, the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTE GDEM) data with a 30-m spatial resolution were collected from the Geospatial Data Cloud (http://www.gscloud.cn),to aid in lake-water extraction.

Flowchart
The flowchart of extracting the lake-water extent in this study is displayed in Figure 3.It mainly consisted of five parts: the collection of feature variables, collection of training samples, image classification based on the random forest model, rule filter, and accuracy assessment.The details are described in the following sections.

Collection of Feature Variables
First, to diminish the negative influence from the images acquired over different locations and times, as well as the solar irradiance, etc. [13,25], we converted digital numbers (DN) to Top of Atmosphere (TOA) reflectance, according to previous research [29].Bands 1-5 and 7 from the TM/ETM+, and bands 2-7 from the OLI, were processed.Second, MNDWI [18], NDWI [19], EVI [30], NDVI [31], and NDMI [32] were calculated, based the TOA reflectance.These spectral indices are common for water extraction [11,33].Finally, a total of six bands and five spectral indexes were used in the model.

Landuse, Precipitation and DEM data
The landuse data acquired in 1995 and 2009 for the Wuhan urban agglomeration was collected, to analyze the impacts of urbanization on the lake area.The data were collected from the LUCC database of China.Moreover, to explore the correlation between the lake-water extent and precipitation, daily precipitation levels from eight stations around the Wuhan urban agglomeration from 1987-2015 were collected from the China Meteorological Data Sharing System.Then, the average precipitations of the eight stations were calculated.Moreover, the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTE GDEM) data with a 30-m spatial resolution were collected from the Geospatial Data Cloud (http://www.gscloud.cn),to aid in lake-water extraction.

Flowchart
The flowchart of extracting the lake-water extent in this study is displayed in Figure 3.It mainly consisted of five parts: the collection of feature variables, collection of training samples, image classification based on the random forest model, rule filter, and accuracy assessment.The details are described in the following sections.

Collection of Feature Variables
First, to diminish the negative influence from the images acquired over different locations and times, as well as the solar irradiance, etc. [13,25], we converted digital numbers (DN) to Top of Atmosphere (TOA) reflectance, according to previous research [29].Bands 1-5 and 7 from the TM/ETM+, and bands 2-7 from the OLI, were processed.Second, MNDWI [18], NDWI [19], EVI [30], NDVI [31], and NDMI [32] were calculated, based the TOA reflectance.These spectral indices are common for water extraction [11,33].Finally, a total of six bands and five spectral indexes were used in the model.

Collection of Training Samples
We collected training samples through the visual interpretation of the Landsat TM, ETM+, and OLI images.In addition, the five calculated spectral indices were also used to improve the confidence.The selected training samples were acquired from different regions, years, seasons, and land-cover types, to reduce possible biases [13,28].We reclassified these training samples into water pixels (1; i.e., reservoirs, rivers, lakes) and non-water pixels (0; i.e., building land, vegetation, cropland, bare land, cloud).The research used 567,320 pixels (46,041 water pixels and 521,279 non-water pixels) from the TM, 235,663 pixels (63,272 water pixels and 172,391 non-water pixels) from the ETM+, and 243,421 pixels (54,300 water pixels and 189,121 non-water pixels) from the OLI (Table 2).Large amounts of training samples can more comprehensively represent the spectral signature of both water and nonwater pixels, contributing to a higher model performance.

Image Classification Based on Random Forest Model
In this study, RF was chosen and implemented for two reasons.RF does not overfit when the number of trees increases and does not need any additional feature selection or model performance

Collection of Training Samples
We collected training samples through the visual interpretation of the Landsat TM, ETM+, and OLI images.In addition, the five calculated spectral indices were also used to improve the confidence.The selected training samples were acquired from different regions, years, seasons, and land-cover types, to reduce possible biases [13,28].We reclassified these training samples into water pixels (1; i.e., reservoirs, rivers, lakes) and non-water pixels (0; i.e., building land, vegetation, cropland, bare land, cloud).The research used 567,320 pixels (46,041 water pixels and 521,279 non-water pixels) from the TM, 235,663 pixels (63,272 water pixels and 172,391 non-water pixels) from the ETM+, and 243,421 pixels (54,300 water pixels and 189,121 non-water pixels) from the OLI (Table 2).Large amounts of training samples can more comprehensively represent the spectral signature of both water and non-water pixels, contributing to a higher model performance.

Image Classification Based on Random Forest Model
In this study, RF was chosen and implemented for two reasons.RF does not overfit when the number of trees increases and does not need any additional feature selection or model performance test method, since a random selection of features and test methods is built into it [34][35][36].Given these advantages, RF has been extensively used to characterize remotely sensed datasets across space and time [3,13,28].We also employed RF as our classifier.The built-in "Out-of-Bag" (OOB) accuracy, which is unbiased and can be used to substitute the cross-validation or independent test datasets [34], was used to assess the performance of water/non-water pixel classification.In the random forest model, the number of classification trees (T) and the number of features (m) at each node for spilling, is critical to the results [36].This study employed an RF classifier with 5,10,15,20,25,30,35,40,45, and 50 classification trees and chose the optimal item based on the OOB accuracy, using water and non-water samples as the input data.Our sensitivity tests indicate that 20 classification trees are best because, beyond that, the classification gains little performance improvement.For the number of features, previous studies suggested m to be the root of the total number of features [35,37].Hence, we set the number of features at four.RF was supported by the scikit-learn package 0.17.

Extracting Lakes Based on Shape Features and Water Inundation Frequency
Our random forest model is not able to distinguish lakes and rivers, and it classified some artificial ponds and paddy fields as water bodies.Given that these features have particular shape features, such as length and roundness (a shape measure that compares the area of the polygon to the square of the maximum diameter of the polygon), the study further adopted rule-based feature extraction for extracted lakes.First, the study composed a maximum water area by stacking all of the classification results.Then, the study built two shape feature rules (length >10,000 m and roundness <0.15) to identify rivers, as well as some artificial ponds and paddy fields.In addition, the study distinguished shadow and water bodies by setting the slope threshold (slope < 15 • ).These thresholds were selected as the most suitable ones after a trial-and-error process.Next, to diminish the effect of "salt and pepper noise" in the classification process, the study calculated the water inundation frequency (WIF, the number of times a pixel was flagged as a water body divided by the number of cloud-free observations per pixel, expressed as 0% to 100%).Those areas with a lower WIF (~5%) were regarded as noise and non-water bodies, although this process may also remove all traces of the low-frequency surface water that characterizes flooding and ephemeral water bodies.On the other hand, a spatial filtering with a 3 × 3 window was also employed, to reduce the effect of "salt and pepper" based on a "majority vote" rule.After that, remaining non-lake pixels are very little, and they are almost paddy fields, which can easily be distinguished from lakes.Therefore, we manually removed these pixels and obtained the composited maximum lake areas.The composited maximum lake extent was used as a mask to extract the lakes in each image.

Accuracy Assessment of Extracted Lakes
Considering that most lakes are located in the tile (path/row 123/39), we selected eight images for this tile and outlined the lake-water extent as being the true lake extent using visual interpretation.The eight selected images consisted of four TM images from 5 December 1995, 13 3).Moreover, these eight images have a good quality due to small influence of clouds.Then, we calculated the percentage of overlapped lake-water extent between the extracted lake-water extent and the true lake-water extent.

Determining Optimal Classification Trees
Figure 4 shows our sensitivity test results regarding the optimal classification trees, based on the OOB error.Overall, the OOB errors are very low, ranging from 0.005 to 0.04.When the number of classification trees is less than 20, the TM RF obtains a lower OOB error, followed by the OLI and ETM+ RF.Their OOB errors decrease quickly as the number of classification trees increases and reach a stable point (T = 20).From then on, increasing the number of classification trees has no significant impact on the OOB error.Therefore, the study set the optimal classification trees at 20 for the TM, ETM+, and OLI RF, when the number of features (m) was four.
Remote Sens. 2017, 9, 270 7 of 16 impact on the OOB error.Therefore, the study set the optimal classification trees at 20 for the TM, ETM+, and OLI RF, when the number of features (m) was four.

Accuracy of the Extracted Lakes
The accuracies of the eight images were 93.17%, 95.14%, 96.04%, 91.71%, 95.53%, 89.68%, 90.17%, and 93.41% (Table 3).The mean value and standard deviation of accuracy were 93.11% and 2.26%, respectively.The results showed that lake-water extraction based on the random forest model is feasible.

Accuracy of the Extracted Lakes
The accuracies of the eight images were 93.17%, 95.14%, 96.04%, 91.71%, 95.53%, 89.68%, 90.17%, and 93.41% (Table 3).The mean value and standard deviation of accuracy were 93.11% and 2.26%, respectively.The results showed that lake-water extraction based on the random forest model is feasible.

Lake Water Extent Changes from 1987~2015 in Different Cities
We further explored the annual maximum lake-water extent across the entire study area from 1987-2015 (Figure 6).Considering that the earlier images are much smaller in number and cannot cover the entire study area, we composited two multi-year maximum lake areas from 1987-1990 and from 1991-1994.The lake in the study area experienced an apparent decrease from 1987-2015.From 1987-2002, the lake area shows a fluctuating tendency.However, it experienced a continuous decrease from 2003-2007, followed by a fluctuating tendency from 2008-2015, but the lake area did not return to the previous condition.The calculated results revealed that the average maximum lake area was 2276.03 km 2 from 1995-1999, and it decreased to 1966.94 km 2 from 2011-2015, with a total loss of 309.09 km 2 and a percentage loss of 13.58%.

Lake Water Extent Changes from 1987~2015 in Different Cities
We further explored the annual maximum lake-water extent across the entire study area from 1987-2015 (Figure 6).Considering that the earlier images are much smaller in number and cannot cover the entire study area, we composited two multi-year maximum lake areas from 1987-1990 and from 1991-1994.The lake in the study area experienced an apparent decrease from 1987-2015.From 1987-2002, the lake area shows a fluctuating tendency.However, it experienced a continuous decrease from 2003-2007, followed by a fluctuating tendency from 2008-2015, but the lake area did not return to the previous condition.The calculated results revealed that the average maximum lake area was 2276.03 km 2 from 1995-1999, and it decreased to 1966.94 km 2 from 2011-2015, with a total loss of 309.09 km 2 and a percentage loss of 13.58%.Among the eight cities, Wuhan experienced the most serious lake-water area loss from 1987-2015.From 1995-1999, the average lake-water area was 803.01 km 2 , and the area decreased to 688.18 km 2 from 2011-2015, with a total loss of 114.83 km 2 or 14.30%.Huanggang, Xianning, and Xiaogan also experienced relatively serious losses, totaling an area of loss of 44.40 km 2 , 45.39 km 2 , and 31.18km 2 , and a percentage loss of 11.83%, 13.16%, and 23.05%, respectively.
Among the eight cities, Wuhan experienced the most serious lake-water area loss from 1987-2015.From 1995-1999, the average lake-water area was 803.01 km 2 , and the area decreased to 688.18 km 2 from 2011-2015, with a total loss of 114.83 km 2 or 14.30%.Huanggang, Xianning, and Xiaogan also experienced relatively serious losses, totaling an area of loss of 44.40 km 2 , 45.39 km 2 , and 31.18km 2 , and a percentage loss of 11.83%, 13.16%, and 23.05%, respectively..

Decreasing Lake Water Extent Due to Urbanization in Urban Areas
From 1995 to 2009, about 193.85 km 2 of lake bodies have been changed into agricultural land (paddy fields), and about 16.60 km 2 of lake bodies have been changed into building land in the Wuhan urban agglomeration, which indicates that a rapid urbanization over the past several years is an important reason for the loss of lake areas, especially for Wuhan [7,8].Therefore, we selected five typical urban lakes (Tangxunhu, Donghu, Yanxihu, Nanhu, and Shahu) in Wuhan to explore the phenomena (Figures 7 and 8, Table 4).These five urban lakes are majorly encroached by building land.Note: Loss rate = (Areai−Areaj)/(j−i).Areai represents lake area in the previous year i, Areaj represents lake area in the later year j.Therefore, the positive value indicates the decrease while the negative value indicates the increase.

Decreasing Lake Water Extent Due to Urbanization in Urban Areas
From 1995 to 2009, about 193.85 km 2 of lake bodies have been changed into agricultural land (paddy fields), and about 16.60 km 2 of lake bodies have been changed into building land in the Wuhan urban agglomeration, which indicates that a rapid urbanization over the past several years is an important reason for the loss of lake areas, especially for Wuhan [7,8].Therefore, we selected five typical urban lakes (Tangxunhu, Donghu, Yanxihu, Nanhu, and Shahu) in Wuhan to explore the phenomena (Figures 7 and 8, Table 4).These five urban lakes are majorly encroached by building land.Note: Loss rate = (Area i −Area j )/(j−i).Area i represents lake area in the previous year i, Area j represents lake area in the later year j.Therefore, the positive value indicates the decrease while the negative value indicates the increase.The study observed that Shahu and Nanhu have been declining from 1987-2015.The total areas of loss of Shahu and Nanhu were 2.23 km 2 and 3.35 km 2 , respectively, from 1995-2015, representing 47.35% and 31.31% of the lake-water areas in 1995.Yanxihu showed a relatively stable tendency before 2003 and then experienced a rapid declining tendency from 2013.The total area of loss of Yanxihu was 0.58 km 2 from 1995-2015, accounting for 4.61% of the lake-water area in 1995.Donghu and Tangxunhu displayed a fluctuating tendency before 2003.From 2003-2010, the two cities experienced a slight declining tendency, followed by a rapid declining tendency after 2010.The total loss areas in Donghu and Tangxunhu were 3.10 km 2 and 2.54 km 2 , respectively, from 1995-2015, accounting for 9.21% and 5.71% of the lake-water areas in 1995.

Lake Water Extent Changes Due to Climate Changes in Non-Urban Areas
Precipitation variation plays a critical role in lake areas.This study chose Liangzihu and Futouhu as non-urban lakes where urbanization has had little impact, to calculate the Pearson's correlation coefficient between lake areas and the annual maximum precipitation from 1987-2015.In addition, the correlation of the entire study area was also calculated.
These two non-urban lakes showed a significant correlation (p-values < 0.05).For Futouhu and Liangzihu, their r-values in area were 0.6771 and 0.6558, respectively.Overall, the lake-water areas of the entire study area showed a significant correlation with precipitation.It is evident from Figure 9 that much rainfall in the study area in 1996, 1998, 2002, and 2010 resulted in floods, which increased the lake-water area (Figure 10), while little rainfall in the study area in 1997, 2001, 2006, and 2011 decreased the lake-water area (Figure 11), because of droughts.From 1999-2001 and from 2002-2006, rainfall in the study area continuously decreased, resulting in a continuous decrease in the lake-water area, whereas, from 2007-2010, rainfall in the study area increased, resulting in an increase in the lake-water area.Moreover, for Futouhu, we found that relatively regular enclosure nets are distributed in Lake Futohu from Figure 5, indicating the influence of human activities.Therefore, the stable and increased extent of Lake Futouhu from 2004 to 2006, while the precipitation was decreasing, is possible.
The study observed that Shahu and Nanhu have been declining from 1987-2015.The total areas of loss of Shahu and Nanhu were 2.23 km 2 and 3.35 km 2 , respectively, from 1995-2015, representing 47.35% and 31.31% of the lake-water areas in 1995.Yanxihu showed a relatively stable tendency before 2003 and then experienced a rapid declining tendency from 2013.The total area of loss of Yanxihu was 0.58 km 2 from 1995-2015, accounting for 4.61% of the lake-water area in 1995.Donghu and Tangxunhu displayed a fluctuating tendency before 2003.From 2003-2010, the two cities experienced a slight declining tendency, followed by a rapid declining tendency after 2010.The total loss areas in Donghu and Tangxunhu were 3.10 km 2 and 2.54 km 2 , respectively, from 1995-2015, accounting for 9.21% and 5.71% of the lake-water areas in 1995.

Lake Water Extent Changes Due to Climate Changes in Non-Urban Areas
Precipitation variation plays a critical role in lake areas.This study chose Liangzihu and Futouhu as non-urban lakes where urbanization has had little impact, to calculate the Pearson's correlation coefficient between lake areas and the annual maximum precipitation from 1987-2015.In addition, the correlation of the entire study area was also calculated.
These two non-urban lakes showed a significant correlation (p-values < 0.05).For Futouhu and Liangzihu, their r-values in area were 0.6771 and 0.6558, respectively.Overall, the lake-water areas of the entire study area showed a significant correlation with precipitation.It is evident from Figure 9 that much rainfall in the study area in 1996, 1998, 2002, and 2010 resulted in floods, which increased the lake-water area (Figure 10), while little rainfall in the study area in 1997, 2001, 2006, and 2011 decreased the lake-water area (Figure 11), because of droughts.From 1999-2001 and from 2002-2006, rainfall in the study area continuously decreased, resulting in a continuous decrease in the lake-water area, whereas, from 2007-2010, rainfall in the study area increased, resulting in an increase in the lake-water area.Moreover, for Futouhu, we found that relatively regular enclosure nets are distributed in Lake Futohu from Figure 5, indicating the influence of human activities.Therefore, the stable and increased extent of Lake Futouhu from 2004 to 2006, while the precipitation was decreasing, is possible.

Long-Term Lake Water Extent Dynamics from Landsat Time Series Data
This study applied all available Landsat time series data from 1987-2015 to generate a comprehensive historical record of lake-water extent dynamics (1987-2015) at a high spatial resolution (30 m), for the entire Wuhan urban agglomeration.The study justified using the random forest model based on the multi-feature variables to extract water bodies from massive remote-sensing data.Tulbure [13] also applied the same method to produce continuous surface water-time series at the subcontinental scale in an Australian semi-arid region, the Murray-Darling

Long-Term Lake Water Extent Dynamics from Landsat Time Series Data
This study applied all available Landsat time series data from 1987-2015 to generate a comprehensive historical record of lake-water extent dynamics (1987-2015) at a high spatial resolution (30 m), for the entire Wuhan urban agglomeration.The study justified using the random forest model based on the multi-feature variables to extract water bodies from massive remote-sensing data.Tulbure [13] also applied the same method to produce continuous surface water-time series at the subcontinental scale in an Australian semi-arid region, the Murray-Darling

Long-Term Lake Water Extent Dynamics from Landsat Time Series Data
This study applied all available Landsat time series data from 1987-2015 to generate a comprehensive historical record of lake-water extent dynamics (1987-2015) at a high spatial resolution (30 m), for the entire Wuhan urban agglomeration.The study justified using the random forest model based on the multi-feature variables to extract water bodies from massive remote-sensing data.Tulbure [13] also applied the same method to produce continuous surface water-time series at the subcontinental scale in an Australian semi-arid region, the Murray-Darling Basin, with an overall classification accuracy of 99.9%.Given that the water-body types and land-use types in the Wuhan urban agglomeration are more complicated than they are in the Murray-Darling Basin, the mean value (93.11%) and standard deviation (2.26%) for this study are reasonable.In addition, this study did not select the same feature variables as those in Tulbure's study.For example, in Tulbure's study, AWEI [20] was used for the TM and ETM+ images.However, in this study, the use of AWEI is not feasible for the OLI images because the SWIR band is different from the TM band.

Decreasing Lake Water Extent Dynamics Under Urbanization Process and Climate Changes
From 1987-2015, the lake-water extent in the Wuhan urban agglomeration suffered from a substantial decrease, for two reasons.One import factor is urbanization.Due to rapid urban development in the Wuhan urban agglomeration over the past several years, the lake was converted into built-up areas, causing a substantial decrease in the lake-water extent (Figure 8).Chen et al. [7] revealed that the population is the main driving factor causing the morphological change of an urban lake in Donghu, Wuhan.Increasing numbers of people flowed into the urban area.Thus, increasing numbers of urban lands were needed, causing land-use changes.However, the transformation did not occur simultaneously.For example, Shahu is close to the city center, so it suffered a decrease in the lake-water extent at the earliest point in time.As the urban area expanded outward, lakes such as Tangxunhu, that were once far from the city's center, began to experience damage due to the conversion into built-up areas.In 2008, the Wuhan urban agglomeration was approved as the national comprehensive reform pilot area for promoting resource-saving and environmentally friendly developments.Increasing attention has been placed on environmental and ecological issues, and lake protection has been contained in development planning.Therefore, the lake's loss rate from 2008-2015 was lower than that from 1987-2015, for Tangxunhu, Yanxihu, Nanhu, and Shahu (Table 4).Another more important factor is climate change, especially the variation in precipitation.The total lake-water extent is strongly correlated with annual precipitation (Figure 9).For example, in 2011, the lake area of the Wuhan urban agglomeration reached the minimum value.At the same time, annual precipitation also reached the relative minimum value (Figure 9).In addition, this study found that annual precipitation showed a decreasing trend, which may explain the decrease in lake area from the perspective of climate change (Figure 9).

Uncertainties and Prospects
Some uncertainties exist in this study's mapping of long-term lake-water extent.One uncertainty is due to cloud cover.Cloud cover not only reduces the number of clear observations, but also causes cloud shadow that is difficult to distinguish from water bodies.The assumption that clouds cannot persist at the same place [38] could diminish the adverse effect by removing low-frequency water bodies.Moreover, atmospheric effects play a very critical role in datasets across time and space.However, due to the high cost of atmospheric correction in time and process, the study didn't pursue atmospheric correction.However, TOA reflectance was used as the feature variable to correct for the difference of the sun-earth distance, exoatmospheric solar irradiance, and solar geometry between images [25].Future study will focus on the impacts of urbanization on the lake landscape and the resulting adverse influences, based on long-term lake mapping.

Conclusions
As sustainable development is vigorously advocated by the government, lakes will play an increasing role in providing urban ecosystem services for the Wuhan urban agglomeration.Therefore, a clear understanding of lake-water extent dynamics is necessary.This research revealed the lake-water extent dynamics in the Wuhan urban agglomeration, based on seasonally continuous Landsat TM/ETM+/OLI data from 1987-2015 and the random forest model.Overall, this study provides the following conclusions:

Figure 1 .
Figure 1.The location of the study area.

Figure 1 .
Figure 1.The location of the study area.

Figure 2 .
Figure 2. Temporal distribution (day of year) of Landsat images used in this study.

Figure 2 .
Figure 2. Temporal distribution (day of year) of Landsat images used in this study.

Figure 3 .
Figure 3. Flowchart of the study.

Figure 3 .
Figure 3. Flowchart of the study.
February 2004, 22 July 2004, and 9 September 2009; two ETM images from 11 January 2000 and 22 July 2000; and two OLI images from 17 September 2013 and 26 November 2015, accounting for four or five year intervals, season changes, and sensor influences (Table

Figure 4 .
Figure 4. Impacts of different numbers of trees in the Random Forest based on OOB error assessment.

Figure 4 .
Figure 4. Impacts of different numbers of trees in the Random Forest based on OOB error assessment.

Figure 8 .
Figure 8. Lake water area changes from 1987~2015.First column represents the water inundation frequency (WIF) from 1987 to 2015 for these five lakes.The other four columns represent the composited maximum lake-water extents of these five lakes at different years: 1995, 2000, 2008, and 2015.

Figure 8 .
Figure 8. Lake water area changes from 1987~2015.First column represents the water inundation frequency (WIF) from 1987 to 2015 for these five lakes.The other four columns represent the composited maximum lake-water extents of these five lakes at different years: 1995, 2000, 2008, and 2015.

Figure 8 .
Figure 8. Lake water area changes from 1987~2015.First column represents the water inundation frequency (WIF) from 1987 to 2015 for these five lakes.The other four columns represent the composited maximum lake-water extents of these five lakes at different years: 1995, 2000, 2008, and 2015.

Figure 9 .
Figure 9. Correlation between annual maximum lake areas and annual precipitation from 1987 to 2015.Figure 9. Correlation between annual maximum lake areas and annual precipitation from 1987 to 2015.

Figure 9 .
Figure 9. Correlation between annual maximum lake areas and annual precipitation from 1987 to 2015.Figure 9. Correlation between annual maximum lake areas and annual precipitation from 1987 to 2015.

Figure 10 .
Figure 10.Increased water area due to floods in four periods: 1996, 1998, 2002, and 2010.The lake water area refers to the composited maximum area in the same year.

Figure 11 .
Figure 11.Decreased water area due to droughts in four periods: 1997, 2001, 2006, and 2011.The lake water areas refer to the composited maximum area in the same year.

Figure 10 .
Figure 10.Increased water area due to floods in four periods: 1996, 1998, 2002, and 2010.The lake water area refers to the composited maximum area in the same year.

Figure 10 .
Figure 10.Increased water area due to floods in four periods: 1996, 1998, 2002, and 2010.The lake water area refers to the composited maximum area in the same year.

Figure 11 .
Figure 11.Decreased water area due to droughts in four periods: 1997, 2001, 2006, and 2011.The lake water areas refer to the composited maximum area in the same year.

Figure 11 .
Figure 11.Decreased water area due to droughts in four periods: 1997, 2001, 2006, and 2011.The lake water areas refer to the composited maximum area in the same year.

Table 1 .
Number of Landsat TM, ETM, and OLI scenes processed per path/row.

Table 1 .
Number of Landsat TM, ETM, and OLI scenes processed per path/row.

Table 2 .
Collected sample size of different sensor types.

Table 2 .
Collected sample size of different sensor types.

Table 3 .
Accuracy of extracted lakes based on the random forest model.

Table 3 .
Accuracy of extracted lakes based on the random forest model.

Table 4 .
Lake-water extent changes of typical urban lakes.

Table 4 .
Lake-water extent changes of typical urban lakes.