Application of Dynamic Threshold in a Lake Ice Detection Algorithm

The traditional method involved in the classification of surface types such as water, ice, and snow rely on thresholds values of spectral properties that are fixed. However, the use of daily fixed thresholds leaves a substantial number of either unclassified and/or misclassified ice and water pixels. In this study, a new dynamic threshold technique is proposed to identify and map lake ice cover in the imagery of GOES-I to P series satellites. In addition, dynamic threshold can be used as an alternative solution to Bidirectional Reflectance Distribution Function (BRDF) models. The technique has been applied using GOES-13 imager data over Lake Michigan, one of five of the Great Lakes. Nine scenes based on an hourly acquisition of a single day are used to visually sample water and ice pixels. A threshold for the visible (0.62 μm) and the reflective component of the mid-infrared (3.9 μm) is determined for each scene by the intersection of the probability distributions of the water and ice samples. The thresholds are used as decision thresholds in a binary test algorithm applied on a per-pixel basis. Both fixed threshold (single scene) and dynamic thresholds (multiple scenes) have been compared. Dynamic threshold shows a significant gain in classified pixels over fixed threshold. A preliminary quantitative assessment is introduced to evaluate the algorithm’s performance using sensitivity and specificity testing. The classification results show a sensitivity of 98% when delineating thick ice and water and 87% when delineating thick/thin ice and water. Implementing a dynamic threshold, can be used in constructing ice maps in applications that benefit from high temporal resolution imagery.


Introduction
With the launch of GOES-16 in November 2016, a new era of remote sensing research over the North and South Atlantic basin and adjacent land masses has begun. The primary GOES-16 instrument is the Advanced Baseline Imager capable of providing high spatial and temporal resolution data in 16 bands. The remote sensing scientific and engineering community has developed various products based on this instrument including ice and snow maps [1]. This particular study is twofold. 1) Present the use of dynamic threshold in developing a sea and lake ice map [2]; 2) Promote further the use of the GOES-13 imager in snow and ice mapping [3,4] and the potential to develop these maps using historical imagery. It is the intention of the authors, that this work may contribute to the study of climate data in the Great Lakes region of which over 30 years of data from NOAA's GOES program are available.
Remote sensing based snow and ice maps for the North American region are generally limited to the short-wave infrared (SWIR) snow/ice detection bands (1.58-1.64 µm window) of certain multispectral sensors typically found on polar orbiting satellites. These include the Moderate Resolution Imaging Spectroradiometer (MODIS) on board NASA's Terra and Aqua satellites, the Advanced Very High Resolution Radiometer (AVHRR) on post NOAA-14 weather satellites, and the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi National Polar-Orbiting Partnership (Suomi NPP) spacecraft. Use of these optica sensors for snow and ice mapping are limited to clear sky conditions. Unlike polar-orbiting satellites, the higher temporal resolution of geostationary satellites offer greater opportunity in monitoring for cloud-free conditions and pixel classification via a daily composite map of cloud-free pixels. Prior to GOES-16, the GOES imagers were limited to 5 bands, none of which were in the SWIR needed for ice mapping. However, research in the use of the mid-infrared (MIR) 3.9 µm band for snow and ice maps have been investigated through the GOES imagers [3,4].
Ice mapping products derived from polar orbiting satellites are also limited in the use of constant threshold in their algorithms [2]. For example, an ice concentration algorithm for VIIRS uses fixed threshold for the visible (vis>0.08) and (NDSI>0.45) for ice identification [5]. A reflectance threshold, as referred to in this study, is a reflectance value shared by two dissimilar surfaces (differences in spectral properties). The ambiguity that exists when classifying these surfaces is predominantly a result of the particular viewing and illumination geometry at the time of data acquisition. For example, snow has a significantly higher albedo than liquid water (difference between albedo and reflectance will be discussed in a moment). However, when sunlight reflects off the surface of a water body at the same angle that a sensor is viewing the surface there is a significant "boost" in the reflectance of water, a phenomena known as sunglint; this can lead to classification errors while distinguishing ice from water. Ice mapping products with fixed or static thresholds that do not account for the variation in illumination geometry during the day are susceptible to these errors. Unlike polar orbiting satellites, which have low frequency of coverage offering at best a daily snapshot of the illumination geometry for a given area, geostationary satellites offer continuous coverage that allows observations across various illumination geometries. This information is useful in building a threshold that adapts to these changing lighting conditions.
The challenge for developing any surface type classification process, accounting for the ambiguity of dissimilar surfaces with overlapping spectral signatures, requires an understanding of multispectral surface reflection which depends upon both illumination and viewing direction, and is described by a particular surface's Bidirectional Reflectance Distribution Function (BRDF). Albedo is defined as the fraction of incident radiation that is reflected by a surface. The albedo of snow and ice is dependent on the solar zenith angle, . Reflectivity of snow and ice increases with increasing [6]. Albedo differs from reflectance in that albedo is the directional integration of reflectance over all possible sun-view geometries. Albedo depends on the spectral and angular distribution of incident radiation and thus is dependent on BRDF [7]. A Lambertian surface (completely diffuse) would have a BRDF plot that is cylindrical in shape with a nearly flat top; however, the reflectance of radiation from most natural surfaces is dependent on the sun-view geometry. Work has been done on laboratory measurements of BRDF [8], including field investigations [9] and the usefulness of these measurements in remote sensing applications. Application of laboratory measurements of BRDF on remote sensing data presents challenges including the broader scale of remote sensing instruments, complexity of physically-based BRDF models in transposing them onto remote sensing data, and the heterogeneity of land cover. These challenges have been investigated by measuring BRDF of various surface types through remote sensing instruments; including airborne [10] and geosynchronous [11] and evaluated against existing BRDF models. Today, natural surface BRDF databases exist that can be used to validate satellite sensors [10].
Dynamic threshold has been proposed in this study as an alternative to investing in a BRDF model for an ice mapping algorithm. In order to build a robust and reliable ice detection algorithm, it is necessary to properly account for the effects of the variable viewing and illumination geometry of observations. Conceptually, a dynamic threshold can be an alternative solution to accounting for these effects over the use of Bidirectional Reflectance Distribution Function (BRDF) models which mathematically correct for surface reflectance anisotropy during the classification process. As a continuation of previous work [4], which provided a preliminary evaluation of the MIR band to lake ice detection using a single threshold, the results of applying a dynamic threshold are presented in this paper. A dynamic threshold in this study is a time series of hourly reflectance thresholds based on ( ∩ ), where I is the probability of a pixel detected as ice and W is the probability of a pixel detected as water. These probabilities have been derived from hourly plots of normal distributions of ice and water pixels and applied to an ice detection algorithm that uses the mid-infrared band to delineate thick ice and water. Dynamic thresholds in the visible have been used in previous studies including ice mapping which has been demonstrated over the Caspian Sea [2] and in brightness temperature for cloud detection [12].

Method
The Great Lakes are a vital freshwater resource for the millions of people who live along the shores. In addition, the whitefish fishery is the most economically valuable commercial fishery in the upper Great Lakes. Climate change and variability may have implications on these resources. The Great Lakes in general are seeing a decrease in annual ice cover due to rising temperatures in the region that by the end of the century, are predicted to increase up to near 11°C in the summer and 0.5 to 9.1°C in the winter. Lake Michigan has been chosen as the study area for lake ice detection. Lake Michigan is one of five of the Great Lakes located in the northeastern Midwest United States along the U.S.-Canadian border The lake has a surface area of approximately 57,800 km 2 and a volume of 1,180 km 3 .

Satellite Channels
The GOES-13 satellite is in geosynchronous orbit, with a nadir located at approximately 0.05°N, 75°W. Data collected by the GOES-13 imager instrument is being used in this study. The GOES imager is a five-channel instrument (one visible, four infrared), as listed in Table 1. All channels except channel 3 are used in this study.

Approach and Algorithm Development
The Normalized Difference Snow Index (NDSI) is a snow index measurement that has traditionally been used in delineating snow and ice from other land surface features. NDSI is the normalized difference between the 0.6 µm visible band and the 1.6 µm SWIR band. An alternative to this, for sensors not operating in the SWIR window, is the ratio of the 0.6 µm visible band and the 3.9 µm MIR band. This method is chosen since GOES-13 does not operate in the SWIR.
Thick ice and snow have relatively high reflectivity in the visible. Ice reflectance substantially increases when ice thickness is above 5 cm, but noticeably lower reflectivity in the MIR. In this study, a new ratio of VIS to MIR called MISI (Mid-Infrared Sea and Lake Ice Index) is proposed to discern thick ice above 5 cm with a high MISI value and thinner or broken ice with a lower MISI value. GOES-13 channels 1 and 2, VIS (0.62 µm) and MIR (3.9 µm), are used in the development of MISI. Channel 4 (10.7 µm) is used to obtain the skin temperature. The reflectance values of channels 1 and 2, MISI, and skin temperature are the basis for ice classification in the algorithm. As this study is a continuation from previous work [4], the method for data acquisition, pre-processing, calibration and conversion, determination of skin temperature, and the derivation of the 3.9 µm reflection component has been implemented. Readers are encouraged to refer to that work for a detailed explanation of the pre-processing and processing procedures used in the algorithm.
As in the previous study [4], data was acquired for February 28, 2015. Acquisition times for that day were every half-hour from 1430 UTC to 2030 UTC. These are the times that GOES-13 is in the continental US (CONUS) extended scan mode with the exclusion of near 18:00 UTC which is during the time GOES-13 is operating in full disc mode. These times also correspond to daytime conditions over the eastern half of the North American continent. The mean solar zenith angle across all scenes was 56.8° with a range from 51.73° to 70.18°.
The distinction between various surface types in a satellite image becomes possible owing to their different spectral response. A threshold-based decision-tree image classification scheme is used to distinguish between water, gray ice, and thick ice. In this study, gray ice is referred to as thin or broken ice. The spectral properties of gray ice and thick ice are significantly different [4]. In this paper, R 1 is referred to the visible and R 2 as the reflective component of the mid-infrared. The R 1 , R 2 , MISI [4] thresholds were determined through a visual inspection of the visible reflectance for a sampling of pixels from southern Lake Michigan. The sampling consists of a mix of water and ice pixels. 50 water pixels and 50 ice pixels were identified from each acquired image. The scene images were rendered in gray scale. Cloud-free pixels relatively dark were visually classified as water, cloud-free moderately bright pixels were classified as ice. Example is shown in Figure 2. The identification of ice was visually validated against ice charts of the Great Lakes from the National Ice Center. Delineating the different ice concentrations as classified by NIC was not conducted. It should be noted that darker ice such as thin or broken ice will have reflectivity in the visible similar to water. The data is next best fitted to a normal distribution. Figure 3 are the distributions for R 1 , R 2 , and MISI at 1430, 1730, 2030 UTC. Blue is water, and red is ice. The number of bins for each distribution is equal to √ . As the sampling size was relatively small, the width of the bins is relatively wide. The R 1 distributions show a clear distinction in the mean value between water and ice. The ambiguity occurs where both distributions overlap. The overlapping indicates that a pixel has probability of being water or ice. For the scope of this paper, the point of intersection ( ∩ , where I is the probability of a pixel detected as ice and W is the probability of a pixel detected as water, is chosen as a threshold value for that particular time. All data values above this point are classified as ice and all data values below this point are classified as water. Since thick ice will have R 1 values significantly higher than water, it is reasonable to suggest that gray ice may occur at or near ∩ .   Table 2 lists the probability threshold values for R 1 (T r1 ), R 2 (T r2 ), MISI (T misi ) for each acquisition time. T r2 is calculated: T r1 /T misi . In addition, observation geometry is also provided: solar zenith angle (θ), solar azimuth angle (ϕ), solar-satellite relative azimuth angle (γ). The observation geometry has been calculated from a point located at 43.450°N, 87.222°W. GOES-13 azimuth and zenith angles at this location are approximately 162.22°, and 51.72° respectively. The fixed threshold values used in the original ice mapping algorithm [4] are replaced by the threshold values from Table  2. Figure 4 shows the revised detection algorithm. Sampling cloud pixels was not conducted as it was not a goal of this study; therefore, the threshold values from the previous study for cloud identification were used. Figure 5 compares the fixed to dynamic threshold values for R 1 . Fixed threshold has been selected at 1730 UTC which is close to local solar noon [4]. There is high variability in the dynamic threshold values near the end of the day, this may be an artifact in the relatively small volume of the training data set.  Figure 6 compares the lake ice map generated for 1430 UTC at a fixed threshold (left) generated from 1730 UT ( Figure 4) and dynamic threshold (right). A visual inspection reveals there are substantially less unclassified pixels. Figure 7 is the resulting gain of classified pixels across the ten acquisition times ( Table 2) for dynamic threshold over a fixed threshold. Gain is defined as the change in the number of pixels with a positive response against the same pixels tested using a static threshold (gain is zero). The most significant gains are near sunrise and sunset. Figure 8 is the final daily composite lake ice map. Application of a dynamic threshold provides less unclassified pixels, resulting in a more complete ice map.

Discussion
At this point, the application of dynamic threshold for this particular scene has resulted in substantially fewer unclassified pixels. The next step is to test the performance of the algorithm in its ability to delineate ice and water. The use of spatial analysis methods commonly used in Geographic Information System (GIS) has been used to quantify this performance. ESRI ArcGIS (10.2.1) is used in validating the classified pixels against Interactive Multisensor Snow and Ice Mapping System (IMS) and NIC ice maps. Vector and raster data are the two basic spatial data types used in GIS. Vector data are comprised of vertices and paths. Vertices can be used to construct lines and polygons (areas). The simplest vector data -vector points -are XY coordinates. Vector data are stored as XY coordinate pairs (latitude and longitude). Raster data are comprised of pixels (grid cells). Each pixel is associated with a value. Discrete rasters have distinct categories. These data may be stored as integer values to represent classes. For example, in land or surface cover, a value of 4 may represent snow, a value of 3 may represent sea ice. Raster data may be stored in the GeoTIFF file format which allows georeferencing information to be embedded within a TIFF file. Raster data can be converted into vector data which may then be stored as a shapefile. The shapefile format is developed and regulated by ESRI and allows geospatial data to be shared across multiple GIS platforms.
The U.S. National Ice Center provide ice and snow products in various formats including geotiff and shapefiles. In this study, the IMS data is obtained as a geotiff and the NIC ice concentration map is obtained as a shapefile. The MISI classification data stored as a matrix in Matlab is converted to vector data for use in GIS.
There are multiple techniques in performing quantitative analysis in GIS. The process of validation in this study involves testing a MISI pixel class that satisfies a boundary condition derived from the IMS and NIC maps that best represents the class. This decision boundary will be referred to as a catchment area. For the IMS, this is done by vectorizing the input raster. The process of vectoring a raster is shown in Figure 9. The polygons conform exactly to the raster's cell edges (non-simplified output). This method should work well as IMS classifications are relatively homogeneous. Table 3 presents the classification scheme of this model (MISI), IMS, and NIC. The MISI values are the classes and the IMS and NIC are values of the catchment areas.  To prepare for analysis in GIS, the classification map was down-sampled to a resolution of (4km) 2 /pixel. The resulting 3440 test pixels (16km 2 /pixel) total area of coverage is 55,040 km 2 (area of Lake Michigan is approximately 58000 km 2 ).
These 3440 test pixels that fall within the Lake Michigan boundary are spatially queried and checked for containment within each catchment area. Figure 10 shows the IMS ice catchment area in red and water catchment area in light blue (a). All MISI pixels are queried within the ice catchment (b); and within the water catchment (c). The percentage of each MISI pixel classification within each IMS catchment area are tabulated (Table 4). 2.56% remain unclassified. For the IMS product, the ice/water numerical fraction is 78.97%/20.35%; MISI product is 70.3%/12.61%.  Here the performance of the classification is discussed. Sensitivity and specificity are statistical measures of the performance of a binary classifier. Sensitivity measures the proportion of actual lake ice pixels which are correctly identified. Specificity measures the proportion water pixels that are correctly identified. Table 5 presents a breakdown of these measurements. From the performance evaluation table, the inclusion of gray ice in the classification results in lower performance (86.52%) than without it (94.01%). This is due to the ambiguity that exists between cold water and gray ice. Both gray ice and water have similar reflectivity and the only ice-water discrimination that occurs in the algorithm is from the skin temperature [4]; therefore, there are a significant number of misclassified water pixels as gray ice (Figure 9c).
NIC provides a more comprehensive ice concentrations delineation. The MISI model is compared. Figure 11 shows the NIC catchment areas for various ice concentrations (a). All MISI pixels are queried within fast sea ice catchment (b); within the medium 1 st year ice catchment (c); and within the ice free catchment.
The percentage of each MISI pixel classification within each NIC catchment area are tabulated ( Table 6). The table shows a possible delineation of thick ice and grey ice outside the medium 1 st year ice catchment area.

Conclusion
A dynamic threshold is an alternative to Bidirectional Reflectance Distribution Function (BRDF) correction which accounts for the biophysical, reflectance, and specular properties of a surface. In this study, a dynamic threshold based on lake ice detection method using an hourly threshold is developed. The threshold values are calculated on an hourly basis and applied to an ice detection algorithm through each iteration. The algorithm has been applied over Lake Michigan, one of five of the Great Lakes. Both fixed and dynamic thresholds have been compared. The resulting composite map using dynamic threshold yields significantly less unclassified pixels then the composite map using fixed threshold. A preliminary quantitative evaluation of the algorithm against IMS reveals good performance when delineating thick ice from water but lower performance when delineating gray ice from water. The use of a dynamic threshold may be used in constructing ice maps in applications that require higher temporal resolution. Future work proposed in this area includes additional yearly test scenes of the same region in identical solar-view geometries to validate the robustness of the dynamic threshold followed by the development of additional dynamic threshold for other solar-view geometries, in the same region.