Using Radiometric and Categorical Change to Create High-Accuracy Maps of Historical Land Cover Change in Watersheds of the Great Lakes Basin

Abstract

Great Lakes Basin landscapes are undergoing rapid land cover and land use (LCLU) change. The goal for this study was to identify changes in land cover occurring in the Great Lakes Basin over three time periods to provide insights into historical land cover changes occurring on a bi-national watershed scale. To quantify potential impacts of anthropogenic changes on important yet vulnerable Great Lakes Wetland ecosystems, the historical changes in land cover over time are assessed via remote sensing. The goal is to better understand legacy effects on current conditions, including wetland gain and loss and the impacts of upland ecosystems on wetland health and water quality. Three key time periods with respect to Great Lakes water level changes and coastal wetland plant invasions were mapped using Landsat-derived land cover maps: 1985, 1995, and 2010. To address change between the three time periods of interest, we incorporate both radiometric and categorical change analysis and open-source tools available for assessing time series data including LandTrendr and TimeSync. Results include maps of annual land cover transition from 1985 to 1995 and 1995 to 2010 basin-wide and by ecoregion and an assessment of the magnitude and direction of change by land cover type. Basin-wide validated change results show approximately 776,854 ha of land changed from c.1980–1995 and approximately 998,400 ha of land changed from c.1995–2010. Both time periods displayed large net decreases in both deciduous forest and agricultural land and net increases in suburban cover. Change by ecoregion is reviewed in this study with many of the change types in central plains showing change in and out of agriculture and suburban land covers, the mixed wood plain ecoregion consisted of a mixture of agricultural, suburban, and forestry changes, and all top five change types in the mixed wood shield consisted of various stages of the forestry cycle for both time periods. In comparison with previous LCLU change studies, overall change products showed similar trends. The discussion reviews why, while most changes had accuracies better than 84%, accuracies found for change from urban to other classes and from other classes to agriculture were lower due to unique aspects of change in these classes which are not relevant for most change analyses applications. The study found a consistent loss in the deciduous forest area for much of the time studied, which is shown to influence the aquatic nitrogen implicated in the expansion of the invasive plant Phragmites australis in the Great Lakes Basin. This underscores the importance of LCLU maps, which allow for the quantification of historical land change in the watersheds of the Great Lakes where invasive species are expanding.

1. Introduction

Vast areas of land in the Great Lakes Basin (GLB) have been undergoing rapid land cover and land use (LCLU) change over the past several decades. LCLU in the watersheds that drain to the coasts have a huge impact on the inflow of nutrients and the biogeochemistry of the coastal zone and coastal wetlands. Despite their importance, more than two-thirds of the wetlands in the GLB have been drained for agriculture and other development [1] since the turn of the last century, and the remaining wetlands are subject to several threats including land conversion, climate change, and invasion by non-native plant species. Land transformations have altered groundwater levels and stream flows, increased nutrient loadings to water bodies, and contributed to anthropogenic contamination that impacts ecosystems at a variety of scales [2]. The Great Lakes are highly susceptible to eutrophication by both nitrogen (N) and phosphorus (P). While coastal wetlands provide a key ecosystem service of retaining inflowing nutrients [3,4], increases in nutrient loading are well-documented to lead to increases in wetland-invasive plants such as phragmites (Phragmites australis) and cattail (Typha angustifolia and T. x glauca) [5,6,7,8,9]. Influxes of coastal invasion are also tied to changing Great Lakes water levels, which have a historical periodicity of about 10 years, but climate change is increasing the variability over shorter time scales. For example, a record low observed in 2013 rebounded to historically high-water levels nearly 18 months later. Water level fluctuations influence coastal wetland extent and function, with implications to the source/sink nature of these invaded ecosystems. Historically, high-water periods (e.g., 1986, 1997) have been followed by waves of invasion into barren, exposed sediments after water levels declined [10,11,12,13]. These invasive species form monocultures in vast areas, reducing biodiversity, changing wildlife and fisheries habitat, and altering ecosystem functions such as nutrient retention and carbon sequestration [8].

In real time, we are seeing changes due to both climate and land use occurring within the GLB, altering ecosystems and affecting the habitat, water quantity and quality—both surface and groundwater [14]. Recent examples include the following: the rapid fluctuations of Great Lake water levels and erosion events [15]; large algae blooms in Lake Erie [16] and other freshwater lakes [17]; contaminated drinking water in Flint and Benton Harbor, which disproportionately affect minorities [18]; oil spills [19]; and urban and agricultural runoff pollution [20]. Timely data and fine resolution imagery are often needed to address the events and understand ecosystem change.

Invasive species represent one of the most significant and immediate threats to the ecosystems and economies of the GLB. Researchers have identified the Great Lakes as one of the most invaded freshwater systems across the world [21,22]. To date, 188 aquatic non-native species have been established in the Great Lakes, of which 34% are considered invasive [23]. Invasive species can alter ecosystem functions including nutrient (nitrogen [N] and phosphorus [P]) retention and carbon sequestration [8,24]. Invasive plants such as Phragmites often form monocultures that typically have greater biomass and productivity than the native vegetation they replace [25], thus altering nutrient and carbon dynamics. Aquatic and terrestrial invasive species in the GLB alter food webs, degrade the habitat, and outcompete and/or cause direct mortality or injury to native species. These impacts have threatened the native biodiversity and abundance, changed the community structure and dynamics of ecosystems in the GLB, and reduced the health and resilience of many of these ecosystems. Because coastal wetlands act as a buffer to the Great Lakes, their invasion can directly influence N, P, and C loading into the Great Lakes. Increases in nutrient loading can lead to increases in wetland-invasive plants such as Phragmites australis and Typha spp. (cattail) [5,6,8,9], as well as harmful algal blooms (HABs) such as Cladophora [26].

Here, we assess the historical changes in LCLU through time via remote sensing towards a goal of better understanding the legacy effects of upland conversions on current coastal conditions in the Great Lakes watersheds of Lakes Huron, Erie and Michigan, the lakes that are most heavily infested with invasive Phragmites. The objective of the study is to quantify basin-wide land cover changes to link these changes to the observed expansion of invasive wetland plants and better understand wetland ecosystem dynamics through process-based modeling efforts [5]. This requires sufficient temporal resolution of change as well as a basin-wide approach. Previous assessments of LCLU change in the Great Lakes have been addressed over portions of the GLB [27,28] or U.S. only [29,30] for various ecological, economic, and social studies. Until recently, changes within the full GLB had never been evaluated [31]. This is partially due to the time and effort required to create land change products over such a vast area from high to moderate resolution data. Cloud computing, such as the Google Earth Engine (GEE), provides the platform and datasets needed to assess inter-annual changes that we and others (e.g., [31]) have recently leveraged for change analysis.

To assess legacy changes in land use and land cover on coastal wetland invasions, Great Lakes water level changes and coastal wetland plant invasions were mapped for three key time periods: 1985, 1995, and 2010. The c.1985 land cover map pre-dates the 1987–1990 declining Great Lake water levels [32] and influx of wetland plant invasion (1987 onward), while c.1995 captures conditions after the 1987–1990 water level decline, during the resurgence in water levels (peak 1997). The c.2010 baseline data capture conditions near the end of a 14-year low-water-level period (1993–2012) when invasive plants rapidly expanded to their current range, and prior to large, localized efforts to control the invasive Phragmites australis in the more southern areas of the Basin. Our goal for this study was to identify the land cover change occurring in the GLB between our three time periods of interest (1985–1995 and 1995–2010) to provide the first insights into historical changes in land cover (urban, forestry and wetlands) occurring on a bi-national watershed scale.

2. Materials and Methods

2.1. Study Region and Approach

The area of interest for this study (Figure 1) consists of all the areas inside the three lake basins most heavily affected by invasive wetland plants (Phragmites and Typha spp.): Huron, Erie, and Michigan. This region consists of 515,118 km² of seven states (Illinois, Indiana, Michigan, New York, Ohio, Pennsylvania, Wisconsin), and the province Ontario, Canada.

To create seamless change maps of urban, forest, and wetlands for the time periods of interest for the United States and Canada with high accuracy, we chose to use open-source tools, LandTrendr [33] and TimeSync [34], combined with the power of cloud computing (GEE), incorporating both radiometric and categorical change analyses to reduce the false change that often occurs when categorical change is used alone [35]. In addition, the new tools allow for the assignment of years of change to each change polygon, which is an important input to the hydrological modeling.

2.2. Categorical Map Data Sources

For the categorical change, we used existing maps. These maps were designed to distinguish important wetland types, including the delineation of the invasive species Phragmites and Typha spp., while also mapping upland land use and land cover. They were produced using a consistent method across the bi-national GLB with accuracies verified through the reservation of field training data via stratified random sampling.

Land cover types were derived from Anderson level II [36] and augmented with wetland classes common in the Great Lakes region. The minimum mapping unit for the data products was 0.2 hectares due to the spatial resolution of the satellite remote sensing data used to derive the maps [37,38].

For the most recent time period of interest, 2010, we used a baseline map shown in Figure 1. This product includes both the United States and Canada and was previously created by our team from multi-date L-band Synthetic Aperture Radar (SAR) and optical-IR reflectance data [37,38]. This map is available upon request from the authors or from a request made through the on-line data viewer (https://geodjango.mtri.org/coastal-wetlands/, accessed on 14 June 2024). Overall accuracy of over 93% was achieved, with high accuracies for upland types and some wetland types showing lower accuracies (

Table 1. Producers and users’ accuracy value for each class in the c.2010 map. Colors of the classes correspond to the classes mapped in Figure 1.

Classified	Producers Accuracy (%)	Users Accuracy (%)
Urban	77.5%	74.2%
Grass	79.0%	50.9%
Agriculture	88.6%	97.6%
Orchard	84.8%	82.7%
Forest	92.0%	88.4%
Pine Plantation	80.9%	89.7%
Shrub	84.2%	80.4%
Barren Light	82.3%	88.1%
Water	99.3%	100.0%
Aquatic Bed	91.3%	61.3%
Marsh	73.1%	58.7%
Typha	87.1%	81.1%
Phragmites	83.7%	69.3%
Open Peatland	65.0%	73.7%
Shrub Peatland	56.5%	65.6%
Treed Peatland	69.2%	65.3%
Wetland Shrub	72.9%	72.5%
Forested Wetland	81.0%	78.5%
Overall Accuracy	93.3%

).

Most past efforts to map LCLU pre-2010 were completed for the United States or Canada separately using unique data sources and methods and did not cross the border. To obtain bi-national categorical land cover for 1985 and 1995, we compiled categorical data from multiple sources derived from Landsat optical-IR imagery (Landsat). For the United States side of the three Lake basins, we used NOAA’s Coastal Change Analysis Program (C-CAP) products [32]. For Canada, Agriculture and Agri-Food Canada (AAFC) [39] land cover products from 1990 and 2000 (Government of Canada 2018) were available, but we needed to match the 1985 and 1995 time period maps of C-CAP. To do this, we utilized a modified technique developed by Kennedy et al. [33] to create interpolated Canadian land cover products for c.1985 and c.1995 (Figure 2). This methodology allowed us to adjust the AAFC maps to be more consistent with the other land cover products as well. Temporally smoothed Landsat composites produced using the LandTrendr tools for c.1985 and c.1995 were classified using visually interpolated training data and a random forest classifier. These temporally smoothed composites emphasize the regions where changes in land cover occurred between the two time periods. The AAFC maps were used to assist the image interpreter in creating the training datasets. While the products were being developed, our team worked directly with Canadian collaborators to ensure that we were interpreting the classes correctly from the source maps and for an initial verification of the new interpolated maps.

All three land cover products were reclassified into common land cover classes to facilitate comparison (

Table 2. Crosswalk of the categorical map classes to desired classes for the hydrological modeling analysis. Colors of the classes correspond to the classes mapped in Figure 1.

CCAP	MTRI	Canadian	Reclassified	Class Number
Developed high	Urban	Settlement	Urban	1
Developed medium	Road
Developed low	Urban grass	Roads	Suburban	2
Developed open space	Suburban
Bare land	Bare land	Other land	Bare land	3
Cultivated crops	Agriculture	Cropland	Agriculture	4
Pasture/hay	Orchard
Grassland/herbaceous	Fallow	Grassland managed	Grasslands	5
Deciduous	Forest	Forest	Deciduous	6
Mixed Forest		Trees
Evergreen	Pine Plantation		Evergreen	7
Shrubs	Shrubs	Grassland unmanaged	Shrubs	8
Palustrine forest	Forested wetland	Forested wetland	Woody wetland	9
Estuarine forest		Treed wetland
Palustrine shrub	Wetland Shrub	Wetland shrub
Estuarine shrub
Palustrine wetland	Emergent wetland	Wetland	Wetland	10
Estuarine wetland	Schoenoplectus	Wetland herb
	Typha
	Phragmites
	Peatlands
Palustrine aquatic bed	Aquatic bed		Aquatic bed	11
Unconsolidated shore
Water	Water	Water	Water	12

). Classification schemes from the different land cover datasets were reclassified into a generalized twelve-class scheme for change analysis; classes were chosen to preserve as much detail as possible while matching the standardized land cover classification system developed by Anderson et al. [36]. Classes were also chosen to best represent hydrological processes occurring on the landscape so they could be used to assist in parameterizing hydrologic models. Reclassified maps were mosaiced together to produce basin-wide land cover products for 1985, 1995, and 2010.

2.3. Change Detection

The more common approach to land cover change assessment uses the difference between two land cover classifications and report on categorical changes; however, the differences between two categorical maps caused by map production methods and varying class definitions can lead to the detection of false change. Integrating radiometric magnitude change information, derived from a change vector analysis (CVA) [35], into traditional categorical change methods can improve accuracy by eliminating some of the false change present [35,40]. Previous studies have integrated radiometric change data into hybrid categorical/radiometric change methodologies through either a uniform threshold approach [41], automatic threshold detection [42], or user-defined change class-specific thresholds [43]. All approaches have their unique strengths and weaknesses that best suit them for specific problem sets. The use of uniform and automatic thresholds is efficient but best-suited for the detection of a single type of land cover change as they are not flexible in capturing a range of spectral responses. Class-specific change thresholds, where the user defines the radiometric threshold expected for each possible land cover change, are time- and labor-intensive but can capture a range of spectral responses to account for all possible changes in land cover observed in an analysis. We chose this latter approach to integrate radiometric change information into our hybrid categorical/radiometric change methodology because our goal was to quantify all land cover change occurring within the GLB, and a single radiometric threshold would not have been sufficient in capturing the regional urban, forestry, and wetland changes.

Change detection methodologies, categorical, radiometric, and hybrid categorical/radiometric, can all be applied on a pixel-by-pixel basis or on image objects, groups of similarly identified pixels. Object-based change detection (OBCD) has been shown to have significant benefits over pixel-based change detection methods including the following: the mitigation of inherent spectral variability caused by inconsistencies in factors such as solar illumination, viewing geometry, and shadowing [44]; the reduced inclusion of vector slivers falsely identified as change caused by the delineation of raster features [45]; and enhanced landscape feature identification [46]. We chose to use a two-step pixel and OBCD methodology. In this approach, categorical change was calculated at the pixel level and individual radiometric class thresholds were applied to image objects defined as contiguous pixels with the same categorical change.

2.4. Categorical Change and Image Segmentation

A pixel-based categorical change was performed between the 1985/1995 and 1995/2010 land cover maps. Change areas were then segmented using the open-source Python package WhiteBoxTools’ clump feature [47]. This feature groups contiguous all pixels that have the same change value and assigns each clump of pixels with a unique identifier. This clump product is similar to that of a shapefile but reduces the need to convert raster data to vector format, which can be expensive in terms of computation and memory when working at landscape scales.

2.5. Change Vector Analysis

LandTrendr tools in the Google Earth Engine [33] were used to preprocess a time series of Landsat TM imagery from 1985–2010. Images were cloud-masked and mosaiced before running the temporal segmentation algorithm of LandTrendr. The advantage of using the temporally smoothed time-series is that spectral variations caused by seasonal changes are reduced while spectral variations caused by land cover changes are preserved. This process helps reduce the amount of false-positive radiometric change in the output map. Tasseled cap brightness (TCB) was used as the input spectral index required to run LandTrendr’s temporal segmentation algorithm. Additional spectral indices, including tasseled cap greenness (TCG), tasseled cap wetness (TCW), and the normalized difference vegetation index (NDVI), were normalized to the temporally segmented TCB image using LandTrendr’s fit-to-vertices tool. Previous studies note that the use of tasseled cap indices helps highlight changes in vegetation such as forest [45,48], urban areas [49], and wetlands [50].

A change vector analysis (CVA) was performed between the 1985–1995 and 1995–2010 change periods using the temporally smoothed TCB, TCG, and TCW as the spectral axes. CVA magnitude images were then downloaded from GEE for further analysis.

2.6. Hybrid Land Cover Change: Categorical and Radiometric Fusion

Categorical change clumps for both change periods were intersected with the corresponding CVA magnitude image. Mean change magnitude was calculated for each unique change clump. User-defined change magnitude thresholds for each change type were used to determine if categorical change was real. Threshold selection was done through an iterative process, where a baseline threshold was set, individual change polygons that had a mean change magnitude above the threshold were categorized as real change, results were visually inspected, and thresholds were adjusted to account for misclassifications that may have been occurring. All 132 possible change combinations were given a threshold to best represent the expected magnitude of spectral change that would occur for that given land cover change. Resulting hybrid categorical and radiometric change can be seen in Figure 3.

2.7. Validation

Land cover change validation was performed using remotely sensed information on a pixel-wise basis using TimeSync v3.0 software [51]. TimeSync is a graphical interface designed to assist image interpreters in determining if a land cover change occurred for a set of given pixels and can be used as an independent validation source to determine land cover change detection accuracy. TimeSync displays the spectral time series of a pixel that has undergone temporal segmentation (Figure 4). A total of 540 points, randomly selected from the 2010 map across the study region, were evaluated using criteria from Congalton and Green.

Validation points were stratified and weighted by land cover class to account for the area each class contributed to the total change, and were well distributed across the study area [52] (Figure 5). Plots were analyzed to validate if a change was occurring and when that change occurred.

2.8. Fusion with LandTrendr

LandTrendr was used to produce temporal change maps which note the likely year a pixel changed in each time series. If a pixel did not likely change, it was assigned a no data value. Individual temporal change maps were created using TCB, TCG, and TCW as the segmentation indices, then subsequently stacked to get a single temporal change map for both change periods of interest by taking the majority occurrence. Each spectral index (TCB, TCG, and TCW) highlights a certain change type (i.e., TCB works well for detecting urban development), and combining the three datasets into one image helps capture the diverse land cover change occurring across the GLB. The hybrid categorical/radiometric land cover change products for each change period were intersected with LandTrendr’s temporal change map product. The majority of the year of change within each clump was considered the year change occurred.

2.9. Comparison with LCMAP

For areas on the United States side of the basin, land cover change products developed for the GLB were compared with USGS’s Land Change Monitoring, Assessment, and Projection (LCMAP) [53] land cover change products as an additional validation assessment. LCMAP was developed to provide continuously updated land cover and change information for the continental United States. Finalized change products for both the 1985–1995 and 1995–2010 time periods were reclassified according to

Table 3. Reclassification of land cover product used in this study for comparison.

MTRI Reclassified	LCMAP
Urban	Developed
Suburban
Bare land	Barren
Agriculture	Agriculture
Grasslands	Grassland/shrubs
Shrubs
Deciduous	Tree cover
Evergreen
Woody wetland	Wetland
Wetland
Aquatic bed
Water	Water

and clipped to the United States only to aid in the comparison with LCMAP. Total net change for all individual land cover classes was then compared for each change time period.

3. Results

3.1. Basin Wide Change

Figure 6 displays forest area gain and loss to visualize the regional-level patterns of changes for this specific type, and Figure 7 displays basin-wide anthropogenic and natural change for two time periods. Approximately 776,854 ha of land changed from c.1980 to 1995 and approximately 998,400 ha of land changed from c.1995–2010 (

Table 4. Net land cover change gain and loss.

Land Cover	c.1980–1995 (ha)	c.1990–2010 (ha)
Agriculture gain	12,023	99,218
Agriculture loss	91,673	102,144
Aquatic bed gain	335	41,801
Aquatic bed loss	2753	3270
Barren gain	52,396	39,779
Barren loss	24,732	52,621
Deciduous gain	180,301	84,831
Deciduous loss	254,070	219,701
Evergreen gain	59,617	70,998
Evergreen loss	48,069	94,324
Grassland gain	92,468	6842
Grassland loss	96,267	159,130
Non-woody wetland gain	44,402	100,901
Non-woody wetland loss	19,459	44,904
Shrubland gain	137,143	290,795
Shrubland loss	163,203	98,202
Suburban gain	130,836	139,238
Suburban loss	6447	46,421
Urban gain	21,155	19,464
Urban loss	5655	11,837
Water gain	9091	2458
Water loss	29,272	56,869
Woody wetland gain	37,087	102,074
Woody wetland loss	35,256	108,977
Total land cover change	776,855	998,400

).

Change validation results are shown in

Table 5. Validation of the timing of land cover change detection.

Years	%
+/−0	29%
+/−1	47%
+/−2	56%
+/−3	64%
+/−4	72%
+/−5	87%
+/−6	82%

, showing the validity of cover change timing change, and in

Table 6. Summary of validation results showing accuracy of change from (a), and to (b), each land cover class.

Change From:	# Plots	Correct	Accuracy (%)		Change To:	# Plots	Correct	Accuracy (%)
Urban	7	4	57		Urban	45	43	96
Suburban	17	15	88		Suburban	168	152	90
Barren	25	21	84		Barren	22	20	91
Agriculture	115	106	92		Agriculture	25	15	60
Grassland	110	94	85		Grassland	28	27	96
Deciduous	117	104	89		Deciduous	66	58	88
Coniferous	30	29	97		Coniferous	35	34	97
Shrubs	56	53	95		Shrubs	87	80	92
Woody wetland	28	26	93		Woody wetland	27	24	89
Non-woody wetland	15	15	100		Non-woody wetland	29	26	90
Aquatic bed	4	4	100		Aquatic bed	6	6	100
Water	16	16	100		Water	2	2	100
Total	540	487			Total	540	487

, summarizing the change occurring from one class into any of the other classes (labeled “Change From”) as well as any class changing to each of the other land cover classes (labeled “Change To”). A total of 487 of the 540 validation plots were correctly identified as real land cover change. All changes from one class to another had an accuracy of 84% or higher, except changes from urban to any other class, which had an accuracy of 57%, and changes to agriculture from any other class, which had an accuracy of 60%. The accuracy of most changes was found to be above 90%.

Figure 8 shows net land cover for each time period for individual land cover classes. Both time periods displayed large net decreases in both deciduous forest, 254,069 ha and 219,700 ha, and agricultural land, 91,672 ha and 102,144 ha, changing for c.1980–1995 and c.1995–2010, respectively. Both time periods also experienced net increases in suburban cover with 137,143 ha and 290,794 ha changing for c.1980–1995 and c.1995–2010, respectively. Largest land cover gains were deciduous forest (180,300 ha), shrublands (137,143 ha), and suburban (130,836 ha) for c.1980–1995, and shrublands (290,794 ha), suburban (139,238 ha), and woody wetland (102,074 ha) for c.1995–2010. Largest land cover losses were deciduous forest (254,069 ha), shrublands (163,202 ha), and grasslands (96,267 ha) for c.1980–1995, and deciduous forest (219,700 ha), grasslands (159,129 ha), and woody wetland (108,976 ha) for c.1995–2010.

Figure 9 and Figure 10 are visual representations of land cover transitions during the c.1985–1995 and c.1995–2010 time periods, respectively. Columns represent the total areas that changed during each time period, while the flow patterns between the columns represent the transitions from each land cover into the others. Land cover patterns such as forestry cycles, where the forest is harvested, leaving barren ground or grasslands, then begins to regrow to shrublands, and finally to a forest, are well demonstrated for both time periods. Several areas across the basin are highlighted in Figure 11 to illustrate the urbanization in cities such as Detroit, MI, and Green Bay, WI, forestry in northern Ontario, and coastal wetland changes in Saginaw Bay, MI, Green Bay, WI, and Lake St. Clair.

3.2. Change by Ecoregion

Resulting LCLU products were intersected with the ecoregion map displayed in Figure 5 [54] and summarized to further observe ecological spatial distribution by ecoregion (Figure 12). The Ecoregions that were intersected are as follows: Central US Plains, the Mixed Wood Plains, and the Mixed Wood Shield. The Central USA Plains are a relatively flat region with a range of soil types, existing in the transition between the largely treeless Great Plains and the forested regions to the east. This region contains the eastern part of the “corn belt”. The Mixed Wood Plains region covers several areas from central Minnesota and Wisconsin, east through Michigan, Ontario, the New York State, and into most of New England, with more hilly topography. The Mixed Wood Shield is a region stretching east–west in the northernmost parts of Minnesota, Wisconsin, and Michigan into Canada and east to about Quebec City. This region represents the portion of the continent covered by the Canadian Shield bedrock that is dominated by mixed (hardwood–coniferous) forest.

For the 1985–1995 change period the changes observed in the central plain, mixed wood plain, and mixed wood shield make up 7% (36,693.37 ha), 30% (159,014.87 ha), and 63% (343,201.31 ha) of the total change, respectively. For the 1995–2010 change period the changes observed in the central plain, mixed wood plain, and mixed wood shield make up 7% (66,913.8 ha), 34% (313,661.12 ha), and 59% (534,924.33 ha) of the total change, respectively. Bar graphs in Figure 12 display the top five change types for each ecoregion. Many of the change types in central plains included transitions in and out of agriculture and suburban land covers; the mixed wood plain ecoregion consisted of a mixture of agricultural, suburban, and forestry changes, and all top five change types in the mixed wood shield consisted of various stages of the forestry cycle for both time periods.

3.3. Comparison to LCMAP

Change results were subset to the United States, reclassified according to Table 2, and compared to LCMAP products from matching dates (Figure 13). Overall change products showed similar trends with decreases in agriculture and grasslands/shrubs, and increases in urban areas for the 1985–1995 period and decreases in agricultural and forest and increases in urban areas for the 1995–2010 period. The largest discrepancy between the two change products occurred in the wetlands for the 1985–1995 time period, with our product showing a net decrease while the LCMAP showed a small net positive change in wetland cover.

4. Discussion

This study builds from three previous mapping endeavors [32,37,38,39] to create continuous and consistent land cover data that crosses the international border in the Great Lakes for three time periods. Relying on existing categorical data for the specific time periods of interest as the basis for the change analysis leverages the previous efforts of the three mapping sources based on their field data and the knowledge of the region from the time periods used to create the original maps. The approach contrasts with other methods of producing new maps from various time periods (e.g., [31]) by finding locations of current field data that have not changed radiometrically over time to use as training for all time periods.

Since the years of the original Canadian maps were not in line with the years of the C-CAP maps for the United States side of the GLB, we adjusted the products, working closely with the AAFC map producers to ensure the new maps were consistent. To address the change between the three time periods of interest, we incorporated both radiometric and categorical change analysis to reduce false change and used open-source tools LandTrendr [33] and TimeSync [34] for assessing time series data. The approach provides a way of assessing change across the region for all change types, providing a comprehensive dataset for use in basin-wide modeling and assessment. But there remain limitations. As with any model, the LandTrendr temporal segmentation approach that underlies the classification is an abstraction of the actual changes in the landscape. The degree to which this abstraction faithfully represents the landscape depends on the parameters and hyperparameters used to apply the algorithm; as developers of the algorithm, we made choices based on the underlying landscape types involved and the expected choices. More importantly, we recognize that our model is not perfect. Challenges include balancing signal and noise in land cover types with substantial inter-annual variability, as well as identifying real change for the first or last observations of the time series, and others. The purpose of the accuracy assessment with TimeSync is to place sideboards on the degree to which these issues introduce false change or, under capture, real change. For this study, LandTrendr was applied broadly to capture all land cover change; this caused a lack of specific change type calibration, leaving many areas denoted as change in the final change maps developed for this study without a noted change date. Because LandTrendr parameters were set to assign a default change date to all pixels, areas without a noted change date were assigned a change date either at the beginning or end of the time period. This caused a lot of the change found to be falsely assigned to either 1985, 1995, or 2010. This could be remedied by incorporating higher resolution temporal data into this research.

It should be noted that the total amount of area changed for each change period is influenced by the input classification maps, as noted in previous studies [44]. Discrepancies in the total change area between the two change periods are partially caused by the direct comparison between optical-based land cover products (C-CAP and Landsat-adjusted AAFC) [32,39] and optical/SAR fusion land cover products, developed for the 2010 baseline maps [39]. Regarding the optical/SAR approach, Bourgeau-Chavez et al. [38] notes the inclusion of SAR for the use of land cover classification, which allows for better detection of woody and non-woody wetlands over traditional optical-based land cover products due to the sensitivity of SAR to vegetation structure and inundation. We believe this more detailed land cover information in 2010, especially that of woody and non-woody wetlands, compared to that of the 1995 product, allowed for the detection of the following: (i) more land cover change in terms of quantity; (ii) more subtle land cover changes (i.e., transitions to/from woody and non-woody wetlands); and (iii) a better understanding of how urban growth may be impacting wetlands at a more detailed level. Additionally, this difference in the earlier categorical mapping methods and the 2010 product is likely the source of the discrepancies found between wetland class changes in this study vs. the LCMAP products (Figure 13b). Finally, it should be noted that direct comparison of land cover products from multiple sources can introduce false positive change into the analysis. Careful consideration of change detection methodology, as outlined in this study, is strongly recommended. Using hybrid approaches of categorical and radiometric change is one method for mitigating the problem.

As shown in Table 6, while most changes had accuracies better than 84%, change from urban to another class had an accuracy of 57%, and change from any class to agriculture had an accuracy of 60%. The low accuracy of these two cases can be explained through the unique processes related to urban and agricultural change. For the case of urban change, new architectural procedures such as new roofing materials, different pavement materials, and other infrastructure updates may show as a change when they are not actually changing classes, which causes accuracy to decrease [55]. Similarly, agriculture has a large amount of spectral variability due to differences in soil moisture and vegetation phenology depending on what has been planted at different times of the year, which complicates change accuracy [56]. In assessing the results of this study, therefore, it is important to note that within-class change, while maybe common, is not relevant for most change analysis applications.

The temporal analysis conducted revealed that aquatic bed gain was occurring for most of the time studied. Aquatic bed includes mostly submerged vegetation, the upper portions of which may float at the surface. Harmful Algal Blooms (HABs) meet this definition. Since 1985, inland water bodies within the GLB have experienced rapid warming, providing longer growth seasons and affording the cyanobacteria that causes HABs an advantage over other taxa because of their higher temperature optima range [57]. This process and increased phosphorus loads from agriculture both contribute to aquatic bed gain. Using LandTrendr to improve spectral information through regressing expected values between inflection years could improve predictions of surface water presence and HABs by removing inter-annual noise [58]. We saw similar success with the LandTrendr’s ability to remove inter-annual noise to create temporally smoothed data in calculating the change vector analysis.

We found a consistent loss in the deciduous forest area for much of the time studied. Nitrogen isotope analyses of archeological fish samples show that the nitrogenous nutrient regimes of even the world’s largest freshwater ecosystems can be highly sensitive to short-term watershed forest cover disturbances [59]. Guiry [59] found that deforestation in the 1930s led to an elevation in aquatic nitrogen sources due to the massive levels of soil denitrification caused by deforestation. This appears to be continuing, based on the levels of forest change from 1985–2010 and the expansion of the invasive plant Phragmites australis during this time period [13] (because this class is only on the 2010 map, this study was not designed to map invasive expansion). Invasive phragmites has been found to outcompete native wetland plants in high N conditions [5]. This underscores the importance of LCLU maps that allow for the quantification of change in the watersheds of the GLB, where invasive species are expanding. We anticipate that our land change maps can inform research on the role of nutrient cycling that may be affecting ecosystems basin-wide, thus allowing for the identification of potential hotspots of nutrient activity for mitigation.

Great Lakes coastal wetlands represent a vital link between land and water, providing key ecosystem services as well as habitat for a wide variety of flora and fauna. The objective of the study was to quantify basin-wide land cover changes to link these changes to the observed expansion of invasive wetland plants and to better understand wetland ecosystem dynamics through process-based modeling efforts. The creation of bi-national change maps provides a dataset that includes continuity in space and time to address this goal. The change mapping from this study is currently being used to parameterize the integration of a hydrological model (LHM, [60]) and a wetland ecosystem model (MONDRIAN [5]) to quantify potential impacts of anthropogenic changes on these very important yet vulnerable ecosystems. Furthermore, the data and results from this work will have value in supporting additional science inquiries and informing policy related to regional land use and climate change [61].

The methods used in this study allow for the continued monitoring of land cover change across the GLB as new and unique land cover classification products become available. The methods used are generic (using established analysis tools) that would have value for similar assessment in other areas with appropriate land cover datasets. This study used multi-spectral sensing (Landsat), which provides a long-time record of landscapes, for radiometric change analysis. SAR, advanced multi-spectral, and imaging spectroscopy sensing systems are being planned, providing new data sources for land monitoring and change studies. As more spaceborne SAR sensors become available, detailed land cover classifications that discriminate wetland classes will improve the monitoring of dynamic wetlands. Similarly, the spectral diversity of new optical sensors, such as those planned for LandsatNEXT [62] and the Surface Biology and Geology (SBG) [63] mission, will allow for the discrimination of agricultural systems, such as crop and harvest patterns. Future studies on how to exploit these capabilities are needed but are optimistic with the increasing availability of small satellites and data analytical resources, such as cloud computing and machine learning, that can efficiently process large datasets.

5. Conclusions

In this study, we demonstrate a methodology for assessing land cover change at specific intervals of interest that allows for the analysis of major ecological events. Studying land cover c.1985 allows for a dataset that was collected prior to the water level decline and influx of wetland plant invasion. Land cover c. 1995 captures mid-invasion conditions, and land cover c. 2010 captures a low-water-level period when the invasive species reached the most expansive area. Through this research, it was possible to create change maps for the time periods of interest for the United States and Canada with high accuracy.

The application of change mapping across full watersheds of Lakes Erie, Huron and Michigan is important because, although significant relationships between land cover (% Natural and % Urban) and water quality are observable at several scales, Croft-White et al. [27] found relationships to be strongest at the watershed scale. We provide a continuous land cover change product for whole watersheds of the GLB that highlight significant land cover classes useful for ecological and hydrological modeling, such as the distinction between non-woody and woody wetlands, and deciduous and coniferous forest. The methodology described in this study demonstrated its effectiveness in capturing the varying spectral characteristics of a large array of land cover changes. This approach provides a rigorous set of data valuable for many applications and policy needs using methods transferable to anywhere where appropriate land cover data are available.