LakeTime : Automated Seasonal Scene Selection for Global Lake Mapping Using Landsat ETM + and OLI

The Landsat series of satellites provide a nearly continuous, high resolution data record of the Earth surface from the early 1970s through to the present. The public release of the entire Landsat archive, free of charge, along with modern computing capacity, has enabled Earth monitoring at the global scale with high spatial resolution. With the large data volume and seasonality varying across the globe, image selection is a particularly important challenge for regional and global multitemporal studies to remove the interference of seasonality from long term trends. This paper presents an automated method for selecting images for global scale lake mapping to minimize the influence of seasonality, while maintaining long term trends in lake surface area dynamics. Using historical meteorological data and a simple water balance model, we define the most stable period after the rainy season, when inflows equal outflows, independently for each Landsat tile and select images acquired during that ideal period for lake surface area mapping. The images selected using this method provide nearly complete global area coverage at decadal episodes for circa 2000 and circa 2014 from Landsat Enhanced Thematic Mapper Plus (ETM+) and Operational Land Imager (OLI) sensors, respectively. This method is being used in regional and global lake dynamics mapping projects, and is potentially applicable to any regional/global scale remote sensing application.


Introduction
Lakes are important features of the landscape at any scale.For example, lakes are crucial as local habitat for aquatic plants and animals, and as a storage component of the global water cycle [1].They are sensitive to global scale changes in climate and localized direct human intervention via damming, irrigation, etc. [2].Each lake is a discreet unit of land cover that could be as small as sub-hectare, but lakes are a global phenomenon existing in all kinds of biomes around the world [3,4].Their diversity and ubiquity make lakes important components of the Earth system and make them difficult to map at the global scale.
In the past, the large data volume and the prohibitive cost of remote sensing data impeded global lake shoreline and surface area mapping projects [5].Advances in technology, however, have made it possible to store and process the massive amount of data required for global, high-resolution lake shoreline mapping.In addition, the release of the full historical Landsat archive and the launch of Landsat 8 have secured the availability of high resolution data for the past, present, and future [6].These advances offer new and exciting opportunities for researchers to study lake area, and other land surface, changes over time at finer resolutions and over larger extents than previously possible.With these new opportunities also come new challenges that are specific to working at the global scale.
Several recent attempts at global scale lake mapping take different approaches with their source data.A longtime standard global lake surface area map is the Global Lakes and Wetlands Database (GLWD), which was produced by combining several regional lake maps and databases into one [7].Because of the multiple non-standard data sources, there are limits to the accuracy of this database, particularly in the Arctic [8].The more recent GLObal WAter BOdies (GLOWABO) database was based on the static GeoCover circa 2000 image dataset mapped lake surface area globally at a single time [9,10].Yamazaki et al. [11] produced a global water body map at 90 m resolution based on the Landsat Global Land Survey (GLS) database.The Landsat GLS database is available for four epochs at circa 1975, circa 1990, circa 2000, and circa 2005, which allow for multitemporal comparison of lake areas [12].Another global lake mapping effort by Pekel et al. [13] used the full USGS archive of all the available Landsat images and develop a pixel-based method that estimates the likelihood that a pixel will be water.None of these approaches explicitly address the seasonality of lake surface area.
One challenge facing global remote sensing is selecting which images to use in multitemporal lake shoreline and area mapping [14,15].Multitemporal image comparison generally requires that an image pair be either acquired from the same season or be somehow corrected for any seasonal change in the observed phenomena [16].If the image pair are not seasonally consistent, then seasonal variability is indistinguishable from longer term, interannual changes observed between the two images.When working at a local spatial scale it is relatively simple to apply a single date range across multiple years to select seasonally consistent images.However, at regional to global scales, there are many different local seasonal cycles across space [17].While a single date range would allow comparisons across time, comparisons across space would then be subject to seasonal interference.For instance, at the same time one region is in the rainy season, a neighboring region may be in the dry season.Therefore, date ranges for image selection must be identified for each location based on that region's local seasonality [18].
Previous methods for observing seasonally dependent variables with Landsat incorporated normalized difference vegetation index (NDVI) data to either select images or calibrate observations.Tatem et al., [14] used Landsat to monitor land use change across Africa.In order to calibrate the NDVI observations from Landsat, they compared concurrent NDVI observations from the NOAA-AVHRR sensor with long term seasonal averages.This was the method that was employed for generating the global orthorectified Landsat dataset [19,20].The Global Land Survey of 2005 (GLS-2005) employed the Landsat Scene Selection Interface (LASSI) tool to assign weights to candidate Landsat images according to variables such as image quality, cloud cover, and NDVI [15].Here, again, seasonal variation in NDVI was used as a proxy for the seasonal vegetation cycle.NDVI, however, may not be an appropriate proxy for seasonal variability in lake surface area, since vegetation greenness is not solely dependent on hydrology.In high latitude regions, for instance, spring green-up is largely driven by temperature change rather than water availability [21].We propose to use a simple water balance, modeled and interpolated from surface observations, to estimate lake surface area change over a season [22].From this dataset, we can identify ideal date ranges and select images from the Landsat archive.
The United States Geological Survey (USGS) Landsat archive includes images that are produced by sensors aboard seven satellites spread over the program's 40+ year history.The combination of a long historical data coverage period, relatively fine spatial resolution, and 16 day repeat time make Landsat the best available option for mapping the extent of lakes at the global scale [5,23,24].Currently, Landsat 8 carries the Operational Land Imager (OLI) sensor, which continues to acquire new images [25].
Here, we propose a method for selecting Landsat images for global multitemporal lake extent mapping based on seasonality in modeled water balance.The algorithm, called LakeTime, informs the selection of images that will be used to map lakes.The algorithm defines the optimal lake surface area mapping period as the month with the smallest net change in surface water storage after the month of maximum rainfall.This definition can be applied globally and is quantifiable using available meteorological data.This paper will consider two periods of Landsat data: the Landsat 7 Enhanced Thematic Mapper Plus (ETM+) SLC-on period from April 1999 to May 2003, and the Landsat 8 OLI period from May 2013 to November 2014.We applied a similar method to historical Landsat images from the circa 1975 Multispectral Scanner (MSS) and circa 1990 Thematic Mapper (TM) periods as well, but global coverage was not possible and they were not included in this paper.Historical Landsat coverage is inconsistent due to changing scheduling priorities, ground station networks, and technical limitation [26,27].Thanks to the application of a consistent long-term acquisition plan for Landsats 7 and 8, we were able to identify seasonally appropriate images from the ETM+ and OLI achieves, while maintaining near-global spatial coverage for circa 2000 and circa 2014 [28].This Landsat image selection method has been used in the production of decadal global lake dynamics products, and is broadly applicable to mapping any seasonally changing land cover features over time with Landsat [29].

Data
The hydrological seasonality was defined by long-term climate data.Precipitation, snow melt, potential evaporation, and temperature data were provided as monthly global gridded long term means  by the Center for Climatic Research at the University of Delaware at 0.5 × 0.5 deg resolution [22].Precipitation and temperature data were gridded from station observations using a spherical spatial-interpolation algorithm and were downloaded as the Terrestrial Air Temperature and Precipitation: Monthly and Annual Climatologies (V.3.02) [30].Potential evaporation and snow melt were modeled independently for each grid cell from the temperature and precipitation fields and came together in the Terrestrial Water Balance Data Archive: Regridded Monthly Climatologies (V.1.02) dataset [31].These global hydrology datasets offer a higher spatial resolution than most climate reanalysis products and because each value represents an interpolated point estimate, it better represents spatial variability in most locations.The drawback of using interpolated values is that the accuracy depends on the spatial coverage of the original station data, which is unevenly distributed and changes over time [32].Using long term monthly means, rather than contemporary data, limits this effect.
To select the specific images to be used for lake mapping, the determined optimal lake mapping month is then compared with the available Landsat images in the entire Landsat archive.The metadata for the full USGS Landsat archive were downloaded from the USGS Landsat Bulk Metadata Service website [33].Landsat metadata, which is updated daily to include new image acquisitions, contains many descriptive fields including acquisition quality, date acquired, cloud cover, which help the algorithm to identify the images to be selected.

Methods (Algorithm Description)
Since hydrology varies from place to place, the optimal lake surface area mapping period must be defined based on local hydrology at each location.A simple water balance was used to model seasonal changes in surface water storage at a monthly interval and at 0.5 × 0.5 deg spatial resolution.The simplified water balance equation was used as follows: where, Net is the change in surface water storage, P is precipitation, SM is snow melt, and E 0 is potential evaporation.Evaporation from a lake is more reasonably represented by potential evaporation than actual evaporation, because the availability of water is essentially infinite and not dependent on soil moisture.Snow melt is particularly important in the high latitudes and high altitudes where a significant percentage of annual precipitation falls as snow.When monthly precipitation and snow melt are equal to potential evaporation, the net change in surface water storage is minimal.At this seasonal point, the model indicates that lake surface areas are not growing or shrinking, and are most likely to be comparable from year to year.To avoid ice-cover periods, months with mean air temperatures near or below freezing must be avoided.Months with a mean monthly air temperature above freezing may still have enough below-freezing-days to freeze lakes, so a threshold of 5 • C ensures that lakes are ice free.The three criteria are prioritized as follows: The season must (a) be ice free; (b) have relatively stable surface water storage; and (c) be following the seasonal month with maximum rainfall.The LakeTime algorithm begins with identifying the month of maximum precipitation (P max Month) for each grid cell from the precipitation climatology dataset (Figure 1).If the change in water storage term (Net) is positive at P max Month, then the algorithm moves forward at monthly time steps until Net becomes negative.That negative Net value is then compared to the previous month's value to identify which has the minimum absolute value.The month with the Net value closest to zero is the zero crossing month (Net zero Month).Similarly, if Net is negative at P max Month, then the algorithm moves forward until Net becomes positive and compares with the previous month's Net value to find Net zero Month.If Net is always positive (i.e., rainforest) or always negative (i.e., desert), then Net zero Month is the month with the Net value closest to zero (lowest absolute value).Next, if the mean surface air temperature for Net zero Month is less than 5 • C, the algorithm moves back in time at monthly time steps until the temperature reaches the 5 • C threshold.The resulting month is the optimal lake area mapping month for that grid cell.The result of this classification performed at each 0.5 × 0.5 deg cell is a gridded map of optimal lake surface area mapping months (Figure 2a).so a threshold of 5 °C ensures that lakes are ice free.The three criteria are prioritized as follows: The season must (a) be ice free; (b) have relatively stable surface water storage; and (c) be following the seasonal month with maximum rainfall.
The LakeTime algorithm begins with identifying the month of maximum precipitation (PmaxMonth) for each grid cell from the precipitation climatology dataset (Figure 1).If the change in water storage term (Net) is positive at PmaxMonth, then the algorithm moves forward at monthly time steps until Net becomes negative.That negative Net value is then compared to the previous month's value to identify which has the minimum absolute value.The month with the Net value closest to zero is the zero crossing month (NetzeroMonth).Similarly, if Net is negative at PmaxMonth, then the algorithm moves forward until Net becomes positive and compares with the previous month's Net value to find NetzeroMonth.If Net is always positive (i.e., rainforest) or always negative (i.e., desert), then NetzeroMonth is the month with the Net value closest to zero (lowest absolute value).Next, if the mean surface air temperature for NetzeroMonth is less than 5 °C, the algorithm moves back in time at monthly time steps until the temperature reaches the 5 °C threshold.The resulting month is the optimal lake area mapping month for that grid cell.The result of this classification performed at each 0.5 × 0.5 deg cell is a gridded map of optimal lake surface area mapping months (Figure 2a).Landsat ETM+ and OLI images are distributed in World Reference System 2 (WRS-2) tiles that contain between 10 and 21, 0.5 × 0.5 deg grid cells each depending on latitude and orbital geometry.In order to determine the optimal lake surface area mapping month for each WRS-2 tile from the gridded map, the grid cells are overlaid by the WRS-2 tiles and are aggregated.The most abundant optimal month from the gridded map is chosen as the optimal month for each WRS-2 tile (Figure 2b).This aggregation method also serves to remove any small anomalies and spatially smooth the map of optimal months.These anomalies were generally the result of single grid cells with high sub-grid variability in the climate data often occurring on the borders of climatological zones or at the edge of mountain ranges.In the WRS-2 optimal month map (Figure 2b), Greenland and Antarctica were excluded from the analysis at this point because they were not intended to be included in the global lake map due to the consistent ice cover.Landsat ETM+ and OLI images are distributed in World Reference System 2 (WRS-2) tiles that contain between 10 and 21, 0.5 × 0.5 deg grid cells each depending on latitude and orbital geometry.In order to determine the optimal lake surface area mapping month for each WRS-2 tile from the gridded map, the grid cells are overlaid by the WRS-2 tiles and are aggregated.The most abundant optimal month from the gridded map is chosen as the optimal month for each WRS-2 tile (Figure 2b).This aggregation method also serves to remove any small anomalies and spatially smooth the map of optimal months.These anomalies were generally the result of single grid cells with high sub-grid variability in the climate data often occurring on the borders of climatological zones or at the edge of mountain ranges.In the WRS-2 optimal month map (Figure 2b), Greenland and Antarctica were excluded from the analysis at this point because they were not intended to be included in the global lake map due to the consistent ice cover.Once the optimal month is identified for each WRS-2 tile, the LakeTime algorithm searches through the USGS Landsat archive, using its metadata, for images that are acquired in the optimal month or one month before and after (to include more quality images), which also satisfy prescribed quality and cloud contamination criteria.For this study, the minimum image quality threshold was 4 (out of 9) and the maximum cloud contamination was 20%, but these thresholds can be easily changed to suit other applications.The scenes that pass the quality control and fall within the time window are then ranked according to their absolute number of days before or after the center of the optimal month.Although the temporal resolution of this study does not justify ranking scenes within the ideal month, the alternatives are clearly worse.Choosing scenes within the ideal month randomly would reduce the repeatability of the process and choosing scenes from the beginning to the end of the month would introduce a systematic bias.As such, the LakeTime algorithm generates a ranked list of the 10 best Landsat images for each WRS-2 tile.The listed images may come from any year Once the optimal month is identified for each WRS-2 tile, the LakeTime algorithm searches through the USGS Landsat archive, using its metadata, for images that are acquired in the optimal month or one month before and after (to include more quality images), which also satisfy prescribed quality and cloud contamination criteria.For this study, the minimum image quality threshold was 4 (out of 9) and the maximum cloud contamination was 20%, but these thresholds can be easily changed to suit other applications.The scenes that pass the quality control and fall within the time window are then ranked according to their absolute number of days before or after the center of the optimal month.Although the temporal resolution of this study does not justify ranking scenes within the ideal month, the alternatives are clearly worse.Choosing scenes within the ideal month randomly would reduce the repeatability of the process and choosing scenes from the beginning to the end of the month would introduce a systematic bias.As such, the LakeTime algorithm generates a ranked list of the 10 best Landsat images for each WRS-2 tile.The listed images may come from any year within the specified input year range (1999-2003, for instance).Some of the tiles have more seasonally acceptable images and some others may have only a few.The result can be exported as a list of scene IDs that are ready to order in bulk from the USGS website.Lastly, a human operator browses the candidate images to remove poor quality or overly cloudy images that are selected due to unreliable information in the metadata.

Results
The identified optimal lake mapping months match the general distribution of diverse climatological zones around the globe.The method performed well in both northern and southern hemispheres, including specific regions with unique seasonal timing (Figure 2a,b).The general pattern of optimal months across the globe follows global atmospheric circulation patterns.Most clearly visible is the effect of the movement of the Inter-Tropical Convergence Zone (ITCZ) summer rains that move across the equator, resulting in the zonal bands of optimal months that stretch across central Africa.High northern latitudes are generally restricted by the temperature threshold to a narrow summer window that moves earlier at higher latitudes.Islands and mountainous regions did not perform as well as larger land masses due to high variability in climate variables within the grid cell.
The high spatial autocorrelation of optimal months (Moran's I value of 0.675 calculated using first order polygon contiguity) indicates that the spatial aggregation of 0.5 × 0.5 deg grid cells to WRS-2 tiles effectively minimized anomalous optimal month classifications and spatially smoothed the dataset.This smoothing effect is shown in the difference between Figure 2a,b.
Ten representative locations across the globe (Figure 3) are presented here to verify the optimal month selection process in various climate conditions.The locations are: (a) central California (37.5  3a) experiences winter rains and a hot, dry summer with temperatures never dropping below 5 • C so the optimal month comes in the early spring.High latitude sites (Figure 3b,d) show the importance of the temperature criteria (green line) and the influence of snow melt (yellow line) on the seasonal water balance.There is a strong spring pulse of melt water that brings Net positive.The zero crossing in high latitudes is often after the temperature drops below the 5 • C threshold, so the optimal month is moved back until the temperature is high enough.Net is always negative for desert sites (Figure 3e,f) and thus the optimal month for both sites is simply the month closest to zero.Wet tropical sites (Figure 3h,i) experience little variation in temperature or evaporation (red line) so Net is driven by the summer rainy seasons.The India and Western Australia sites (Figure 3c,g) experience monsoon rainfall patterns with clearly defined wet and dry seasons, which make optimal month classification relatively simple.Three sites form a North-South transect across Africa (Figure 3e,f,j).The Sahara Desert (3e) and South Africa (3f) sites each have one very clearly defined, and opposite, seasonal cycle, while the equatorial Africa site (3j) has two rainy seasons with relatively dry periods between them.This is caused by the seasonal movement of the ITCZ, which crosses the equator twice each year, bringing rainfall with it.The LakeTime algorithm identifies the rainy season with the highest peak rainfall as the primary rainy season and works forward from that month.Available Landsat images were identified from the USGS Landsat data archive within 45 days of the middle of the optimal month for each tile.Thresholds for maximum cloud contamination (20%) and minimum image quality (4) were also considered.Global coverage of tiles with one or more image varied by region and by sensor, but was almost complete for both sensors.Global coverage was 86% for both ETM+ (1999ETM+ ( -2003) ) and OLI (2013-2014) (Table 1).The United States was given the highest priority in ETM+ acquisitions and the sensor acquisition setting was "always on", regardless of season or cloud cover (Arvidson et al. 2001).Thus, ETM+ coverage was best in the United States with 93% of ETM+ tiles and 91% of OLI tiles covered.ETM+ coverage was relatively poor in coastal equatorial regions, central China, and the Arabian Peninsula (Figure 4a).OLI (2013-2014) generally had fewer images per WRS-2 tile than ETM+ (1999-2003) (Figure 4b).
Available Landsat images were identified from the USGS Landsat data archive within 45 days of the middle of the optimal month for each tile.Thresholds for maximum cloud contamination (20%) and minimum image quality (4) were also considered.Global coverage of tiles with one or more image varied by region and by sensor, but was almost complete for both sensors.Global coverage was 86% for both ETM+ (1999ETM+ ( -2003) ) and OLI (2013-2014) (Table 1).The United States was given the highest priority in ETM+ acquisitions and the sensor acquisition setting was "always on", regardless of season or cloud cover (Arvidson et al. 2001).Thus, ETM+ coverage was best in the United States with 93% of ETM+ tiles and 91% of OLI tiles covered.ETM+ coverage was relatively poor in coastal equatorial regions, central China, and the Arabian Peninsula (Figure 4a).OLI (2013-2014) generally had fewer images per WRS-2 tile than ETM+ (1999-2003) (Figure 4b).Because neighboring tiles overlap each other, the area coverage was generally higher than the coverage based on WRS-2 tiles.Global area coverage was 98% for both ETM+ and OLI.This overlapping was particularly evident at high latitudes where the orbits are close enough that most of the land areas are covered by at least two tiles, helping to counteract the persistent cloudiness in the region.Coastal tiles which contain only a small fraction of land were also generally redundant and had little effect on the area coverage.These coastal tiles account for most of the difference between the % of tiles with ≥1 image and the % of area covered columns in Table 1.
Lastly, the output scene list from LakeTime for ETM+ was analyzed and downloaded by a human operator, who inspected the browse image for each of the ranked images and selected the best candidate for download and use in lake mapping.The operator downloaded images for all of the land surface areas south of 45 • N and logged the rank of the downloaded image for each tile.The mean rank of the downloaded image for tiles with coverage of one or more images was 2.46.

Discussion
Although developed independently and based on different criteria, the LakeTime optimal month results tend to fall within with the seasons prioritized by the Landsat Long-Term Acquisition Plan (LTAP).The LTAP was developed for Landsat 7/ETM+ and considers vegetation seasonality and cloud persistence when prioritizing images to be acquired each day [27,28].A new plan, developed for Landsat 8, called LTAP-8, was based on the original LTAP and maintains similar priorities [34].The acquisition plan assigns high priority to leaf-on seasons where green vegetation would be clearly visible and to the wettest seasons to offset image loss to cloud contamination [27].It is not enough, however, to rely on the LTAP alone to select images for hydrological applications, because it also acquires images from all the seasons for balance.The requirements for long term lake surface area mapping are far stricter than the LTAP, which must acquire images for all the potential users.The ultimate benefit is improved coverage, even after the strict parameters of the LakeTime algorithm are applied.
LakeTime output coverage is fundamentally limited by the images available in the USGS Landsat archive.This is particularly evident in coverage of historical Landsat images from the MSS and TM periods, which had significantly lower global coverage than contemporary sensors, and therefore were not included in this study.TM coverage in the 1990s was particularly poor in Siberia, due to the lack of a ground station within line-of-sight of the sensor.More modern sensors have improved onboard storage capacity to mitigate this problem.The priorities for Landsat image acquisitions before the implementation of the LTAP were also driven by political and bureaucratic forces, which changed over time.As the Landsat program was shifted from agency to agency, the coverage shifted and changed to suit the priorities and capabilities of the managing agency or company [27].Sharing of data from international collaborators has been inconsistent and many of these historical images remain unavailable through the USGS archive [26].Little can be done to improve the historical archive but the consistent application of LTAP, and LTAP-8 will build a future archive with maximum coverage and scientific utility.As the archive becomes larger in the future, coverage will improve and the image quality and cloud thresholds could be adjusted to be more selective.Additionally, the window of seasonally appropriate images could be shrunk from the current size of 90 days to a 60 day or 30 day window.
This hydrology based method for seasonal Landsat image selection is complimentary to other image selection methods focused on land cover or vegetation [14,15].Studies that involve remote sensing and mapping of terrestrial water in rivers, lakes, floodplains, etc. can use this or a similar method for minimizing the seasonal signal from long term change.River width, for example, can be measured reliably from Landsat imagery and be used to model discharge in remote rivers without in-situ measurements [35,36].Seasonality must be taken into account in these observations, and a method similar to LakeTime could be employed for selecting Landsat images.
The choice to use the Landsat WRS tiles as the fundamental geographical unit rather than either the 0.5 • climate data grid cells or some other predefined climate or bio region comes from the necessity of using Landsat images for the broader lake mapping study.This presupposes the treatment of Landsat scenes as single units that are either used or excluded as a whole.Other studies, like Pekel et al., have taken a pixel-based approach that would allow for a sub-image breakdown along edge features [13].This is a potentially valid alternative, although it pushes the limits of the resolution of climate reanalysis data [22].Indeed, mapping projects that take this approach do not generally treat seasonality explicitly as in this study [13,37].
We compared the ideal month map, aggregated to WRS tile, to the Köppen Geiger Climate Classification Map, and saw some agreement in rough delineation of broad climate zones [38].However, the climate map falls short of our needs in that seasonality is only treated locally and in a relative sense.For instance, zone BWh (hot desert climate) is the dominant climate zone for both the Sahara Desert and Central Australia because both of these regions have high annual mean temperatures and low annual total precipitation.However, as these regions are in the Northern and Southern Hemispheres, respectively, their hottest and driest months will be totally different.This study needs to have a result that ties the climate to a calendar month that can be used to select images.Ultimately, the best reason to use Landsat WRS tiles is that the map would eventually need to be applied at that fundamental geographical unit anyways for scene selection.
The water balance equation that is used in this study was simplified, missing runoff and ground storage terms entirely, due to the lack of global hydrology datasets.While this would be an oversimplification for a water balance model, the representation of surface water is broadly adequate for selecting Landsat images.The most important factor here is seasonality, so we focused on the water balance terms that experience large seasonal changes.Furthermore, the absolute values of the water balance terms are not as important because we are only interested in seasonality.

Conclusions
This paper addresses critical remote sensing image selection issues in global scale, high-resolution remote sensing.The proposed lake mapping image selection method is able to consistently identify an optimal lake surface area mapping season that (a) is ice free; (b) has relatively stable surface water storage; and (c) is after the rainy season.Candidate images selected from within this optimal season from the USGS Landsat archive for ETM+ and OLI sensors provide sufficient coverage for global lake surface area mapping.Contemporary lake mapping studies using ETM+ and OLI imagery will have 98% areal coverage, and OLI imagery will likely have 100% coverage soon if the LTAP-8 plan continues.Unfortunately, multitemporal studies based on historical Landsat imagery from TM and MSS sensors will have poor areal coverage due to large spatial gaps.
Newly available data from the Sentinel satellites, operated by the European Space Agency (ESA), and "flocks" of microsats, launched by the privately run Planet Labs offer huge potential for transformative research [39,40].This is one example of the kind of applications that are needed to extract valuable information from the vast quantities of remote sensing data.This work has been successfully implemented in regional and global lake mapping projects [29].The principles discussed here are broadly applicable to global and regional mapping projects beyond lakes and beyond the Landsat series.

Figure 1 .
Figure 1.This schematic illustrates the optimal month selection method.The algorithm detects the zero-crossing month on the net water balance curve (black dashed line) after the rainy season, during the ice free period.

Figure 1 .
Figure 1.This schematic illustrates the optimal month selection method.The algorithm detects the zero-crossing month on the net water balance curve (black dashed line) after the rainy season, during the ice free period.

Figure 2 .
Figure 2. Optimal lake surface area mapping month is shown (a) for each 0.5 × 0.5 degree grid cell and (b) aggregated to each WRS-2 tile.The largely ice covered regions in Greenland and Antarctica were excluded from this analysis.

Figure 2 .
Figure 2. Optimal lake surface area mapping month is shown (a) for each 0.5 × 0.5 degree grid cell and (b) aggregated to each WRS-2 tile.The largely ice covered regions in Greenland and Antarctica were excluded from this analysis.
• N, 120.5 • W) for the Mediterranean hot climate; (b) Barrow Alaska (71.5 • N, 156.5 • W) for the Arctic tundra climate; (c) central India (22 • N, 77.5 • E) for the Tropical savanna climate; (d) central Siberia (63.5 • N, 131 • E) for the Subarctic wet climate; (e) Sahara Desert (26 • N, 20.5 • E) for the Subtropical desert climate in the Northern Hemisphere; (f) South Africa (28.5 • S, 19.5 • E) for the Subtropical desert climate in the Southern Hemisphere; (g) Western Australia (21 • S, 125.5 • E) for the Subtropical desert climate in the Southern Hemisphere; (h) Brazilian Amazon (8 • S, 56 • W) for the Tropical wet climate; (i) Sumatra (3.5 • S, 104.5 • E) for the Tropical wet climate; and (j) Equatorial Africa (0 • N, 18.5 • E) for Equatorial wet climate.Each location plot shows the simplified water balance data and optimal month from a single 0.5 × 0.5 deg grid cell.The temperate California site (Figure

Table 1 .
Coverage of Seasonally Appropriate Landsat Scenes.