Characterizing the Landscape Structure of Urban Wetlands Using Terrain and Landscape Indices

: Several studies have shown human impacts on urban wetlands. These impacts are mostly studied at broad scales, which may generalize and aggregate important information needed for landscape quantiﬁcation or terrain analysis. This situation can weakly or inappropriately address the structure of wetland landscapes, thus a ﬀ ecting the assessment of the quantities and qualities of terrestrial wetland habitats. To address these issues for urban wetland dynamics, this study proposes the use of landscape and terrain indices to characterize the landscape structure of urban wetlands at a ﬁne scale in order to assess its usefulness in contributing to wildlife sustainability. To achieve this goal, secondary terrain attribute data are integrated with an object-based satellite image classiﬁcation at the wetland and watershed level. The result reveals a general swell in wetland coverage at the watershed level. Further analysis shows the size and shape complexities, and edge irregularities are increased signiﬁcantly at the patch level but slightly at the watershed level. Terrain analysis further reveals a potential increase in wetness and decrease in stream power vulnerability for most of the major wetlands under study. These results suggest that terrain and landscape indices are e ﬀ ective in characterizing the structure of urban wetlands that supports socio-ecological sustainability.


Introduction
Urban wetlands are an important part of the global ecological system, which have been affected by various human activities, mostly urban land developments [1]. Land-use change has been recognized as one of the biggest anthropogenic influences on global ecosystems [2], including by altering landscapes at various scales. The altering of landscape structures not only affects the local habitat and the quality of organisms, but also the composition and configuration of the surrounding landscape [3]. As a result, the resulting landscape patterning has largely determined the quality of resources, structure, and function of many urban ecosystems, including most urban wetlands [4][5][6]. To quantify these land cover changes, mapping and monitoring of urban wetland landscape dynamics within the terrestrial habitats are needed. Remote sensing-based techniques have been proven to be a valuable tool in this regard [7][8][9]. While many of these studies have produced interesting results on urban land cover dynamics, the landscape structure of urban terrestrial wetlands resulting from land-use and land-cover (LULC) changes have not been fully studied. We believe that understanding human impacts on the structure of urban landscapes, which could result in fragmentation of habitat loss, is essential for the sustainability of the urban wetland ecosystem. The concept of landscape structure has been defined differently in different studies [10,11]. However, regardless of how it has been described, assess its usefulness in contributing to wildlife sustainability.

Study Area
The three major watersheds in the Kansas City metropolitan area and ten of its major wetlands were included in this study. The Kansas City metropolitan area is located along the eastern boundary of Kansas and the western boundary of Missouri in the central United States [35]. The general topography of the area is characterized by rolling hills with open plains [8]. The predominant land-cover types are grasslands, forests, and croplands [35]. Figure 1 shows the three major watersheds in the Kansas City area, which include the Blue River, Little Blue River, and Shoal Creek-Missouri River watersheds. Ten major wetland areas in the three major watersheds were also included ( Figure 1). These wetlands are the Missouri River, Lake Tapawingo, Blue Springs Reservoir, East Lake Wood, West Lake Wood, Lake Jacomo, Prairie Lee Lake, Longview Lake, Heritage Park Lake, and Loch Lloyd Lake.
In recent decades, population growth, increased economic activities, and the continued expansion of urbanized areas in the watersheds and wetland areas have affected the natural environment of the area [6,8,36]. This continued depletion of the natural area provided a basis for investigating the terrestrial wetland landscape of our study area. While the landscape dynamics and trends of the area have been well studied, with examples in references [6,8], less has been done to identify and quantify the landscape structure of the terrestrial habitat surrounding the urban wetland landscape.  In recent decades, population growth, increased economic activities, and the continued expansion of urbanized areas in the watersheds and wetland areas have affected the natural environment of the area [6,8,36]. This continued depletion of the natural area provided a basis for investigating the terrestrial wetland landscape of our study area. While the landscape dynamics and trends of the area have been well studied, with examples in references [6,8], less has been done to identify and quantify the landscape structure of the terrestrial habitat surrounding the urban wetland landscape.

Methods
The methodology used in this study is divided into four parts (see Figure 2). The first part employed an object-based image analysis (OBIA) approach for image classification. Two classification algorithms: support vector machine (SVM) and K-Nearest Neighbor (K-NN) were used. Both methods involved Land 2020, 9,29 4 of 25 feature extraction using an example-based method for segmentation and merging. Furthermore, Change detection statistics (CDS) was performed on the classified images resulting from the application of both algorithms. This was done to understand the class change among the LULC in different classes. CDS was used specifically for the calculation of wetland coverage change at the watershed level in this study. The second part involved deriving secondary terrain attributes using the compound topographic index (CTI) and stream power index (SPI) from the United States Geological Survey (USGS) digital elevation model (DEM). These indices were used for wetness and stream power quantification around the terrestrial wetlands. The third part involved metric calculations, which were implemented based on the spatial patterning the structure of landscape heterogeneity in relation to some aspect of ecological function (e.g., a direct or indirect role). The major metrics employed for analysis were related to the edge, diversity, shape, and size of wetland patches. Fourthly, in the methodology, a 30-m inward buffer was performed on the extracted wetlands to serve as estimated core areas. These were then clipped to derive secondary terrain attributes (CTI and SPI) at 340-m buffers, which served as the terrestrial wetland habitat area (see Section 2.2.5).

Methods
The methodology used in this study is divided into four parts (see Figure 2). The first part employed an object-based image analysis (OBIA) approach for image classification. Two classification algorithms: support vector machine (SVM) and K-Nearest Neighbor (K-NN) were used. Both methods involved feature extraction using an example-based method for segmentation and merging. Furthermore, Change detection statistics (CDS) was performed on the classified images resulting from the application of both algorithms. This was done to understand the class change among the LULC in different classes. CDS was used specifically for the calculation of wetland coverage change at the watershed level in this study. The second part involved deriving secondary terrain attributes using the compound topographic index (CTI) and stream power index (SPI) from the United States Geological Survey (USGS) digital elevation model (DEM). These indices were used for wetness and stream power quantification around the terrestrial wetlands. The third part involved metric calculations, which were implemented based on the spatial patterning the structure of landscape heterogeneity in relation to some aspect of ecological function (e.g., a direct or indirect role). The major metrics employed for analysis were related to the edge, diversity, shape, and size of wetland patches. Fourthly, in the methodology, a 30-m inward buffer was performed on the extracted wetlands to serve as estimated core areas. These were then clipped to derive secondary terrain attributes (CTI and SPI) at 340-m buffers, which served as the terrestrial wetland habitat area (see Section 2.2.5).

Image Classification techniques
The classification techniques employed for the study involved testing two algorithms to determine which of the two provides better accuracy ( Figure 3). The spectral analysis was performed in ENVI before the image extraction and classification. SPOT 2 sensors with the data collected in three relative broad multispectral bands and one panchromatic band were compared to SPOT 7 that contains four similar bands. In this study, the two images were compared using the spectral signature to distinguish the types of ground cover or objects performed in ENVI software. The processing level 1A, which involves correction by normalizing the charge-coupled device (CCD) response to compensate for radiometric variations due to detector sensitivity, was applied to SPOT 2 image 1992. On the other hand, the 2017 SPOT 7 image was pre-processed on the 2A level, meaning scenes were rectified to match a standard map projection (UTM WGS 84) without using ground control points. Level 2A is the entry-level map product for the SPOT image. This means no geometric correction was needed for SPOT 7, as the systematic distortion effects and transformations were compensated for by projecting the image in a standard map projection UTM WGS 84.

Image Classification Techniques
The classification techniques employed for the study involved testing two algorithms to determine which of the two provides better accuracy ( Figure 3). The spectral analysis was performed in ENVI before the image extraction and classification. SPOT 2 sensors with the data collected in three relative broad multispectral bands and one panchromatic band were compared to SPOT 7 that contains four similar bands. In this study, the two images were compared using the spectral signature to distinguish the types of ground cover or objects performed in ENVI software. The processing level 1A, which involves correction by normalizing the charge-coupled device (CCD) response to compensate for radiometric variations due to detector sensitivity, was applied to SPOT 2 image 1992. On the other hand, the 2017 SPOT 7 image was pre-processed on the 2A level, meaning scenes were rectified to match a standard map projection (UTM WGS 84) without using ground control points. Level 2A is the entry-level map product for the SPOT image. This means no geometric correction was needed for SPOT 7, as the systematic distortion effects and transformations were compensated for by projecting the image in a standard map projection UTM WGS 84.
The OBIA was then applied, which involved image segmentation and merging processes being performed on two archived SPOT imageries. Both SPOT 2 (29 January 1992; 20 m × 20 m resolution) and SPOT 7 (22 October 2017; 6 m × 6 m resolution) were classified using K-Nearest Neighbor (K-NN) and Support Vector Machine (SVM) algorithms. Both datasets were resampled to a uniform pixel size, Land 2020, 9, 29 5 of 25 with the adoption of a representative minimum mapping unit (6 × 6 m). Jensen [37] reiterated that this does not have much of an impact on the results, because the resampled data can never be greater than the instantaneous field of view (IFOV). He mentioned that additional information should not be expected from a resampled image. In addition, Jensen suggested that object-based classification may be a good alternative to the traditional pixel-based methods. This is to overcome the high-resolution problem and the salt and pepper effect. The object segmentation considers group of pixels characterized and organized into objects and treats each object as a minimum classification unit [38]. The OBIA was then applied, which involved image segmentation and merging processes being performed on two archived SPOT imageries. Both SPOT 2 (29 January 1992; 20 m × 20 m resolution) and SPOT 7 (22 October 2017; 6 m x 6 m resolution) were classified using K-Nearest Neighbor (K-NN) and Support Vector Machine (SVM) algorithms. Both datasets were resampled to a uniform pixel size, with the adoption of a representative minimum mapping unit (6 x 6 m). Jensen [37] reiterated that this does not have much of an impact on the results, because the resampled data can never be greater than the instantaneous field of view (IFOV). He mentioned that additional information should not be expected from a resampled image. In addition, Jensen suggested that object-based classification may be a good alternative to the traditional pixel-based methods. This is to overcome the high-resolution problem and the salt and pepper effect. The object segmentation considers group of pixels characterized and organized into objects and treats each object as a minimum classification unit [38].
The images were selected with the consideration of their availability and the type of classification methods used to ensure that they will not affect the accuracy of the result. O'Neill et al. [39] proposed that to avoid bias in calculating landscape metrics, the grain size of the map should be 2-5 times smaller than the spatial features being analyzed and the map extent should be 2-5 times larger than the largest patch. In addition, the classification method adopted for segmentation and merging focused more on the spatial, spectral and textual information of the images than alternative methods. Segmentation is the process of partitioning an image into objects by grouping neighboring pixels with common values [40]. Merging, on the other hand, combines the adjacent segments with similar spectral attributes textural information together based on image resolution [40]. This process is example-based, which is also referred to as supervised classification. It is used to obtain training data as objects and then assigns them to one or more known features. Archived SPOT imageries with The images were selected with the consideration of their availability and the type of classification methods used to ensure that they will not affect the accuracy of the result. O'Neill et al. [39] proposed that to avoid bias in calculating landscape metrics, the grain size of the map should be 2-5 times smaller than the spatial features being analyzed and the map extent should be 2-5 times larger than the largest patch. In addition, the classification method adopted for segmentation and merging focused more on the spatial, spectral and textual information of the images than alternative methods. Segmentation is the process of partitioning an image into objects by grouping neighboring pixels with common values [40]. Merging, on the other hand, combines the adjacent segments with similar spectral attributes textural information together based on image resolution [40]. This process is example-based, which is also referred to as supervised classification. It is used to obtain training data as objects and then assigns them to one or more known features. Archived SPOT imageries with a temporal resolution of 25 years between the two images were classified. Overall, a segmentation process using feature extraction and machine learning techniques were sequentially applied in the classification process. For most of the analysis in the study, the SVM classification result was used because it provided higher accuracy than the other methods (see Figure 3).

Classification Scheme
This study adopted a classification scheme that grouped the different classes on the images into four (Table 1). This is mostly due to the issue of spectral signature similarity experienced during the process of image segmentation when generating the training sites. The four categories as shown in Table 1 include Impervious Surfaces (IS) which represent urbanized areas; Forestland (FL), which covers all areas with a collection of trees. Further, all areas with brush, crops and non-forest vegetation/grassland were classified as Farmland/grassland (FGL), and Open water bodies were classified as Wetlands (WL). The output was a classified raster grid, which was then converted to vector polygons to perform landscape metrics calculation on them. Patch Analyst 5, an extension of ArcGIS 10.6 software, was used for the landscape metrics calculations. ENVI software confusion matrix was used to calculate the accuracy of the Object-Oriented Classification, with emphasis on wetland classification. Fifty regions of interest (ROI) polygons for each class were selected and referenced using Google earth archival imageries and other classified maps of the Kansas City area. The report of the calculation pairs ROIs with the classes of classified images to show what percentage of the ROI pixels that were or were not contained in a resulting class. The overall accuracy was calculated by summing the number of correctly classified values and dividing by the total number of values. Table 2 shows the accuracy assessment reports with a higher overall accuracy percentage for SVM classification for both SPOT 1992 and 2017 images and the overall accuracy for K-NN SPOT 1992 and 2017. SVM-based classification for SPOT 2017 revealed an 89.14% accuracy while the K-NN-based classification on the same image resulted in a 79.54% total accuracy. The SVM-based classification for SPOT 1992 resulted in 63.84% while the K-NN approach for the same year gave an overall accuracy of 61.42%. Overall, little disparity was revealed for both algorithms when results for each year were compared. However, the SVM accuracy results were better when compared to that of K-NN ( Table 2). Despite that two different classification algorithms were used on the same image, their accuracy reports revealed a consistent trend. Similar trends in landscape patterns were obvious using the two classification algorithms (see Table 2). Table 2. Accuracy assessment with emphasis on wetland accuracy from both classifications. Accuracy Assessment for both Classification algorithms: support vector machine (SVM) and K-Nearest Neighbor (K-NN).

Terrain Analysis
Terrain analysis was performed with the raster DEM data of 1/9 arc-second acquired from the USGS National Elevation Dataset (NED), which has an approximately 10-m horizontal resolution. These datasets from around 2008 through 2013 became the seamless DEM layer under the 3DEP Land 2020, 9, 29 7 of 25 program [41]. The NED is derived from diverse sources of data and were later processed into a common coordinate system and unit of vertical measure. The geographic coordinates system (GCS) is the North American Datum of 1983 (NAD83). The elevation values are provided in units of meters and referenced to in North American Vertical Datum of 1988 (NAVD88) over the conterminous United States. The pre-processing of the DEM was performed in ArcGIS 10.6 software, and included pit and sink filling. This is a procedure to fill depressions with hypothetical water flows and force drainage to the lowest possible outlets [42]. Also, primary terrain attributes were calculated, which included slope, flow direction, and flow accumulation. Secondary terrain attributes were derived from the primary terrain attributes for SPI and CTI using the raster calculator in ArcGIS 10.6 software. SPI is calculated as the product of the natural log of both slope and flow accumulation, while CTI is the quotient of both slope and flow accumulation [43]. The former measures the power of streamflow while the latter measures potential wetness locations. where: Ln represents natural logarithm; A represents the catchment area per pixel; β refers to the slope in degrees.
The two secondary terrain indices for wetness and streamflow were adopted based on Wilson and Gallant [43].

Urban Wetland Terrestrial Habitat Buffer
To understand the landscape structure surrounding the terrestrial wetland habitat necessary for the urban wetland ecosystem, a buffer of 340-m from the core area was adapted for the study [44,45]. To achieve this, a 30-m inward buffer was used on the extracted wetland files. This served as the aquatic buffer zone from the core wetland. Another 260-m outward buffer was created from the extracted wetland files, which served as the core habitat area. The aquatic buffer is also part of the core habitat area; an additional 50-m buffer recommended by Murcia [44] to protect core habitat from edge effects was used (see Figure 4a,b). This implementation is similar to a study by the Environmental Law Institute (2008), in which the area served as a terrestrial wetland area sufficient for the movement of wetland wildlife (e.g., Amphibians and Reptiles) through the landscape. The buffered terrestrial wetland area files were clipped to derive secondary terrain attributes, CTI and SPI. The clipped portion was reclassified into wetness and non-wetness areas for CTI, and stream power and non-stream power areas for SPI using the reclassify tool in the ArcGIS 10.6 software. This was done for the ten major wetlands in the study area. The percentages of vulnerability are the potential area of gully erosion for CTI, and potential area of increased stream power for SPI. These were calculated using the pixels spilling into the different segments of the reclassified and clipped terrestrial wetland habitats. Following this process, the study estimated the percentage vulnerability for both years studied and the results compared.
The change in vulnerability estimate for CTI The change in vulnerability estimate for SPI %∆vE = y 1 Land 2020, 9, 29 8 of 25 %∆vE-percentage change in vulnerability estimate; χ 1 -wetness area for CTI; y 1 -stream power area for SPI; χ 0 -non-wetness area for CTI; y 0 -non-stream power area for SPI.  . Buffered protection zones: wetlands (a) and streams (b). Source: Adapted for this study based on Semlitsch and Bodie [45], with an additional 50-m buffer to protect the core habitat from edge effects recommended by Murcia [44].

Landscape Metrics Calculations
Landscape metrics quantification was performed at two major levels and loosely grouped according to the heterogeneity and spatial pattern of urban wetlands in the study area. These two levels are landscape level (three major watersheds), and patch level (individual wetland core area). Spatial statistics in Patch Analyst 5, an extension to the ArcGIS 10.6 software for spatial analysis of landscape patches, and modeling was used for the analysis. The default map unit for area is in meters (m) and was later stated in hectares (ha), which is an option provided by the Patch Analyst tool. The choice of these metrics (Table 3) was based on the fact that most landscape metric results are correlated and statistically redundant. That is, they may quantify a similar or identical aspect of landscape patterns [10,11]. However, some metrics may be empirically redundant not because they measure same aspect of the landscape pattern, but they are alternative ways of representing similar information for the landscape under investigation [10]. The metrics used in this study took into consideration the scale of change patterns of the variables under study. This was done with specific consideration to the spatial pattern of the entire landscape and the patch types, bearing in mind the heterogeneity of the terrestrial wetland area. The major components (extent, subdivision, geometry, isolation, and connectedness) of habitat landscape composition and configuration were also considered when interpreting the landscape structure of our study area. For the three major watersheds, "landscape-centric" metrics were considered for measurement at this level. These are metrics that summarized patches using the mean and the area-weighted mean [10]. Also, for the individual patches such as wetlands, "patch-centric" quantifications were performed which measured or described the spatial context of an isolated patch [10]. Landscape metrics quantification was performed on the derived watershed polygons for the three major watersheds. This study followed the recommendation of reference [45] to preserve wetland core habitats such as amphibians and reptiles. Table 3 depicts the metrics used in this study, including acronyms, descriptions, and justification.  [45], with an additional 50-m buffer to protect the core habitat from edge effects recommended by Murcia [44].

Landscape Metrics Calculations
Landscape metrics quantification was performed at two major levels and loosely grouped according to the heterogeneity and spatial pattern of urban wetlands in the study area. These two levels are landscape level (three major watersheds), and patch level (individual wetland core area). Spatial statistics in Patch Analyst 5, an extension to the ArcGIS 10.6 software for spatial analysis of landscape patches, and modeling was used for the analysis. The default map unit for area is in meters (m) and was later stated in hectares (ha), which is an option provided by the Patch Analyst tool. The choice of these metrics (Table 3) was based on the fact that most landscape metric results are correlated and statistically redundant. That is, they may quantify a similar or identical aspect of landscape patterns [10,11]. However, some metrics may be empirically redundant not because they measure same aspect of the landscape pattern, but they are alternative ways of representing similar information for the landscape under investigation [10]. The metrics used in this study took into consideration the scale of change patterns of the variables under study. This was done with specific consideration to the spatial pattern of the entire landscape and the patch types, bearing in mind the heterogeneity of the terrestrial wetland area. The major components (extent, subdivision, geometry, isolation, and connectedness) of habitat landscape composition and configuration were also considered when interpreting the landscape structure of our study area. For the three major watersheds, "landscape-centric" metrics were considered for measurement at this level. These are metrics that summarized patches using the mean and the area-weighted mean [10]. Also, for the individual patches such as wetlands, "patch-centric" quantifications were performed which measured or described the spatial context of an isolated patch [10]. Landscape metrics quantification was performed on the derived watershed polygons for the three major watersheds. This study followed the recommendation of reference [45] to preserve wetland core habitats such as amphibians and reptiles. Table 3 depicts the metrics used in this study, including acronyms, descriptions, and justification.

Results
The analysis was based on the wetland (patch) and watershed (landscape) quantification in terms of the level of heterogeneity of the study area. At the watershed level, the landscape level analysis in this study included change detection statistics (CDS) analysis and landscape metrics. The CDS was done to get a clearer view of the overall change in wetland coverage in the three major watersheds. The landscape metric quantification was performed for the three major watersheds. The other part of the analysis was at the wetland level that describes the patches. These include calculating and deriving secondary terrain attributes and integration with extracted wetland files on one hand. On the other hand, the quantification of patch metrics was performed to assess the ten major wetlands in the three major watersheds in the study area.

Change Detection Statistics (CDS)
At the landscape level, the CDS was used to measure the changes between a pair of images that present the initial state (SPOT 1992) and the final state (SPOT 2017). This was performed for three watersheds to quantify the total wetlands coverage between 1992 and 2017 in the study area. The CDS results will guide the interpretation of the metrics and terrain analysis, as it provides prior knowledge of the general wetland level. The CDS as summarized in Table 4 for SVM algorithms classification reveals an increase in wetland coverage from 1992 to 2017. The result showed an increase of 21% class change from the initial state (1992) to the final state (2017) for the three watersheds. The CDS in Table 5 Land 2020, 9, 29 10 of 25 for K-NN algorithms revealed an 18% class change. This result is similar to the change observed for SVM algorithm classification for the three watersheds between 1992 and 2017.

Landscape Level Metric Calculation
The result of the metrics calculation examined at the landscape level for watersheds is shown in Figure 5. For the landscape metric calculation, at the watershed level, the three major watersheds in the Kansas City area were quantified. Indices used were SDI, SEI, AWMSI, MSI, MPFD, ED, and MPS. The identified diversity indices used are SDI and SEI, which measure at the landscape level. For the shape metrics MSI, we used MPFD, ED, and AWMSI to quantify the shape complexity of the watershed at the landscape level. Further, MPS was used to quantify the size irregularities between the study periods.

Terrain calculation
The analysis at this level deals with individual wetland patches; ten major wetlands were calculated in the study area using digital terrain analysis. The terrain analysis involves integrating extracted wetland files, buffered at a distance of 340-m from the core area with secondary terrain attributes. The results are summarized in Figure 6a-d. All the major wetlands except the Heritage Park Lake showed a slight reduction in the core area for 2017. Similarly, nine out of the ten major wetlands examined for wetness vulnerability using the CTI reveal potential wetness in wetlands within the study area in 2017 (see Appendix A: Table A1). The only wetland that did not show a slight wetness potential is Heritage Park Lake, with -0.16%. On the other hand, the SPI calculated revealed a relatively no change in stream power, except for three averagely small wetlands. The three wetlands slightly increased with a net change of 3.31% for East Lake Wood, 2.97% for Lake Jacomo, and 2.71% for Prairie Lee Lake in 2017 (see Appendix B: Table A2).

Terrain Calculation
The analysis at this level deals with individual wetland patches; ten major wetlands were calculated in the study area using digital terrain analysis. The terrain analysis involves integrating extracted wetland files, buffered at a distance of 340-m from the core area with secondary terrain attributes. The results are summarized in Figure 6a-d. All the major wetlands except the Heritage Park Lake showed a slight reduction in the core area for 2017. Similarly, nine out of the ten major wetlands examined for wetness vulnerability using the CTI reveal potential wetness in wetlands within the study area in 2017 (see Appendix A: Table A1). The only wetland that did not show a slight wetness potential is Heritage Park Lake, with -0.16%. On the other hand, the SPI calculated revealed a relatively no change in stream power, except for three averagely small wetlands. The three wetlands slightly increased with a net change of 3.31% for East Lake Wood, 2.97% for Lake Jacomo, and 2.71% for Prairie Lee Lake in 2017 (see Appendix B: Table A2).

Terrain calculation
The analysis at this level deals with individual wetland patches; ten major wetlands were calculated in the study area using digital terrain analysis. The terrain analysis involves integrating extracted wetland files, buffered at a distance of 340-m from the core area with secondary terrain attributes. The results are summarized in Figure 6a-d. All the major wetlands except the Heritage Park Lake showed a slight reduction in the core area for 2017. Similarly, nine out of the ten major wetlands examined for wetness vulnerability using the CTI reveal potential wetness in wetlands within the study area in 2017 (see Appendix A: Table A1). The only wetland that did not show a slight wetness potential is Heritage Park Lake, with -0.16%. On the other hand, the SPI calculated revealed a relatively no change in stream power, except for three averagely small wetlands. The three wetlands slightly increased with a net change of 3.31% for East Lake Wood, 2.97% for Lake Jacomo, and 2.71% for Prairie Lee Lake in 2017 (see Appendix B: Table A2).

Patch Level Metric Calculation
Similar to the terrain analysis, the patch metrics were calculated for the ten major wetlands in the study area. Patch-centric metrics were used to quantify the patch in terms of shapes and sizes. Indices used include MPE, CA, MSI, SI, and TE. As summarized in Figure 7a-e, CA reveals a general decrease for all ten wetlands for 2017. Heritage Park Lake and Loch Lloyd Lake revealed the least decrease in terms of core area change for 2017. Heritage Park Lake showed 108666.04 ha for 1992 and 103089.57 ha for 2017, with a net loss of 5576.467 ha between 1992 and 2017. Similarly, the Loch Lloyd Lake showed 287215.62 ha for 1992 and 271206.34 ha for 2017, with a net loss of 16009.28 ha. Blue Springs Reservoir revealed the most decreased core area change for 2017 with 3107875.70 ha for 1992 and 2338980.60 ha for 2017, a net loss of 768895.01 ha. On the other hand, the SI and MSI indices used to quantify the wetland shape revealed increased shape complexity and irregularities for all the ten wetlands in 2017. This was not the case for the TE and MPE used to quantify the wetland edge. While nine of the ten wetlands show an increase in the edge difference, Blue Springs Reservoir did not increase edge difference for 2017, as compared to 1992.

Patch level metric calculation
Similar to the terrain analysis, the patch metrics were calculated for the ten major wetlands in the study area. Patch-centric metrics were used to quantify the patch in terms of shapes and sizes. Indices used include MPE, CA, MSI, SI, and TE. As summarized in Figure 7a-e, CA reveals a general decrease for all ten wetlands for 2017. Heritage Park Lake and Loch Lloyd Lake revealed the least decrease in terms of core area change for 2017. Heritage Park Lake showed 108666.04 ha for 1992 and 103089.57 ha for 2017, with a net loss of 5576.467 ha between 1992 and 2017. Similarly, the Loch Lloyd Lake showed 287215.62 ha for 1992 and 271206.34 ha for 2017, with a net loss of 16009.28 ha. Blue Springs Reservoir revealed the most decreased core area change for 2017 with 3107875.70 ha for 1992 and 2338980.60 ha for 2017, a net loss of 768895.01 ha. On the other hand, the SI and MSI indices used to quantify the wetland shape revealed increased shape complexity and irregularities for all the ten wetlands in 2017. This was not the case for the TE and MPE used to quantify the wetland edge. While nine of the ten wetlands show an increase in the edge difference, Blue Springs Reservoir did not increase edge difference for 2017, as compared to 1992.

Discussion
Firstly, the study revealed a similar trend for class change for the two algorithms used, with SVM having a better accuracy assessment report compared to K-NN. The CDS was performed with two algorithms and the change detection revealed K-NN with a value of 21% and SVM with a value of 18%. This indicates that the total percentage change in pixels from other classes to wetland increased for K-NN by 21% but not as much for SVM at 18%. The two algorithms were compared to achieve a better result for the image classification, which avoids the possible misinterpretations that could have resulted from using different multispectral images [46]. Regardless of the class change in the CDS, results from using both algorithms revealed a swell in the wetland coverage for the study area. This swell may be a result of increased precipitation in the past decade [6], and improved streamside ordinance protection for wetlands [1] in the study area. In addition, this result is similar to a study by Ji [1] which showed that larger wetlands accumulated more precipitation, while the smaller wetlands were prone to human impacts. Similarly, Zubair et al. [6] in their study of urban wetland modeling in the study area observed that wetlands increased in two of the major watersheds historically, but reduced as a result of urban expansion in one. These results generally support our findings that wetlands swelled in the study area within our study period, and this can be associated with the effects of human and natural factors as suggested by previous studies, such as reference [1].
Secondly, the resulting indices for terrain analysis that showed an increase in potential wetness for nine out of ten wetlands studied may indicate the level of activities around these wetlands. This may be associated with the size of the wetland core area. Heritage Park Lake is one of the smallest of the ten major wetlands, and the effect of natural and human activities may have more influence on this wetland. This is similar to the study by Tomaselli et al. [32] where the small landscape that could harbor biodiversity were highly fragmented in small patches. On the other hand, wetlands with a slight increase for CTI can be an indirect influence on soil moisture and provide information about the soil condition near the wetland. O'Neill et al. [39] in a similar study revealed how the topsoil is lost through erosion because of fragmentation due to distance between patches. This can also be an indirect indication of changes in the type of vegetation surrounding the terrestrial wetland areas over the years. Furthermore, Weilert et al. [9] when studying the effects of a streamside ordinance protecting the riparian vegetation within the ordinance area, found that the streamside ordinance can effectively reduce human impacts within the protected area. Most vegetated areas are close to the wetlands, and the availability of moisture in a protected wetland can improve the growth rate of trees and reduce erosive power. The condition of the wetland terrestrial soil and vegetation can impact on biodiversity, especially for reptiles and amphibians. The slight stream power increase for the averagely small wetlands can result in vulnerability of the terrestrial wetlands. These impacts may indirectly increase the erosive potential for the terrestrial wetland and expose the locations to severe gullies [13]. In addition, the results of the topographic indices can be used to assess the vulnerability of the wetland habitat condition available for biodiversity sustainability, as shown in the studies conducted by Murcia [44], and Semlitsch and Bodie [45].
Thirdly, the landscape structure was also quantified at the patch and landscape level. It is important to interpret each metric in a manner appropriate to its scale [1,14,39]. At the patch levels, which describe wetlands in this study, the indices used for the metric quantification are MPE, CA, MSI, SI, and TE. In this study, decreases in the core area for all the ten wetlands were revealed for the studied period between 1992 and 2017. According to McGarigal [47], all other things held constant, increasing shape complexity decreases the core area, and increasing patch area increases the core area. This decrease in the core area implies that wetlands during the period of study may have an increased patch area. The complexities and irregularities in the core area may be a result of human activities (e.g., agriculture etc.) around most of the terrestrial wetland areas [6,9]. The revealing effect of MSI and SI on shape increased for 2017. This implies a decrease in the core area. In addition, the increase edge effect revealed by the TE and MPE may also result in a decreased core area. The edge effect differs among organisms at different ecological habitats [48]. This implies that an increase in edge effects may impact the wildlife and biodiversity around the terrestrial wetland areas, particularly the amphibians and reptiles. Also, core areas are a much better predictor of habitat quality than patch areas [49]. Some human-induced activities and can result in habitat loss and may have a considerable impact on the terrestrial habitat. These effects are bound to affect most species with low mobility, a narrow feeding niche, and low reproduction [50] around the terrestrial habitat. These effects can be very useful to assess when monitoring the landscape structure surrounding wetland and riparian habitats. They provided a good link between LULC and ecological modeling for predictions or informed management. Fourthly, at the landscape level quantification, the diversity indices used SDI and SEI revealed little or no proportional diversity within the study period. Similarly, the shape indices MSI, MPFD, and AWMSI revealed slight shape complexity at the watershed level. An increase in edge density (ED) in 2017 relative to the landscape area in 1992 was exhibited. This increased edge effect may reflect a reduced core area resulting from more activities around the watershed in 2017 compared to 1992. This indicates that some part of the wetland landscape areas might have been converted to other usages as a result of urbanization, similar to in a study conducted by Liu et al. [51]. Generally, the landscape level quantification revealed changes in ED and a slight increase in shape complexities and irregularities. This might have resulted from human activities in the form of urbanization or agriculture impacting on the area during the period of study. This is similar to another study by Zubair et al. [36] revealing the depletion of wetland landscape at the sub-watershed level in this area due to agriculture activities. This shows that response of species richness for agriculturally induced fragmented wetlands for wildlife such as amphibians and reptiles can be monitored, even at a very large-scale map (small area), similar to the study by Kolozsvary and Swihart [52].

Conclusions
Firstly, the general status of the coverage of the wetlands in the three major watersheds was revealed statistically, using an object-based classification method. Secondly, the study has been able to reveal some of the dynamics of urban wetlands as it relates to the terrestrial portion of the landscape. The impact has been found to be more on the smaller wetlands as compared to the larger ones, which could influence some of the wildlife. Thirdly, the impact on wetland core area in terms of potential wetness and stream power was observed as it relates to terrestrial habitat. In addition, the impacted areas that could affect wildlife such as amphibians and reptiles showed increased edge effects and shape complexities, particularly at the patch level.
At the patch level, the analysis is effective in revealing the fine-scale structure of urban wetlands. The general change pattern of wetland coverage can be better observed individually at the patch level than at the landscape level. As shown in the result of the terrain analysis, most major wetlands in the study area experienced an increase in wetness potential, which may have contributed to the swell of wetland at the landscape level. The dynamics of the wetlands were better depicted at the patch level using the terrain attribute data. In addition, the effects on size, shape, and edge that were generalized at the landscape level could be better identified at the patch level. Increased edge effects and shape complexity were more obvious at the landscape level, and this may be associated with various urban development activities in the study area. Though it may not be practical when monitoring dynamic changes, acquiring remote sensing data with a similar anniversary date will assist in improving the change detection results. In addition, the ability to capture landscape functionality and spatial heterogeneity will further improve future studies.