WETLAND MAPPING IN NEW BRUNSWICK, CANADA WITH LANDSAT5-TM, ALOS-PALSAR, AND RADARSAT-2 IMAGERY

Several maps of wetland areas in central New Brunswick, Canada, were produced by applying the Random Forests classifier to different combinations of optical Landsat-5 TM images, dual-polarized (HH, HV) Radarsat-2 C-band and Alos-1 PalSAR L-band Synthetic Aperture Radar (SAR) images and digital elevation data. The resulting maps were compared to 199 GPS wetland sites that were visited between 2012 and 2018 as well as to a combination of two wetland maps currently used by the Province of New Brunswick. The number of correctly identified GPS wetland sites was the highest when both the Alos-PalSAR and Radarsat-2 images are used (97.9%). This percentage of correctly identified sites were well above the accuracy of the official New Brunswick wetland maps (44.7 %). With the best-classified image, the misidentifications were due to wetlands not being classified in the right wetland class, and just one case was a wetland site being classified in a non-wetland class. For the NB wetland map, about a quarter of the wetland validation sites were classified in a non-wetland class, and about the same number of sites were classified in the wrong wetland class.


INTRODUCTION
Wetlands are important ecosystems that perform a variety of services that are beneficial to society and the environment. They are crucial for groundwater recharge, flood control, water quality improvement, and mitigation of erosion (Li, Chen, 2005). Understanding the distribution and dynamics of wetlands is essential for understanding ecosystem diversity and function and how human practices and global changes impact it. There is increasing global importance of wetlands, as there has been an extensive wetland loss in the last half-century due to anthropogenic pressures for land-use changes and development (Henderson, Lewis, 2008). A detailed and accurate wetland inventory would be a useful tool to protect them (Reimer, 2009).
Developing a mapping method to extract information about wetland areas from satellite imagery is essential for mapping large scale regions such as the Province of New Brunswick (Canada). Satellite remote sensing has several advantages over other methods for large area mapping, such as aerial photograph interpretation and ground surveys, as they provide multitemporal data with large footprints and moderate resolution. Most satellite images are readily available from multiple dates, allowing multi-seasonal analysis, which improves mapping accuracy and is more cost-effective than other mapping methods (Li, Chen, 2005). Satellite images provide a practical approach to mapping wetlands in New Brunswick, given the remoteness of some parts of the province. Optical images like those acquired by Landsat, SPOT or Sentinel-2 satellites were already tested for mapping wetlands (e.g., Amani et al. 2017). However, accuracy was quite low. Optical imagery has the additional inconvenience of having an availability, which is limited to clear sky conditions. By contrast, Synthetic Aperture Radar (SAR) images are of longer wavelengths (cm-scale vs. nm-scale), which penetrate clouds and, therefore, can be acquired whatever the weather conditions. SAR imagery has long been known as suitable for wetland mapping because radar waves can more easily penetrate the vegetation canopy for the detection of flooded areas and are sensitive to moisture conditions (Brisco et al., 2011;2013b;White et al., 2014;Mahdianpari et al., 2017). Both L-and C-band images were shown to be appropriate for wetland mapping. Lband can penetrate the forest canopies and detect standing water (Henderson, Lewis, 2008). C-band data have also been useful in detecting standing water in the case of short vegetation (Li, Chen, 2005;Henderson, Lewis, 2008), and in some forested wetlands with low-density canopies or leaf-off conditions (Townsend, 2002). The best approach will be to combine both optical and SAR imagery. There are several studies on wetland mapping that use this approach (e.g., Li, Chen 2005;Bourgeau-Chavez et al. 2009;2015;2016;Corcoran et al. 2012Corcoran et al. , 2013LaRocque et al. 2014, Amani et al. 2018Jahncke et al. 2018;Mahdianpari et al. 2020), but the accuracy achieved by these studies was rarely above 90%.
In this study, a combination of optical Landsat5-TM images with Radarsat-2 C-band and/or the Alos-PalSAR L-band SAR dualpolarized (HH, HV) images was tested for mapping wetlands in central New Brunswick, Canada. To take account of the local topography, a digital elevation model (DEM) was also used in the classification. The maps were produced using the Random Forests (RF) classifier applied to various image combinations. By contrast to the standard supervised maximum likelihood classifier (MLC), RF is a supervised non-parametric classifier that does not require normal data distribution (Breiman, 2001). Also, compared to the classical MLC, RF was showed to outperform in several land cover studies (Pal, 2005;Gislason et al., 2006;Waske, Braun, 2009;LaRocque et al., 2014). For assessing the map accuracy, the resulting maps were compared to GPS field data as well as to a combination of two wetland maps used by the Province of New Brunswick.

Study area
The Greater Fredericton area is in the heart of the Province of New Brunswick (Canada), lying between Lat 45° 40' and 46° N and Long 66° 15' and 66° 50' W ( Figure 1). This area (approximately 2,700 km 2 ) is relatively easy to access and contains all the wetland classes that can be found in Eastern Canada, such as peatlands (which include bogs and fens), marshes, aquatic beds, shrub wetlands, and forested wetlands. The drainage network is dominated by the St. John River, with some major tributaries. As indicated by the 1:50,000 DEM (Figure 1), the elevation in the study area ranges from 1 m above mean sea level (AMSL) in the downstream portion of the St-John River to 225 m AMSL in the highlands of the western part of this region. The topography is rolling in most of this area, except in the lowest portion of the St. John River valley and the Grand Lake basin (east part of the study area) where the surface is mostly flat. Sixteen land cover classes were considered in this study; namely, (1) Urban dense, (2) Urban sparse, (3) Cultivated, (4) Pasture,

Data
Three types of satellite imagery were used in the classification (Table 1) . Images used in this study were acquired during high-water levels (May to early June) and low-water levels in the wetlands (Mid-August to Mid-September) ( Table 1). All images were collected during ice-free and snow-free periods, and the soil was completely unfrozen.
Six Alos-1 PalSAR L-band dual-polarized images were acquired during three different dates (Table 1), with an incident angle at about 38°, from an ascending orbit looking east. For each date of acquisition. Three Radarsat-2 dual-polarized images were also acquired (Table 1) with a Standard beam mode (S6), and an ascending (A) or descending (D) orbit. The S6 beam mode corresponds to incident angles ranging from 41° to 46°. The images were acquired with an ascending east-looking orbit and with a descending west-looking orbit. The ascending orbit occurs during the evening hours and the descending orbit during the early morning. Both types of SAR imagery were acquired under different moisture conditions, as shown by the precipitation data recorded at the Fredericton Airport weather station ( Finally, the wetland map extracted from the classified image was compared to the wetland map derived from the combination of two maps currently in use by the Province of New Brunswick: (1) the wetland map from the Department of Environment, and (2) the wetland classes in the forest map of the Department of Natural Resources. Both maps were downloaded in a shapefile format from the GeoNB website (www.snb.ca/geonb1/e/DC/catalogue-E.asp). The map produced from the combination of these two maps will be called hereafter "NB wetlands map". The image-based maps and the NB wetland map were compared to GPS that were collected during summers between 2012 and 2018. 830 field sites that were spread throughout the study area were visited (Table 2). Less than half of these sites (340) were used to delineate training areas for the image classification. The remaining sites (490) were used as validation sites to assess the mapping accuracy of each classified image and the NB wetland map. These sites were selected according to various criteria, among them a relatively large extent (at least 10 pixels or 9000 m 2 ) of the related land cover class.  Among all the visited sites, 346 sites (147 training sites and 199 validation sites) were considered as a wetland site, according to Tiner (1999)'s criteria: the water table was close to (less than 10 cm) or at the surface, or we found indicator plants, soil hydromorphy, or other evidence of an area that is very often saturated with water. Most of the wetland sites were found by interpretation of aerial photographs, and some of them were previously mapped as a wetland on the NB wetlands map. On each field site, GPS location, elevation, class identification, and ground photographs were recorded.

Image processing
Most of the image processing was performed in PCI Geomatica® (PCI Geomatics, 2018). The DEM tiles were first imported and checked for holes, where holes were filled using interpolated values and then mosaicked together using the PCI OrthoEngine® module of PCI Geomatica®. The digital numbers of the Landsat-5 TM images were first converted into top of atmosphere (TOA) reflectance values using the equations of Chander et al. (2009). Such a conversion also removes some of the atmospheric interference.
The Radarsat-2 dual-polarized (HH/HV) images were filtered for removing speckle, using a Gaussian filter developed by Eric Grunsky (Geological Survey of Canada, Ottawa), as in LaRocque et al. (2012). Speckle can be considered as noise, and its intensity must be attenuated to enhance fine details on SAR images (Goodman, 1976). Each image was then orthorectified with the "Radarsat-2 Rational Function Model" in the PCI OrthoEngine® module of PCI Geomatica® using the DEM and ground control points (GCP's). Between 100 and 200 GCP's (depending on the image) were extracted from orthorectified Landsat-5 TM data (NAD83 -UTM, zone 19, row T) for geo-correction. On average, the orthorectification process achieved a root mean square error (RMSE) of less than 1 pixel in the x-direction and less than 1 pixel in the y-direction.
All the Alos-1 PalSAR images were imported, calibrated and orthorectified using the ASF MapReady® software (Alaska Satellite Facility Engineering Group, 2013). Each orthorectified image was then imported into PCI Geomatica® and filtered for speckle reduction, using a Gaussian filter, as for the Radarsat-2 images. The filtered and orthorectified Alos-1 PalSAR images were then mosaicked together using the "Automatic Mosaicking" processing in PCI Geomatica OrthoEngine® that requires the use of the DEM. Radiometric differences between adjacent images were minimized using a tonal-balancing histogram-matching method that compares the histogram of the image to be inserted in the mosaic with the mosaic histogram for each polarization.
Finally, all the datasets (Landsat-5 TM, Alos-1 PalSAR, Radarsat-2, and DEM) were re-projected to the New Brunswick Double Stereographic NAD83 (CanNBnad83) datum. All the input data were converted to the pixel spacing of 30 m of the Landsat5-TM imagery, allowing a direct superposition of all layers.

Image classification
Representative training areas of each of the sixteen land cover classes, based on 340 training sites (Table 2), were delineated over the orthorectified Landsat-5 TM. Training areas included at least ten sites per class for adequate class representation and were delineated within relatively large features, having exactly 10 pixels each in size. The training areas were randomly drawn throughout the study area, but the number of training areas by class reflected the relative frequency of the different land cover classes in the study area.
A non-parametric decision tree type classifier, Random Forests (RF), was applied to various combinations of DEM, Landsat-5 TM images, HH and HV polarization SAR images either from Alos-1-PalSAR L-band and/or Radarsat-2 C-band. RF does not assume a normal distribution of the input data (Breiman, 2001), The algorithm used in this study is an adaptation of the code of Horning (2010) by the first author. We used as settings a forest of 500 independent decision trees with the default value for the mtry variable, which is the number of variables randomly sampled as candidates at each node split. The default value includes all input features, i.e., all pixels are randomly sampled as candidates at each node split. RF randomly selects two third of the training data (that are referred as "In Bag" data) to develop one decision tree. This tree is then validated using the remaining third of the training data referred to as "Out of Bag data". The process is repeated for the 500 individual decision trees and produces 500 independent classifications. These independent classifications are then combined into the final classification (Waske, Braun 2009). When there are relatively limited training data for some classes, RF allows bootstrap aggregating of the "In Bag" data to increase the number of training pixels. RF is not sensitive to noise or over-classifying. RF gives an estimate of the importance of each input image for the classification (Pal, 2005;Gislason et al., 2006;Waske, Braun, 2009) under the form of a "Mean Decrease Accuracy" plot. The higher the image is on the "Mean Decrease Accuracy" plot Y-axis, the more useful the image was in the classification (Strobl et al., 2008;Louppe et al., 2013).

Accuracy assessment
Classification accuracy was assessed first by comparing training areas with the equivalent classified land use in the imagery. Such comparison was performed under the form of a "confusion" or "error matrix", where each cell expresses the number of pixels classified with the class defined by the training areas. The confusion matrix also allows computing the overall accuracy and the Kappa coefficient. The overall accuracy is the average of individual User's class accuracy, i.e., the probability that a pixel in the classified image belongs to the right class, weighted by the size of the class in the classified or reference image (Congalton, 1991). Finally, the Kappa coefficient of Cohen (1968) is defined as a weighted measure of agreement between the number of wellclassified pixels, while a value close to 0% corresponds to poor classification accuracy and a value closer to 100% indicates a good classification accuracy.
The image classification accuracy only assesses the classified image accuracy, which is different than the field accuracy. A more robust accuracy assessment is to compare the resulting classified image with an independent set of field validation sites that have a GPS location. If the image returns the same class as the one observed at the validation sites, then the classified image is considered as correct, and the pixel related to this validation site is associated with a value of 1; otherwise, its value will be zero. A percentage of correct identification was then computed as a function of the total number of validation sites. Such a comparison was not made using a classical confusion matrix because we only use the validation sites for the wetland classes, and some misidentifications are due to confusion with nonwetland classes.

Effect of wetland water level on the images
The effect of the different water levels (high-water and low-water levels) in the wetlands is visible in the true-colour composite of the Landsat-5 TM images (Figure 2), as well as in the false-colour composites of the Radarsat-2 ( Figure 3) and of the Alos-PalSAR images ( Figure 4).   Also, during the flood event, forested wetlands appear distinctively bright (light green) on both SAR images ( Figure 5), as the interaction between the water and the trees produce a double-bounce scattering resulting to a bright return. In the case of the Landsat5-TM, no usable image acquired during a flood event was available because of the presence of an intermediate to high cloud cover during this event. Figure 5. False-colour composite of the SAR images acquired during a flood event (a) Radarsat-2 dual-polarized C-HH and C-HV image, and (b) Alos-1 PalSAR dualpolarized L-HH and L-HV image. Each composite is made as follows: HH in red, HV in green, and HH in blue.

Image classification
The overall accuracy and the Kappa coefficient are compared for the various image combinations of optical and both SAR images (Table 3). This table shows that, whatever the combination, the image classification accuracies are very often better when the images acquired during a flooding event are added to the image combination. The highest classification overall accuracy is achieved by using either the Landsat-5 TM and DEM data (94.9%) or when both Radarsat-2 and/or Alos-1 PalSAR dualpolarized SAR images are added to the classification (94.3%).  In RF, the importance of each input data in the classification using all the dataset is ranked on the "Mean Decrease Accuracy" plot ( Figure 6). It shows that the DEM is the most important input data in the image classification. The middle infrared (TM5) and the near-infrared (TM4) bands of Landsat-5 TM come in second place, the optical images acquired during low-water-levels being more important than the ones taken during high-water levels. The Radarsat-2 C-band images are in the top ten of the most effective input datasets. The other optical Landsat-5 TM bands and the Alos-1 PalSAR L-band images have lower importance in the image classification. However, SAR images acquired during a flooding event are more important in the classification than the other SAR images.
The classified images produced with the data combination having the two best classification accuracies are presented in Figure 7 and Figure 8. While both classifications have similar overall accuracies, using only the Landsat-5 TM and DEM data produces an image where there is an apparent large extent of the forested wetlands across the study area (Figure 7).  By contrast, the image produced with the addition of both SAR image types shows a less extent for the forested wetland ( Figure   8). Therefore, the overall classification accuracy is not a good indicator of the accuracy, and it is the need to validate the classification using data acquired over GPS sites that were not used in the classification. Figure 8. Classified image produced with the full combination of Landsat-5 TM, Alos-1 PalSAR, and Radarsat-2, and DEM.

Validation accuracies
Analyzing the performance of each image combination based solely on the image classification accuracies is not enough, and it is necessary to compare the classified images with independent validation data sets. This validation accuracy should come from the comparison between the classified image and independent sites identified on the base of field observations.
Out of the 199 wetland sites, only 74.1% were correctly identified from the Landsat-5 TM and the DEM classified image (Table 5). Adding SAR images in the classification, the number of correctly identified wetland validation sites is higher, particularly when the image acquired during a flooding event is also considered, the difference being higher when the Alos-1 PalSAR images are used alone or in combination with the Radarsat-2 images. The number of correctly identified sites is the highest when all the dataset is used (195 sites or 97.9%). When only one type of SAR images is used, the number of correctly identified sites is higher with Alos-1 PalSAR (180 sites or 90.5%) than with the Radarsat-2 images (176 sites of 88.2%). In all cases, these percentages of correctly identified wetland sites are well above the 89 correctly identified sites wetland sites (44.7%) using the NB wetlands map coming from the two maps produced by the Province of New Brunswick. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-3-2020, 2020 XXIV ISPRS Congress (2020 edition) Table 5. Overall statistics about the correct identification of the 199 wetland validation sites on the classified images, related with each image combination, or on the NB wetlands map When either Radarsat-2 and/or Alos-1 PalSAR data are used in the image classification, most of the misidentifications are related to wetland sites not being classified in the correct wetland class, and very few are wetland sites being classified into a non-wetland class (Table 6).  Table 6. Distribution of the poorly identified wetland validation sites as a function of the source of error When all the dataset is used, only four sites are poorly identified. Three sites are not in the proper wetland class, and the last one is mapped as a non-wetland class. Their rather small width characterizes these four sites, close the spatial dimension of one pixel. For the NB wetlands map, about half of the 100 poorly identified sites are associated with wetland validation sites not being mapped as a wetland, the remaining half being sites that are not classified in the right wetland class.

DISCUSSION
In the present study, the RF classifier was applied to multitemporal SAR and optical images combined with DEM data. Our study also showed that combining optical, C-band and/ or L-band SAR, and DEM data strongly improved the wetland mapping accuracies, particularly when using validation sites. The benefit to use images from multi-sensors for wetland mapping was already shown by Li and Chen (2005), Bourgeau-Chavez et al. (2009;2015;2016), Corcoran et al. (2012Corcoran et al. ( , 2013, LaRocque et al. (2014), Amani et al. (2018), Jahncke et al. (2018) and Mahdianpari et al. (2020). SAR and optical images are complementary. SAR has a unique ability to detect surface texture and provides information on scattering mechanisms that are related to surface roughness (and thus to the presence or absence of vegetation as well as to the vegetation type) and soil moisture content. Optical images, such as Landsat-5 TM, allow acquiring information on the reflective properties that are related to the presence or absence of vegetation, the vegetation type, and the surface moisture content if the canopy is sparse enough.
The variable importance plot ( Figure 6) shows that the DEM is the most important input data in the classification. A similar result was obtained by Jahncke et al. (2018). It can be explained by the fact that wetlands are usually found in lowlands. The second and third most important input variables are the Landsat-5 TM shortwave-infrared (TM5) and near-infrared (TM4) images. Harris et al. (2006) and Meingast et al. (2014) already showed that these bands are suitable for mapping wetlands, particularly those with low vegetation. However, both SAR images are less important in the classification than the Landsat-5 TM-4 and TM-5 bands. Finally, the Alos-1 PalSAR bands are less important in the classification than the Radarsat-2 images. The highest validation accuracy is achieved for all the wetland classes when both C-band and L-band images are used (Table 5). If only one of the two types of SAR images is used in the dataset, the number of well-identified GPS wetland validation sites is higher when the Alos-1 PalSAR L-band images are included in the image combination. Both SAR-type images used in the study have HH and HV polarizations. Figure 6 indicates that the crosspolarization (HV) is more important than the HH polarization for the Radarsat-2 C-band images, but not for the Alos-1 PalSAR Lband images. Therefore, the volume scattering, as measured by the HV polarization, seems to be more critical for the wetland mapping with the C-band. The highest penetration of the L-band Alos-1 PalSAR image induces a double-bounce scattering from the wetland areas that can be captured by the HH polarization. Such double-bounce scattering makes wetlands very distinct on the images (Figures 4a and 5b). L-HH images were already found to be suitable for mapping forested wetlands (Hess et al., 1990;Whitcomb et al., 2009) and peatlands (Yamagata, Yasuoka, 1993).
Our study used not only a multi-sensor approach but also a multitemporal data approach, as the satellite imagery was acquired during three different periods of water level inside wetlands (high-water level, low-water level and flooding) as well as in two different seasons (spring and summer). Combining multitemporal and multi-sensor approaches have already been shown to produce highly accurate wetland maps (Bourgeau-Chavez et al., 2009;2015;Corcoran et al., 2013). Such an approach allows capturing seasonal differences in vegetation and water level conditions, helping for better discrimination between different wetland types. Figure 6 also shows that the Landsat TM images acquired during low-water levels are more important than the ones acquired during high-water levels. This is probably because the information input in the classification from the optical images is more related to the vegetation status than to the water level in the wetlands. Whatever the SAR sensor, the SAR images acquired during the flooding event are more important than those acquired at another time, indicating that the presence of water in the wetlands as observed on the SAR images is critical for the classification.
This study produced a classified image using multi-sensor and multi-temporal satellite images. By contrast to the NB wetlands map, all the classified images produced are in better agreement with the validation sites. The few misidentifications of the GPS wetland validation sites are mainly due to some wetland sites not being classified in the right wetland class, and there are very few confusions with the non-wetland classes. In the case of the NB wetlands map, there is a lot of misidentifications of the GPS wetland validation sites. These worst results are coming mostly from the interpretation of panchromatic air photos where wet areas are difficult to be detected. About half of the misidentifications are associated with sites mapped into a nonwetland class, the remaining half being sites not classified in the right wetland class.

CONCLUSION
This study showed that using a dataset including Radarsat-2 (C-HH, C-HV), Alos-PalSAR (L-HH, L-HV) and Landsat-5 TM images acquired during different water levels in the wetlands, and a DEM improves the mapping of wetland areas in the Greater Fredericton area. This good result is confirmed by the evaluation of the classified images produced with the validation sites. The few misidentifications of the GPS wetland validation sites are mainly due to wetlands not being classified in the right wetland class. There are very few confusions with non-wetland classes.
In this study, the selection of satellite images considers the hydrological conditions in the wetlands (high-water level, lowwater level, and flooding), and the leaf-on conditions. Corcoran et al. (2012) previously showed that leaf-on/leaf-off conditions could influence the mapping accuracy of wetlands. Further work is needed with images acquired on leaf-off conditions, but the good accuracies obtained in this study may indicate that there would be perhaps not the need to add such images.
This study was limited to SAR images acquired using large incident angles (38° for the Alos-1 PalSAR and between 41° and 46° for the Radarsat-2 S6 images). Further work should be done by testing steeper incident angles since these angles allow better penetration of the wetland vegetation to detect flood conditions because of the shorter path-length, as already shown in Sokol et al. (2004). Also, using multi-beam mode images has been shown to provide better classification results than single beam modes (LaRocque et al., 2012;. In this study, SAR images were also limited to two polarizations (HH and HV). Both Radarsat-2 and Alos-1 PalSAR polarimetric images are available that allow advanced polarimetric analysis, including target decomposition techniques, to better model the various scattering mechanisms.
The availability of the new Sentinel-1 C-band SAR satellite and the Sentinel-2 optical satellite should facilitate the mapping of land cover and wetlands, thanks to their better spatial resolution. Such data can be useful to produce a map for the wetland inventory for a large country like Canada, as suggested by Mahdianpari et al. (2020). Nevertheless, such mapping needs to consider also the non-wetland classes as well as wetlands can be different from one region to another, as well as the comparison between resulting the classified image with validation sites.