Semi-Automated Classiﬁcation of Lake Ice Cover Using Dual Polarization RADARSAT-2 Imagery

: Lake ice is a signiﬁcant component of the cryosphere due to its large spatial coverage in high-latitude regions during the winter months. The Laurentian Great Lakes are the world’s largest supply of freshwater and their ice cover has a major impact on regional weather and climate, ship navigation, and public safety. Ice experts at the Canadian Ice Service (CIS) have been manually producing operational Great Lakes image analysis charts based on visual interpretation of the synthetic aperture radar (SAR) images. In that regard, we have investigated the performance of the semi-automated segmentation algorithm “glocal” Iterative Region Growing with Semantics (IRGS) for lake ice classiﬁcation using dual polarized RADARSAT-2 imagery acquired over Lake Erie. Analysis of various case studies indicated that the “glocal” IRGS algorithm could provide a reliable ice-water classiﬁcation using dual polarized images with a high overall accuracy of 90.4%. However, lake ice types that are based on stage of development were not effectively identiﬁed due to the ambiguous relation between backscatter and ice types. The slight improvement of using dual-pol as opposed to single-pol images for ice-water discrimination was also demonstrated.


Introduction
Lake ice is an important component of the cryosphere for several weeks to several months of the year in high-latitude regions [1]. Ice cover on lakes is known to have an effect on both regional climate and weather events (e.g., thermal moderation and lake-induced snowfall) and to be a robust indicator of climate change and variability due to its sensitivity to changing atmospheric conditions, air temperature, in particular [1][2][3]. The Laurentian Great Lakes have the world's largest freshwater surface covering an area of 245,000 km 2 . Ice cover of the Great Lakes has a strong impact on the regional climate, navigation, economic activities, and public safety [4]. Knowledge of ice conditions and variability on the Great Lakes is especially important for the shipping industry and marine resource management.
Monitoring detailed ice conditions on large lakes requires the use of satellite-borne Synthetic Aperture Radar (SAR) data that provide all-weather sensing capabilities, high resolution, and large spatial coverage. For most ice centers, SAR is considered to be the most satellite data source suitable for operational ice mapping to support navigation safety [5]. C-band satellite SAR data from ENVISAT,

Study Area
The selected lake for this study, Lake Erie, is the fourth largest of the North American Great Lakes by surface area (Figure 1). Lake Erie lies between 78° and 94°W with a surface area of 25,655 km 2 and an average depth of 19 m, making it the shallowest and smallest by volume of the Great Lakes. Lake depth determines the amount of heat stored in the lake body and thus is an important factor for ice freeze-up [31]. Deeper lakes have greater heat storage capacity so that they take longer to cool down and eventually freeze [32]. Therefore during winter, the Great Lakes do not usually freeze over completely, except the shallowest Lake Erie [33].  Over the past five years (2013-2017), the mean winter air temperature (average daily air temperature of December, January, and February) and mean annual temperature recorded at the New Glasgow weather station, near the northern shore of Lake Eire (Figure 1), was −0.64 • C and 8.83 • C, respectively (Table 1). The period of analysis is the cold winter of 2013-2014, which experienced an average winter temperature of −6.36 • C. During that winter, the constant low temperature caused 92.5% of the Great Lakes to become frozen by early March, which is the second largest ice coverage since 1973 [34]. The daily minimum and maximum air temperatures in the area of Lake Erie during this winter are shown in Figure 2. The hourly air temperature at the time of each RADARSAT-2 image acquisition used in the study is also shown. Ascending and descending scenes are acquired at 6 p.m. and 6 a.m. local time, respectively. Over the past five years (2013-2017), the mean winter air temperature (average daily air temperature of December, January, and February) and mean annual temperature recorded at the New Glasgow weather station, near the northern shore of Lake Eire (Figure 1), was −0.64 °C and 8.83 °C, respectively (Table 1). The period of analysis is the cold winter of 2013-2014, which experienced an average winter temperature of −6.36 °C. During that winter, the constant low temperature caused 92.5% of the Great Lakes to become frozen by early March, which is the second largest ice coverage since 1973 [34]. The daily minimum and maximum air temperatures in the area of Lake Erie during this winter are shown in Figure 2. The hourly air temperature at the time of each RADARSAT-2 image acquisition used in the study is also shown. Ascending and descending scenes are acquired at 6 p.m. and 6 a.m. local time, respectively.

Synthetic Aperture Radar
In this study, 26 dual-polarized (HH+HV) RADARSAT-2 ScanSAR images acquired over Lake Erie from January to April 2014 were used for lake ice classification ( Figure 3). The SAR acquisition dates for all the scenes are listed in Table 2. They have a nominal pixel spacing of 50 m × 50 m and a scene size of 500 km × 500 km in both range and azimuth directions. With a wide spatial coverage, ScanSAR Wide A (SCWA) mode is the main mode that is used by the CIS for operational lake ice monitoring. The incidence angle of the images ranges from 20° to 49°. The temporal resolution is around 1-3 days, with varying incidence angle ranges.

Synthetic Aperture Radar
In this study, 26 dual-polarized (HH+HV) RADARSAT-2 ScanSAR images acquired over Lake Erie from January to April 2014 were used for lake ice classification ( Figure 3). The SAR acquisition dates for all the scenes are listed in Table 2. They have a nominal pixel spacing of 50 m × 50 m and a scene size of 500 km × 500 km in both range and azimuth directions. With a wide spatial coverage, ScanSAR Wide A (SCWA) mode is the main mode that is used by the CIS for operational lake ice monitoring. The incidence angle of the images ranges from 20 • to 49 • . The temporal resolution is around 1-3 days, with varying incidence angle ranges.  Table 2. List of synthetic aperture radar (SAR) acquisition dates, time, ascending/descending mode of the RADARSAT-2 scenes used in the study and the incidence angle range for Lake Erie in each scene. The average wind speed across the whole lake within half an hour of the SAR acquisition time that are acquired from the National Oceanic and Atmospheric Administration (NOAA) Great Lakes Coastal Forecasting System (GLCFS) are also provided for the scenes with open water.

SAR Acquisition
Date (M/D/Y) Shaded dates are scenes when Lake Erie is fully ice covered.  Shaded dates are scenes when Lake Erie is fully ice covered.

Ascending (A)/Descending (D) Mode
An example of a RADARSAT-2 scene acquired on 19 January 2014 during the ice freeze-up period is shown in Figure 4. The scene covers Lake Erie (bottom), part of Georgian Bay (top), and Lake Ontario (right). Various ice types on Lake Erie are well captured at HH polarization. As shown, HV polarization is less sensitive to incidence angle and wind-induced open water surface roughness in contrast to HH polarization. The original images were down sampled from approximately 10,000 × 10,000 pixels to 5000 × 5000 pixels by performing a 2 × 2 block averaging (100 m pixel spacing). This reduces the image size by a factor of 4, which increases the computational efficiency. This is the same procedure followed at the CIS for ice charting from SAR imagery.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 27 An example of a RADARSAT-2 scene acquired on 19 January 2014 during the ice freeze-up period is shown in Figure 4. The scene covers Lake Erie (bottom), part of Georgian Bay (top), and Lake Ontario (right). Various ice types on Lake Erie are well captured at HH polarization. As shown, HV polarization is less sensitive to incidence angle and wind-induced open water surface roughness in contrast to HH polarization. The original images were down sampled from approximately 10,000 × 10,000 pixels to 5000 × 5000 pixels by performing a 2 × 2 block averaging (100 m pixel spacing). This reduces the image size by a factor of 4, which increases the computational efficiency. This is the same procedure followed at the CIS for ice charting from SAR imagery.

CIS Image Analysis Charts
Since in situ observations on large lakes like the Great Lakes is impractical and expensive, we used the Great Lakes image analysis charts that are provided by the CIS who has decades of experience interpreting SAR imagery as the reference data. They are produced through visual interpretation of the ice conditions from SAR imagery in near real time (within 4 h of acquisition). Visual interpretation of SAR imagery is conducted by an experienced ice analyst through a digital image display and vector drawing tool. In addition to SAR imagery, historic ice patterns and meteorological conditions are also analyzed to support the ice charting. Ice analysts manually draw polygons on the SAR image, and for each polygon, the ice types, and their estimated concentrations are reported in an oval shape "egg code" ( Figure 5). The egg code is the World Meteorological Organization (WMO) standard. The standard thickness-based lake-ice stage of development coding is given in Table 3. More information about the CIS image analysis charts can be found in [35].
Image analysis charts are color-coded using the international WMO code for total concentration (CT) or stage of development (SD) of ice in the area. The total ice concentration is often monitored for weather forecasting while stage of development is more useful for operational purposes, such as ship navigation. Color code CT represents the estimated total ice concentration (in tenths) in each manually drawn polygon. Fast ice that is fastened to the shore is 100 percent concentration by definition. Color code SD is the ice type that has highest partial concentration in that polygon. The color code for CT is used when concentration is more variable than stage of development, while the color code for SD is intended to be used when stage of development is more variable than concentration [36]. An example of the Eastern Great Lakes image analysis chart color coded as CT and SD drawn from the RADARSAT-2 19 January 2014 scene ( Figure 4) is shown in Figure 6.

CIS Image Analysis Charts
Since in situ observations on large lakes like the Great Lakes is impractical and expensive, we used the Great Lakes image analysis charts that are provided by the CIS who has decades of experience interpreting SAR imagery as the reference data. They are produced through visual interpretation of the ice conditions from SAR imagery in near real time (within 4 h of acquisition). Visual interpretation of SAR imagery is conducted by an experienced ice analyst through a digital image display and vector drawing tool. In addition to SAR imagery, historic ice patterns and meteorological conditions are also analyzed to support the ice charting. Ice analysts manually draw polygons on the SAR image, and for each polygon, the ice types, and their estimated concentrations are reported in an oval shape "egg code" ( Figure 5). The egg code is the World Meteorological Organization (WMO) standard. The standard thickness-based lake-ice stage of development coding is given in Table 3. More information about the CIS image analysis charts can be found in [35].
Image analysis charts are color-coded using the international WMO code for total concentration (CT) or stage of development (SD) of ice in the area. The total ice concentration is often monitored for weather forecasting while stage of development is more useful for operational purposes, such as ship navigation. Color code CT represents the estimated total ice concentration (in tenths) in each manually drawn polygon. Fast ice that is fastened to the shore is 100 percent concentration by definition. Color code SD is the ice type that has highest partial concentration in that polygon. The color code for CT is used when concentration is more variable than stage of development, while the color code for SD is intended to be used when stage of development is more variable than concentration [36]. An example of the Eastern Great Lakes image analysis chart color coded as CT and SD drawn from the RADARSAT-2 19 January 2014 scene ( Figure 4) is shown in Figure 6.

"Glocal" Iterative Region Growing with Semantics Classification
To map lake ice types and open water, image segmentation was performed on the 26 RADARSAT-2 scenes. The hierarchical region-based classification called "glocal" IRGS method [11] was used in the study. This is an unsupervised classification method that can identify homogeneous regions with arbitrary class labels. It has been integrated into the MAp-Guided Ice Classification (MAGIC) System that was developed by the Vision and Image Processing Research Group at University of Waterloo [37].
The "glocal" method used in the study was built upon the IRGS algorithm developed by [28]. It is an unsupervised segmentation algorithm proposed for SAR sea-ice mapping. IRGS is an edge-based method which uses edge penalty functions and a region growing technique [28]. Basically, IRGS first breaks the scene into small homogeneous regions using a watershed segmentation algorithm and assigns each region a class label using an MRF model [27]. Then, adjacent regions with the same assigned class are successively merged until the system energy is minimized. The classification and merging process are combined in a iterative manner until the merging cannot be performed further or the maximum number of iterations is reached [38]. This algorithm incorporates the edge strength, which improves the adaptivity of the spatial context model to non-stationarity situations [28,38]. Since IRGS is performed on regions rather than pixels, long range (tens or hundreds of pixels) of the spatial interactions between pixels are efficiently accounted [38].
For large lakes, such as the Great Lakes, large SAR scenes on the order of 500 km wide are necessary for mapping lake ice. However, such large scenes introduce statistical nonstationarities across the image as a result of the incidence angle effects. To minimize the effects of incidence angle, the hierarchical "glocal" IRGS classification approach that combines the "high-detail local" and "large-scale global" information was introduced [23]. The local step divides the image into a number of autopolygons and IRGS classification is performed separately for each autopolygon. This step over-segments each autopolygon into a number of homogeneous regions. The class statistics for each autopolygon can be considered stationary [23]. The global step then glues the regions from the local autopolygon classification and forms a classification across the whole scene. This hierarchical method identifies homogeneous regions in the full scene with arbitrary class labels.
The "glocal" IRGS classification is performed in the MAGIC system and the flowchart is shown in Figure 7. The inputs to the system are HH and HV polarizations of a RADARSAT-2 image and its corresponding landmask that masks out land and other lakes. First, a number of autopolygons (10 × 10 grid) are created on the HV image while using the watershed segmentation [39]. The HV polarization is chosen because it is less sensitive to incidence angle and wind-induced surface roughness. An example of autopolygon segmentation for the 19 January 2014 scene (see original imagery in Figure 4) is presented in Figure 8a. Each autopolygon is then segmented to five arbitrary classes using IRGS classification that is based on both HH and HV bands. This step over segments the image to ensure each class is almost homogeneous. Figure 8b shows this classification for each autopolygon. Following this, the global step glues the regions of the local autopolygon classification using IRGS and then completes the full scene classification with arbitrary class labels. Here, different number of arbitrary classes are segmented for the ice-water classification and lake ice type classification.
(1) Ice-water classification: The global step glues the regions into eight classes across the entire scene. Each of the eight classes is considered to be homogeneous and contains only water or only ice. An example of the full-scene "glocal" classification is presented in Figure 9a.
(2) Lake ice type classification: The global step glues the regions into four distinct classes across the whole scene. Each of the four classes is considered as one of the WMO standard thickness-based ice types.
The last step is to manually label the classes as ice types or water based on visual interpretation of SAR imagery. A labelled ice-water classification of the example scene is shown in Figure 9b.
Experimentation was accomplished to ensure the parameters of "glocal" classification (10 × 10 grid, five classes for local, four or eight classes for global) used in the study performed well for the Remote Sens. 2018, 10, 1727 9 of 27 26 tested scenes. To over segment the most homogenous regions, best ice-water classification results are obtained when eight classes are segmented. The number of classes that are higher than eight does not provide more homogeneous segmentation results. For ice type classification, only four classes are segmented because it is difficult to visually label different ice types. The performance of the algorithm was not sensitive to the minor adjustment of grid size and number of classes for local step.

Return:
Labelled lake ice classification 26 tested scenes. To over segment the most homogenous regions, best ice-water classification results are obtained when eight classes are segmented. The number of classes that are higher than eight does not provide more homogeneous segmentation results. For ice type classification, only four classes are segmented because it is difficult to visually label different ice types. The performance of the algorithm was not sensitive to the minor adjustment of grid size and number of classes for local step. Labelled lake ice classification

Dual-Pol vs. Single-Pol
The advantage of utilizing dual polarization as compared to single polarization for ice observation is evident particularly for ice-water classification. One of the main challenges in discriminating open water from ice is the wind-affected water in co-polarized images. Windroughened open water has similar backscatter to ice in a single co-polarized image particularly at near range, whereas cross-polarization is less sensitive to wind effects and it can be useful for icewater discrimination [40]. However, there is limited ice-water contrast in cross-polarized backscatter because of the relatively low signal-to-noise ratio [26]. Therefore, a combined use of co-and crosspolarized data should improve the potential for ice-water classification. To test the advantages of using the dual-polarization (HH+HV) opposed to single HH polarization, "glocal" IRGS ice-water classification was also performed on only HH-polarized images.

Accuracy Assessment
Classification results were assessed against the CIS image analysis charts. The charts are manually drawn directly from the same RADARSAT-2 scenes that were used in the study on the same date by trained ice analysts at the CIS. They are currently the most reliable reference data available for the ice cover over the Great Lakes. However, the quantified image analysis charts contain uncertainty from the subjectivity of visual interpretation by ice analysts and its coarse representation [8,20]. This study assumes that the CIS image analysis charts represent the "true" lake ice conditions.
For ice-water classification, the absolute difference between the charts and classification results was used as the classification error. Image analysis charts that were color coded in CT were reclassified into lake ice and open water using a total ice concentration threshold of 10%. This is the lowest ice concentration identified by ice analysts. Polygons with an ice concentration smaller than 10% were reclassified as open water, otherwise they were considered as ice. For each scene, open water and lake ice correspondence as well as the overall accuracy were obtained. Since image analysis charts are region-based products that have greater level of generalization than our pixel-based classification, per-pixel accuracy assessment was also conducted using 400 randomly selected sample points per scene labelled based on visual interpretation of SAR imagery. To provide an insight on how our results compare directly to the total ice concentration in the egg codes, we also provided an example for a polygon-based comparison for one of the scenes (1 April 2014).

Dual-Pol vs. Single-Pol
The advantage of utilizing dual polarization as compared to single polarization for ice observation is evident particularly for ice-water classification. One of the main challenges in discriminating open water from ice is the wind-affected water in co-polarized images. Wind-roughened open water has similar backscatter to ice in a single co-polarized image particularly at near range, whereas cross-polarization is less sensitive to wind effects and it can be useful for ice-water discrimination [40]. However, there is limited ice-water contrast in cross-polarized backscatter because of the relatively low signal-to-noise ratio [26]. Therefore, a combined use of co-and cross-polarized data should improve the potential for ice-water classification. To test the advantages of using the dual-polarization (HH+HV) opposed to single HH polarization, "glocal" IRGS ice-water classification was also performed on only HH-polarized images.

Accuracy Assessment
Classification results were assessed against the CIS image analysis charts. The charts are manually drawn directly from the same RADARSAT-2 scenes that were used in the study on the same date by trained ice analysts at the CIS. They are currently the most reliable reference data available for the ice cover over the Great Lakes. However, the quantified image analysis charts contain uncertainty from the subjectivity of visual interpretation by ice analysts and its coarse representation [8,20]. This study assumes that the CIS image analysis charts represent the "true" lake ice conditions.
For ice-water classification, the absolute difference between the charts and classification results was used as the classification error. Image analysis charts that were color coded in CT were reclassified into lake ice and open water using a total ice concentration threshold of 10%. This is the lowest ice concentration identified by ice analysts. Polygons with an ice concentration smaller than 10% were reclassified as open water, otherwise they were considered as ice. For each scene, open water and lake ice correspondence as well as the overall accuracy were obtained. Since image analysis charts are region-based products that have greater level of generalization than our pixel-based classification, per-pixel accuracy assessment was also conducted using 400 randomly selected sample points per scene labelled based on visual interpretation of SAR imagery. To provide an insight on how our results compare directly to the total ice concentration in the egg codes, we also provided an example for a polygon-based comparison for one of the scenes (1 April 2014).
For ice type classification, the image analysis charts that were color coded in SD were used for visual comparison. The ice type classification results were not quantitatively assessed because image analysis charts are coarse and do not contain homogeneous ice type polygons. The lake ice types in the image analysis charts are based on stage of development, and each category contains a wide range of ice thicknesses (Table 3). Although the CIS provides the stage of development color code, most of the polygons are heterogeneous and contain multiple ice types without identifying the exact locations of the different ice types. In most cases, multiple ice types are mixed together in the ice analysis polygons and small features are not included. Therefore, it is visually difficult to recognize different ice types from the SAR imagery without ancillary data.

Overall Results
The semi-automated "glocal" IRGS classification was tested on 26 ScanSAR Wide RADARSAT-2 scenes and the results were compared against the CIS image analysis charts. For ice-water classification, only 12 scenes were included because Lake Erie was fully ice covered for the remaining 14 scenes (Table 2). Figure 10 shows an example of the validation process for both ice-water and ice type classification for the 19 January scene. The HH and HV images are shown in Figure 9a,b. The ice-water classification result (Figure 10e) was compared with the binary image analysis chart (Figure 10d) that was reclassified from the original image analysis chart color coded in CT (Figure 10c). To perform a quantitative accuracy assessment, error maps based on pixel-by-pixel differences between the reclassified image analysis charts and "glocal" IRGS results were generated (Figure 10f). The lake ice type classification result (Figure 10h) was visually compared with the image analysis chart color coded in SD (Figure 10g). For ice type classification, the image analysis charts that were color coded in SD were used for visual comparison. The ice type classification results were not quantitatively assessed because image analysis charts are coarse and do not contain homogeneous ice type polygons. The lake ice types in the image analysis charts are based on stage of development, and each category contains a wide range of ice thicknesses (Table 3). Although the CIS provides the stage of development color code, most of the polygons are heterogeneous and contain multiple ice types without identifying the exact locations of the different ice types. In most cases, multiple ice types are mixed together in the ice analysis polygons and small features are not included. Therefore, it is visually difficult to recognize different ice types from the SAR imagery without ancillary data.

Overall Results
The semi-automated "glocal" IRGS classification was tested on 26 ScanSAR Wide RADARSAT-2 scenes and the results were compared against the CIS image analysis charts. For ice-water classification, only 12 scenes were included because Lake Erie was fully ice covered for the remaining 14 scenes (Table 2). Figure 10 shows an example of the validation process for both ice-water and ice type classification for the 19 January scene. The HH and HV images are shown in Figure 9a,b. The icewater classification result (Figure 10e) was compared with the binary image analysis chart ( Figure  10d) that was reclassified from the original image analysis chart color coded in CT (Figure 10c). To perform a quantitative accuracy assessment, error maps based on pixel-by-pixel differences between the reclassified image analysis charts and "glocal" IRGS results were generated (Figure 10f). The lake ice type classification result (Figure 10h) was visually compared with the image analysis chart color coded in SD (Figure 10g). Overall, there is generally a good agreement between the image analysis charts and our icewater classification results. The average overall accuracy of the ice-water classification is 85.1% and the highest overall accuracy reaches 93.8% (Table 4). Most of the errors come from misclassifying ice as open water. Although the HH backscatter for open water is highly dependent on incidence angle as well as wind speed and directions [20], the "glocal" IRGS algorithm can provide a good ice-water discrimination (overall accuracy > 85%) in most conditions. The random pixel sampling shows an overall accuracy of 90.4%. Examples of different cases are analyzed in Section 5.2.1. Table 4. Accuracy assessment of ice-water classification for 12 scenes.  Overall, there is generally a good agreement between the image analysis charts and our ice-water classification results. The average overall accuracy of the ice-water classification is 85.1% and the highest overall accuracy reaches 93.8% (Table 4). Most of the errors come from misclassifying ice as open water. Although the HH backscatter for open water is highly dependent on incidence angle as well as wind speed and directions [20], the "glocal" IRGS algorithm can provide a good ice-water discrimination (overall accuracy > 85%) in most conditions. The random pixel sampling shows an overall accuracy of 90.4%. Examples of different cases are analyzed in Section 5.2.1. Table 4. Accuracy assessment of ice-water classification for 12 scenes.

Pixel-by-Pixel Difference Against Original SAR Images (400 Randomly Selected Pixels per Scene)
Overall Accuracy For comparison between classification using single-pol and dual-pol, the results showed that there is generally a slight improvement of using dual-pol rather than single HH polarization ( Figure 11). When compared with the image analysis charts (pixel-by-pixel difference), the average overall accuracy using single-pol for 12 scenes is 83.4%, which is around 3% lower than the accuracy for dual-pol. It is worth noting that for scenes that showed little difference (<5%) between using single-pol and For scenes with a good ice-water contrast in the HH polarization (1/12/2014, 1/14/2014, 1/19/2014, and 3/25/2014); however, the contribution of the HV polarization is often negligible. Overall, introducing HV polarization slightly improved the ice-water classification. Example of scenes for comparison will be discussed in Section 5.2.2.
For ice type classification, visual comparison shows that "glocal" IRGS algorithm can generally identify different ice types with single-pol and dual-pol imagery, but only when the surface conditions are different for different ice types. New lake ice is often confused with open water due to their similar low backscatter range. In most cases, thick lake ice and very thick lake ice appear to be similar in the SAR imagery, so that they were not discriminated from each other in most cases. Therefore, the algorithm has limited capabilities classifying different ice types. Specific scenes are analyzed in Section 5.2.3.
The total time for ice classification of Lake Erie in a 5000 × 5000 pixels scene is under 10 min, which is much faster than the manual ice charting at the CIS. The CIS ice analyst takes around 2 h to produce a Great Lakes image analysis chart from a large swath SAR image. The local autopolygon classification takes around 1 min, the global IRGS gluing takes 30 s, and the manual labelling takes about 2-5 min, depending on the complexity of the scene. For comparison between classification using single-pol and dual-pol, the results showed that there is generally a slight improvement of using dual-pol rather than single HH polarization ( Figure  11). When compared with the image analysis charts (pixel-by-pixel difference), the average overall accuracy using single-pol for 12 scenes is 83.4%, which is around 3% lower than the accuracy for dualpol. It is worth noting that for scenes that showed little difference (<5%) between using single-pol and dual-pol images, this difference can be ignored because some inconsistencies can be attributed to the coarse spatial resolution of the image analysis charts. The inclusion of HV polarization can improve the ice-water classification accuracy, especially at near range where wind-roughed open water tends to be as bright as lake ice (1/15/2014 in Figure 12a, 1/18/2014, 3/28/2014, and 4/4/2014). For scenes with a good ice-water contrast in the HH polarization (1/12/2014, 1/14/2014, 1/19/2014, and 3/25/2014); however, the contribution of the HV polarization is often negligible. Overall, introducing HV polarization slightly improved the ice-water classification. Example of scenes for comparison will be discussed in Section 5.2.2.
For ice type classification, visual comparison shows that "glocal" IRGS algorithm can generally identify different ice types with single-pol and dual-pol imagery, but only when the surface conditions are different for different ice types. New lake ice is often confused with open water due to their similar low backscatter range. In most cases, thick lake ice and very thick lake ice appear to be similar in the SAR imagery, so that they were not discriminated from each other in most cases. Therefore, the algorithm has limited capabilities classifying different ice types. Specific scenes are analyzed in Section 5.2.3.
The total time for ice classification of Lake Erie in a 5000 × 5000 pixels scene is under 10 min, which is much faster than the manual ice charting at the CIS. The CIS ice analyst takes around 2 h to produce a Great Lakes image analysis chart from a large swath SAR image. The local autopolygon classification takes around 1 min, the global IRGS gluing takes 30 s, and the manual labelling takes about 2-5 min, depending on the complexity of the scene.

Analysis of Specific Cases
The performance of lake ice classification can vary with incidence angle, wind speed, and direction, as well as ice conditions. Here, a sample of scenes exhibiting different ice and open water conditions are analyzed and discussed.

Ice-Water Classification
The scene on 15 January, presented in Figure 12, is a typical case where Lake Erie is located in the near range (incidence angle <35 • ) with the smallest incidence angle in the east of the image. This scene is acquired during the ice freeze-up and it mainly contains thin lake ice and medium lake ice. A Landsat-8 image acquired one day before is shown in Figure 12f. The seasonally-predominant southwesterly winds push the ice to the northeastern part of the lake. The HH backscatter of open water is highly dependent on the incidence angle with a large increase in backscatter towards smaller incidence angles, as shown in Figure 12a towards the right portion of the image. This is because at near range, open water with a high dielectric constant has a stronger surface scattering of the radar energy than that of ice [4]. With wind-induced surface roughness, open water can have even higher backscatter than the lake ice. Using this ice/water contrast, Leshkevich and Nghiem [4] successfully separated lake ice and open water using a simple threshold of HH backscatter at near range. This scene achieves an overall accuracy of 89.5%. Some errors appear at the ice-water boundaries where the contrast of the edges of open water and ice is less evident.
The scene on 21 February is a case during the spring melt when wet ice and snow or melt ponds on top of the lake ice are present ( Figure 13). This scene was acquired when the minimum air temperature was higher than the freezing point during that day (Figure 3). The decaying lake ice in this scene has relatively low backscatter because melting ice and wet snow on the ice surface absorbs most of the radar signal [41]. Since this is a near-range scene, there is generally a good contrast between the bright wind-roughed open water and dark lake ice. However, the bright brash ice with high angular topography in area A of the scene can barely be distinguished from open water at HH polarization, but HV can provide a good contrast. As a result, this scene achieves a high overall accuracy of 96.0%. Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 27 The next scene shown in Figure 14 was acquired 12 h after the 21 February scene in Figure 13. This is a scene acquired in the morning around 6 am local time (descending mode) of 22 February, whereas the last scene is acquired the night before at around 6 pm (ascending mode). When compared to the 21 February scene, the ice in this scene appears much brighter because the melted ice refroze as the temperature started to drop at 4 a.m. and reached the freezing point at 6 a.m. As shown in Figure 15, there was a significant melting episode before the 21 February scene was acquired as the air temperature for the previous 24 h was greater than 1.5 °C and reached 6.3 °C at 1 a.m. A Landsat-8 image acquired around 6 h after this SAR scene is shown in Figure 14e. Melting ice or melt ponds on the surface of the ice are likely present on the western part of the lake where the signal from the ice underneath is obscured. The low radar return of ice cover makes it difficult to discriminate from open water in the near-range. As shown in Figure 14d, ice in the western part was misclassified as open water and the overall accuracy is 81.0%. The next scene shown in Figure 14 was acquired 12 h after the 21 February scene in Figure 13. This is a scene acquired in the morning around 6 am local time (descending mode) of 22 February, whereas the last scene is acquired the night before at around 6 pm (ascending mode). When compared to the 21 February scene, the ice in this scene appears much brighter because the melted ice refroze as the temperature started to drop at 4 a.m. and reached the freezing point at 6 a.m. As shown in Figure 15, there was a significant melting episode before the 21 February scene was acquired as the air temperature for the previous 24 h was greater than 1.5 • C and reached 6.3 • C at 1 a.m. A Landsat-8 image acquired around 6 h after this SAR scene is shown in Figure 14e. Melting ice or melt ponds on the surface of the ice are likely present on the western part of the lake where the signal from the ice underneath is obscured. The low radar return of ice cover makes it difficult to discriminate from open water in the near-range. As shown in Figure 14d, ice in the western part was misclassified as open water and the overall accuracy is 81.0%.  The final classification result for the 1 April scene is shown in Figure 16. This is a simple scene for ice-water discrimination because there is a good backscatter contrast between open water and lake ice in both HH and HV polarizations. Although this scene is acquired during the break-up season where the melting ice weakens the backscatter, open water in the far range scatters away most of the radar signal and appears much darker than the decaying lake ice. This scene achieves a high correspondence of 93.8% with the image analysis chart. Differences between the classification results and image analysis charts are found at the ice-water boundary and some cracks between ice floes. This can be attributed to the coarse resolution of the image analysis charts. Errors also appear at the northeastern part of Lake Erie (circled area in Figure 16d), where some fast ice is misclassified as open water. A polygon-based comparison between the image analysis chart and the IRGS classification result (Figure 16d) is shown in Table 5. For each image analysis polygon in Figure 16c, the difference between the total ice concentration from IRGS result and egg code ice concentration is mostly within 20%. Since the ice concentration in egg codes is a coarse approximation that can be rounded off to 10% or higher, this comparison showed a good agreement between the image analysis chart and IRGS result. Table 5. Comparison between total concentration (CT) from egg code and IRGS classification result for each image analysis polygon in Figure 16c.

Polygon ID Egg Code CT (%) IRGS CT (%)
A  The final classification result for the 1 April scene is shown in Figure 16. This is a simple scene for ice-water discrimination because there is a good backscatter contrast between open water and lake ice in both HH and HV polarizations. Although this scene is acquired during the break-up season where the melting ice weakens the backscatter, open water in the far range scatters away most of the radar signal and appears much darker than the decaying lake ice. This scene achieves a high correspondence of 93.8% with the image analysis chart. Differences between the classification results and image analysis charts are found at the ice-water boundary and some cracks between ice floes. This can be attributed to the coarse resolution of the image analysis charts. Errors also appear at the northeastern part of Lake Erie (circled area in Figure 16d), where some fast ice is misclassified as open water. A polygon-based comparison between the image analysis chart and the IRGS classification result (Figure 16d) is shown in Table 5. For each image analysis polygon in Figure 16c, the difference between the total ice concentration from IRGS result and egg code ice concentration is mostly within 20%. Since the ice concentration in egg codes is a coarse approximation that can be rounded off to 10% or higher, this comparison showed a good agreement between the image analysis chart and IRGS result. Table 5. Comparison between total concentration (CT) from egg code and IRGS classification result for each image analysis polygon in Figure 16c.   Of the 12 scenes, the scene on 4 April is a difficult scene during spring melt ( Figure 17). The ice in this scene is mainly very thick lake ice (>70 cm) mixed with 20% to 30% of thick lake ice (30-70 cm). Since this scene is acquired during the melt season and the melting lake ice has a medium grey tone backscatter, which results in a very low contrast with open water. Some open water areas were misclassified as lake ice, but the overall accuracy achieved 87.0%. The algorithm also identified cracks in the northeastern part of the lake that were not delineated in the image analysis chart.

Polygon ID Egg Code CT (%) IRGS CT (%)
Remote Sens. 2018, 10, x FOR PEER REVIEW 20 of 27 Of the 12 scenes, the scene on 4 April is a difficult scene during spring melt ( Figure 17). The ice in this scene is mainly very thick lake ice (>70 cm) mixed with 20% to 30% of thick lake ice (30-70 cm). Since this scene is acquired during the melt season and the melting lake ice has a medium grey tone backscatter, which results in a very low contrast with open water. Some open water areas were misclassified as lake ice, but the overall accuracy achieved 87.0%. The algorithm also identified cracks in the northeastern part of the lake that were not delineated in the image analysis chart.

Dual-Pol vs. Single-Pol
The advantage of using dual-pol as opposed to single-pol images is more evident for scenes acquired at near range. As mentioned earlier in Section 5.2.1, the scene on 15 January is a near-range scene in which wind-roughed open water can appear just as bright as lake ice at HH polarization (area A in Figure 12a). Therefore, using only a HH image cannot discriminate some of the rougher ice types from open water (Figure 12e). However, HV polarization (Figure 12b) is less sensitive to incidence angle with open water appearing darker than the ice in area A. This provides a better backscatter contrast at the ice-water boundary. As shown in Figure 12d, the inclusion of both HH and HV polarizations in the algorithm produced a generally better ice-water classification result than using single polarization alone.
Most of the scenes analyzed showed that classification while using single-pol or dual-pol results in similar classification accuracy. For example, the scene of 1 April shown in Figure 16 is a case when using dual-pol is slightly better than using HH-pol only. Although the dual-pol result (Figure 16d

Dual-Pol vs. Single-Pol
The advantage of using dual-pol as opposed to single-pol images is more evident for scenes acquired at near range. As mentioned earlier in Section 5.2.1, the scene on 15 January is a near-range scene in which wind-roughed open water can appear just as bright as lake ice at HH polarization (area A in Figure 12a). Therefore, using only a HH image cannot discriminate some of the rougher ice types from open water (Figure 12e). However, HV polarization (Figure 12b) is less sensitive to incidence angle with open water appearing darker than the ice in area A. This provides a better backscatter contrast at the ice-water boundary. As shown in Figure 12d, the inclusion of both HH and HV polarizations in the algorithm produced a generally better ice-water classification result than using single polarization alone.
Most of the scenes analyzed showed that classification while using single-pol or dual-pol results in similar classification accuracy. For example, the scene of 1 April shown in Figure 16 is a case when using dual-pol is slightly better than using HH-pol only. Although the dual-pol result (Figure 16d) shows a better correspondence with the image analysis chart (Figure 16c), the result of using single HH-pol (Figure 16e) can identify more small cracks between ice floes, which were not delineated by the ice analyst.

Ice Type Classification
The lake ice type classification result for the 14 January scene is shown in Figure 18. There are three ice types including thin lake ice, medium lake ice, and new lake ice in the scene. The classified ice types (Figure 18d) have much different distribution as compared to the image analysis chart (Figure 18c). Some of the inconsistencies do not necessarily mean the results are erroneous because the ice expert at the CIS does not take into account small features in the image. Moreover, many of the image analysis polygons are not homogenous. They contain more than one ice type, as shown in the egg code (Figure 18c). New lake ice is misclassified as open water because its smooth surface results in similar low radar return as calm open water (area C in Figure 18c). Overall, most of the thin and medium lake ice is correctly discriminated from open water. shows a better correspondence with the image analysis chart (Figure 16c), the result of using single HH-pol (Figure 16e) can identify more small cracks between ice floes, which were not delineated by the ice analyst.

Ice Type Classification
The lake ice type classification result for the 14 January scene is shown in Figure 18. There are three ice types including thin lake ice, medium lake ice, and new lake ice in the scene. The classified ice types (Figure 18d) have much different distribution as compared to the image analysis chart (Figure 18c). Some of the inconsistencies do not necessarily mean the results are erroneous because the ice expert at the CIS does not take into account small features in the image. Moreover, many of the image analysis polygons are not homogenous. They contain more than one ice type, as shown in the egg code (Figure 18c). New lake ice is misclassified as open water because its smooth surface results in similar low radar return as calm open water (area C in Figure 18c). Overall, most of the thin and medium lake ice is correctly discriminated from open water.  The scene of 21 March is a complicated scene that contains various ice types ( Figure 19). The ice thickness in this scene ranges from less than 5 cm to more than 70 cm. Generally, the algorithm accurately identified new lake ice, medium lake ice, and very thick lake ice in the east part of the lake. However, thick lake ice cannot be differentiated from very thick lake ice due to their similar backscatter intensities. The scene of 21 March is a complicated scene that contains various ice types ( Figure 19). The ice thickness in this scene ranges from less than 5 cm to more than 70 cm. Generally, the algorithm accurately identified new lake ice, medium lake ice, and very thick lake ice in the east part of the lake. However, thick lake ice cannot be differentiated from very thick lake ice due to their similar backscatter intensities.

Classification Errors
Lack of ground truth data is one of the main challenges in evaluating ice classification algorithms [27]. The coarse representation of image analysis charts imposes limitations in the validation process. The resolution of image analysis charts is significantly coarser than our pixel-based classification results. In order for clients to easily read the charts, ice conditions in image analysis charts are expressed for large polygons that are at least tens of square kilometers in size. Furthermore, since the image analysis charts are produced near real-time on almost a daily basis, SAR image analyses are performed quickly by ice analysts and they do not take into account small features such as cracks between ice floes. Therefore, an absolute quantitative comparison between image analysis charts and classification results is not possible.
To compare our pixel-based ice-water classification results with the region-based image analysis charts, we reclassified the charts into ice and open water using an ice concentration threshold of 10%. Here, polygons with total ice concentration larger than 10% are regarded as lake ice, otherwise are reclassified as open water. This means that even regions with 20% of ice are considered fully ice covered, which explains why most errors come from misclassifying ice as water (Table 4). This subjective assumption contributes to some of the discrepancy and reduces the accuracy. Considering the limitations of the validation procedure, the classification accuracy is likely to be more reliable than the calculated pixel-by-pixel differences with the image analysis charts. Therefore, we have also conducted the per-pixel validation using random sample points from visual interpretation of SAR imagery and the overall accuracy reached 90%.
Although image analysis charts are produced by ice experts, they may contain some human-made errors. The total concentration, number of classes, and ice types are assigned subjectively from visual interpretation. A quantitative assessment of the CIS visual analysis that was done by [20] has indicated the existence of anomalies in image analysis charts. To ensure navigation safety, ice analysts are generally biased toward overestimation of ice concentration and assigning thicker ice types. Additionally, the interpretation of SAR images may vary between ice analysts depending on their level of expertise. Some analysts may include more details. A study by [8] showed that manually drawn ice charts by two different ice analysts at the Norwegian Ice Service for the same scene can be inconsistent to some degree. Even given the same polygon, significant differences of the ice concentration estimations were also found between different ice analysts at the Finnish Meteorological Institute [42].

Ice-Water Classification
There are two limitations of the "glocal" IRGS ice-water classification used in this study. One limitation is the difficulty in discriminating calm open water from new lake ice ( Figure 18). New lake ice is a thin layer of ice that forms at the beginning of the freeze-up period. It generally has a smooth surface/ice-water interface and high optical depth allowing for light to transmit to the subsurface water, which makes it visually transparent [43]. Therefore, new lake ice can be challenging to recognize from open water in either optical or SAR imagery. Ice experts at CIS often have to rely on weather conditions to determine the presence of new lake ice. Although a robust ice-water discrimination is desirable, the very thin new lake ice is less of an issue from an operational perspective [27].
The presence of water on ice due to rainfall or melt is another challenge for ice-water classification. Melting snow on ice can absorb or melt ponds reflect away most of the radar signal resulting in similar low backscatter to open water. Melt conditions are also a challenging situation for ice analysts interpreting SAR images.

Ice Type Classification
The correlation between backscatter and different ice types is ambiguous [18,19]. Each ice type (stage of development) can depict a large variability of backscatter that overlaps each other. Depending on the roughness at the surface and the ice-water interface, the same ice type can display very high backscatter with heavy ridging or low backscatter with a smooth surface. As shown in Figure 20c, the CIS ice analyst provided a ridging description of the same type of ice. As they get to thicker ice types (thick and very thick lake ice), ice analysts often have to rely on freezing degree days, as well as history and weather of the last few days to determine the stage of development. Therefore, it is difficult in practice to identify different ice types and estimate ice thickness based only on SAR backscatter. For further improvement, additional image characteristics, such as texture features, could be included [10,12,23,42].

Ice Type Classification
The correlation between backscatter and different ice types is ambiguous [18,19]. Each ice type (stage of development) can depict a large variability of backscatter that overlaps each other. Depending on the roughness at the surface and the ice-water interface, the same ice type can display very high backscatter with heavy ridging or low backscatter with a smooth surface. As shown in Figure 20c, the CIS ice analyst provided a ridging description of the same type of ice. As they get to thicker ice types (thick and very thick lake ice), ice analysts often have to rely on freezing degree days, as well as history and weather of the last few days to determine the stage of development. Therefore, it is difficult in practice to identify different ice types and estimate ice thickness based only on SAR backscatter. For further improvement, additional image characteristics, such as texture features, could be included [10,12,23,42].

Conclusions
The semi-automated "glocal" IRGS classification was tested on 26 dual polarized RADARSAT-2 imagery over Lake Erie. This is hierarchical region-based approach that divides the image into small polygons and minimizes the effect of incidence angle. This unsupervised classification method identifies homogeneous regions with arbitrary classes followed by manual labeling. Validation of the classification was done via pixel-by-pixel comparison with the CIS image analysis charts. The algorithm achieved an overall accuracy of 90.4% for ice-water classification. Analysis of various case studies indicates that the "glocal" IRGS algorithm can provide a reliable ice-water classification while using dual polarized images. The algorithm has difficulty in distinguishing calm open water from new lake ice and decaying ice. However, the misidentification of new lake ice is not a significant drawback for operational purposes. For lake ice type classification, most thin ice types were effectively identified but thick and very thick lake ice were often confused due to the ambiguous relation between backscatter and ice types. GLCM texture features could be included for further

Conclusions
The semi-automated "glocal" IRGS classification was tested on 26 dual polarized RADARSAT-2 imagery over Lake Erie. This is hierarchical region-based approach that divides the image into small polygons and minimizes the effect of incidence angle. This unsupervised classification method identifies homogeneous regions with arbitrary classes followed by manual labeling. Validation of the classification was done via pixel-by-pixel comparison with the CIS image analysis charts. The algorithm achieved an overall accuracy of 90.4% for ice-water classification. Analysis of various case studies indicates that the "glocal" IRGS algorithm can provide a reliable ice-water classification while using dual polarized images. The algorithm has difficulty in distinguishing calm open water from new lake ice and decaying ice. However, the misidentification of new lake ice is not a significant drawback for operational purposes. For lake ice type classification, most thin ice types were effectively identified but thick and very thick lake ice were often confused due to the ambiguous relation between backscatter and ice types. GLCM texture features could be included for further improvement. Overall, our "glocal" IRGS classification results are close to visual interpretation by ice analysts and would have expected to be closer if they could draw ice charts at a more detailed level. The testing of dual-and single-pol images demonstrated the slight improvement of ice-water discrimination utilizing dual polarized data as opposed to single polarization, particularly for near-range scenes.
Although the "glocal" IRGS algorithm needs manual labelling, the whole classification process takes less than 10 min for each scene, which is likely to be operationally useful. The algorithm can be fully automated by implementing automatic labelling while using trained models, as demonstrated by [23] for application to sea ice mapping.