A snow cover mapping algorithm based on a multitemporal dataset for GK-2A imagery

ABSTRACT Snow plays a crucial role in the climate. In particular, snow controls the surface temperature and stabilizes the energy budget due to its spectral characteristics. However, snow can be affected by climate change. If the snow cover extent (SCE) diminishes due to global warming, the rate of introducing radiance to the surface increases because non-snow areas can absorb more solar radiance. Therefore, it is highly important to accurately detect snow-covered areas. In this paper, we develop an algorithm for deriving a daily snow cover map of East Asia using Geostationary-Korean Multi-Purpose SATellite-2A (GEO-KOMPSAT-2A, GK-2A) imagery. After processes for cloud masking and misclassification removal are conducted, a threshold with consideration for the characteristics of a geostationary satellite is applied to the normalized difference snow index (NDSI), and a daily snow cover map is generated utilizing a stacking process. In the case of the cloud mask, the angle-time variation and displacement by clouds proposed in this study were used to enhance the accuracy of cloud detection. For the quantitative validation, the F1 score was 0.89 for Landsat-8 and 0.79 for the interactive multisensor snow and ice mapping system (IMS). In addition, when the snow cover extent calculated from the IMS was compared, the correlation reached 0.91 from December 2020 to January 2021.


Introduction
Snow is a mixed material containing air and ice crystals. Approximately 98% of the Earth's snow cover is located in the Northern Hemisphere (NSIDC 2020). Snow plays a crucial role in the climate, hydrological cycle, and ecosystems, as it affects the energy budget at Earth's surface (Awasthi and Varade 2021), thermal regimes (Callaghan et al. 2011), and weather (Wang et al. 2017). Specifically, snow has high shortwave albedo and facilitates Earth's energy balance by reflecting solar energy back into space, thus helping to cool the planet (Callaghan et al. 2011). The snow advance index (SAI), which is calculated from the average rate of change in the snow cover extent (SCE) derived from the daily snow cover map, can explain about 75% of the variance in Arctic oscillation (AO) in the winter season (Cohen and Jones 2011). AO refers to the variability in sea level pressure between the Arctic and mid-latitudes with opposite phases. In general, low pressure in the Arctic and high pressure in the mid-latitudes exist in winter. When this pattern is reinforced, it is called the AO's positive phase, and the opposite is defined as the negative phase. If the negative phase occurs in winter, the Arctic jet stream meanders, and cold air from the Arctic moves south, increasing the possibility of a cold wave in the Republic of Korea. Regarding the water cycle, snowmelt can provide water to approximately 17% of the world population (Hock et al. 2006); e.g. the snowgenerated runoff in the Arctic drainage basin contributes up to 75% of the total annual flow in some northern regions of Siberia and North America (Callaghan et al. 2011). Furthermore, snowmelt contributes to hydropower production, irrigation, and industry (Bormann et al. 2018).
Over the past 30 years, the monthly SCE of the Northern Hemisphere has decreased by 1.3% per decade (Barry et al. 2009). This reduction contributes to global warming by decreasing the surface albedo and affecting ecosystems (Barry et al. 2009). In addition, due to the variance in SCE, streamflow in highmountain Asia dropped by an average of 16% between 197916% between -199916% between and 199916% between -201916% between (Kraaijenbrink et al. 2021). In the northern area of Pakistan, the maximum and minimum temperatures increased by approximately 0.1°C per year from 1995 to 2010, and the minimum snow cover area shows a decreasing trend (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) (Qureshi et al. 2022). Therefore, monitoring the spatio-temporal variability in SCE is crucial.
Various snow cover map algorithms have been developed using satellite imagery; these algorithms are based on optical sensors, passive microwave sensors, or the integration of both.
In the case of optical sensors, the normalized difference snow index (NDSI) has been widely used to understand the presence of snow in pixels (de Wildt et al. 2007;Hall et al. 2001;Key et al. 2013;Lee et al. 2017;Wang et al. 2017;Zhu and Woodcock 2012). NDSI is based on the spectral characteristics of snow, reflecting its very high visible reflectance and very low reflectance in shortwave infrared (SWIR). Snow ranges are theoretically possible for pixels with NDSI values between 0 and 1 (Riggs, Hall, and Román 2016). Typically, an NDSI over 0.4 has been widely used as a threshold in cases where obvious snow-covered pixels need to be captured (Zhang et al. 2020). However, if only 0.4 is used as the NDSI threshold when the forest density is over 50% in the pixel, the classification accuracy decreases by approximately 25% because optical sensors detect viewable snowcovered areas (Klein, Hall, and Riggs 1998). The forest area is obscured by canopy or tall understory shrubs, resulting in less detectable snow cover (Rittger et al. 2020). In addition, if the view zenith angle (VZA) is high, the oblique perspective is introduced, and the detected area is stretched (Rittger et al. 2020). Therefore, Klein, Hall, and Riggs (1998) proposed a method consisting of the NDSI and the normalized difference vegetation index (NDVI) to accurately classify snow cover in forest areas. The method is used for pixels with NDSI from 0.1 to 0.4. In the case of the function of mask (Fmask) and the algorithm of de Wildt et al. (2007), 0.15 was used as the NDSI threshold to include pixels with snow coverage below 50% and snow-contaminated forest areas (de Wildt et al. 2007;Zhu and Woodcock 2012). Zhu and Woodcock (2012) found that the NDSI was below 0.15 in a clear area (without snow and clouds) in Landsat images. In addition, the normalized difference forest snow index (NDFSI), which is composed of a combination of the 0:86 μm and 1:6 μm bands, was used to supplement the NDSI in forest areas because the NDFSI had a more stable classification accuracy than the NDSI, with a variation of 0.05 in forest areas (Wang et al. 2019(Wang et al. , 2018. For most optical sensor-based algorithms, the cloud mask is first applied to find clearsky areas. However, some clouds and water regions misclassified as snow pixels may remain. Therefore, the snow cover mapping algorithms of the moderate resolution imaging spectroradiometer (MODIS) and visible infrared imaging radiometer suite (VIIRS) utilize the thresholds corresponding to the visible (0:65 μm), near-infrared (NIR, 0:86 μm), SWIR (1:6 μm), and brightness temperature (BT, 11:5 μm) bands to minimize misclassification in water, tropical forests, and mixed areas (Hall et al. 2001;Key et al. 2013). In some studies, variations in the BT, ratio of SWIR and NIR, and spectral-contextual snow index (SCSI) which is calculated by integrating the local standard deviation of the green band and NDSI were used to enhance the algorithm's accuracy (Qiu, Zhu, and Binbin 2019;Romanov, Gutman, and Csiszar 2000;Siljamo and Hyvärinen 2011).
Passive microwave sensors can efficiently generate snow cover maps because they are not sensitive to atmospheric components and can pass through rain, except for extremely heavy rain. Generally, snow cover maps based on passive microwave sensors are generated by using the variation in emissivity along with the frequency over snow and other land cover types (Grody 2008). In addition, since passive microwave can penetrate snow cover, the snow depth and the snow water equivalent can be acquired (Xinghua et al. 2019;Varade and Dikshit 2020). However, in vegetation areas covered by snow, microwavescattering signals are attenuated and can affect the accuracy of snow products (Xiao et al. 2018). Therefore, to decrease misclassification, the forest fraction in the pixel can be added into the equation; alternatively, various modified algorithms have been developed in accordance with various land cover types (Goïta, Walker, and Goodison 2003;Romanov and Tarpley 2007;Xiao et al. 2018). A disadvantage of microwave sensor-based snow cover products is their low spatial resolution due to the low amounts of emitted radiance.
Therefore, methods that integrate optical and microwave sensors have been developed to compensate for their limitations (Bergeron et al. 2014;Foster et al. 2011;Jinyuan et al. 2015). These methods were initially based on optical sensor-based snow cover. If cloud pixels exist in the optical sensor-based snow cover, these pixels are then compared with the microwave sensor-based snow cover map and replaced with snow or non-snow areas accordingly. Therefore, the accuracy of the optical sensor-based snow cover as a base map is critical.
If a snow cover map is generated using an optical sensor onboard a satellite with a sun-synchronous orbit, limitations associated with cloud cover can occur because the number of images observed each day is restricted. This problem has been discussed and mitigated by various researchers (Hall et al. 2010). The methods for generating cloud-free MODIS images are composed of spatial, temporal, and spatio-temporal algorithms without the usage of other types of sensors (Xinghua et al. 2019). However, each of the above-mentioned methods has limitations. In the case of spatial methods, the preceding conditions are complex, and the algorithm can be applied in areas with a small proportion of clouds (Dietz, Kuenzer, and Conrad 2013;Gafurov and Bárdossy 2009;Xinghua et al. 2019). With regard to temporal methods, two kinds of satellites with the same sensors (Aqua and Terra) are utilized or the values estimated from images observed on different days are conjugated (Hori et al. 2017;Xinghua et al. 2019;Mazari et al. 2013). However, these methods are less accurate if the clouds are stalled in the area for a long time (Xinghua et al. 2019). In addition, because the occurrence of abnormal weather conditions caused by Arctic amplification is becoming more frequent (Cohen et al. 2014), the probability of obtaining a clear area is lowered. The spatio-temporal method uses the aforementioned spatial and temporal methods, and it is divided into two kinds of methods, which are used step by step or simultaneously (Dariane, Khoramian, and Santi 2017;Xinghua et al. 2019;Xia et al. 2012). However, the spatio-temporal method cannot completely eliminate the effect of clouds, as it does not solve the aforementioned limitations of the spatial and temporal methods.
In the case of geostationary satellites, owing to their high temporal resolution, it is possible to significantly reduce the effect of clouds when generating snow cover maps. In particular, since the observation time of recently launched geostationary satellites is significantly shorter than that of past satellites, there is a high possibility of acquiring more surface information (de Wildt et al. 2007). However, because the geostationary satellite is located near the equator, the area where most of the snow occurs is observed with a high VZA. As mentioned above, in high VZA, the proportion of viewable snow cover obstructed by trees becomes smaller than in the nadir. In addition, a slight variation in top of atmosphere (TOA) reflectance caused by the change of solar zenith angle (SZA) over time exists. For this reason, the usage of stationary thresholds, which were generally used in snow cover mapping and cloud mask algorithms, is limited for use in geostationary satellite images.
Therefore, in this study, we propose a GK-2A multitemporal imagery-based snow cover mapping algorithm. To improve the cloud mask accuracy over a snow-covered area, the angle-time variation and the displacement by clouds (dClouds) were introduced. In the case of snow detection, pixels that were falsely detected as snow cover were removed with dynamic thresholds. In addition, the NDSI threshold was 0.1, and the ability to detect snow cover was improved through optimal stacking. This paper is organized as follows: Section 2 presents the dataset and the method of accuracy assessment, Section 3 describes the methodology of the proposed algorithm, and Section 4 provides the qualitative and quantitative validation results. Section 5 presents a discussion about the determination of the NDSI threshold and characteristics of the proposed algorithm. In Section 6, a summary and conclusion are presented.

GK-2A
GK-2A, launched in December 2018, is a next-generation geostationary satellite that succeeds meteorological observations from the Communication, Ocean, and Meteorological Satellite (COMS) in the Republic of Korea. GK-2A is operated at an altitude of about 36; 000km and 128:2 � E on the equator ). GK-2A has an advanced meteorological imager (AMI) sensor . The AMI has 16 spectral channels, including 6 reflective solar channels and 10 emissive channels. Table 1 represents the channels and spatial resolution of the AMI .
Observed images can be divided into full disk every 10 minutes, extended local area every 2 minutes, and local area every 2 minutes. We used images covering East Asia that were clipped from the full disk images with lambert conformal conic (LCC) projection, as shown in Figure 1. In Figure 1, blue points indicate the location of the automated synoptic observing system (ASOS), which is a weather station in the Republic of Korea. As an accuracy assessment, the results of the proposed algorithm are compared with the ASOS data, and the results are presented in Section 4.2.3.

Validation dataset
For the validation process, the daily snow cover products generated by Landsat-8, VIIRS, the Interactive Multisensor Snow and Ice Mapping System (IMS), and in situ data were used. All the products based on satellite images are converted to LCC projection for the generation of spatio-temporally collocated data, as shown in Figure 2.

Landsat-8
Landsat is an Earth observation satellite project jointly hosted by the United States Geological Survey (USGS) and the National Aeronautics and Space Administration (NASA). In 1972, Landsat-1 was launched, followed by Landsat-8 in 2013 and Landsat-9 in 2021. Among the Landsat series, we used the Landsat-8 image, which passes through the same spot every 16 days. Landsat-8 carries an operational land imager (OLI) with nine reflective solar channels and a thermal infrared sensor (TIRS) with two thermal infrared (TIR) channels. The spatial resolutions of the OLI and TIRS, excluding panchromatic, are 30 m and 100 m, respectively. In the case of images captured by TIRS, they are provided at 30 m resolution through interpolation. For validation, approximately 2,500 Landsat-8 images with below 30% cloud cover were used. Landsat-8 images were classified into clouds, cloud shadows, snow, and water using Fmask, which has an overall accuracy of 96.41% for the cloud class (Frantz et al. 2018). Fmask uses an object-based cloud and cloud shadow detection algorithm, and it can automatically classify clouds, cloud shadows, snow, and water for single images (Baetens, Desjardins, and Hagolle 2019;Zhu, Wang, and Woodcock 2015). We divided the Fmask results into snow and non-snow areas. For the spatio-temporal collocation of Landsat-8 and GK-2A data, the ratio of snow cover by collecting Landsat-8 snow pixels included in the 2 km spatial resolution was calculated and introduced into the image, as shown in Figure 2b.

VIIRS
The VIIRS sensor, mounted on the Suomi National Polar-orbiting Partnership (Suomi-NPP) and the National Oceanic and Atmospheric Administration-20 (NOAA-20), has 16 moderate-resolution bands, five imaging bands, and a day-night band from 0:41 μm to 12:5 μm (Cao et al. 2013). VIIRS has two spatial resolutions 375 m for imaging bands and 750m for moderate-resolution bands and can detect the surface twice each day crossing the equator. Products of VIIRS use the Sinusoidal grid tiling system, constituting 10 degrees by 10 degrees at the equator. VIIRS products are distributed by the Land Processes Distributed Active Archive Center (LP DAAC).
VNP10A1 is a daily snow cover product generated by VIIRS images, and it is used to validate the results of the proposed algorithm. VNP10A1 has a 375 m spatial resolution, and it is constructed using four imaging bands and a VIIRS cloud product. Due to the absence of the 0:5 μm channel in the imaging bands, I1 (centered at 0:645 μm) was replaced to calculate the NDSI image. For the construction of spatial-temporal collocation data, we calculated the proportion of snow in each GK-2A pixel, as shown in Figure 2c.

IMS
The IMS, which is produced by the National Oceanic and Atmospheric Administration/National Environmental Satellite, Data, and Information Service (NOAA/NESDIS), provides snow cover and snow ice maps for the Northern Hemisphere from 1997 to the present. The IMS is generated from a large variety of products, including satellite images and in situ data. The IMS utilizes polar-orbiting, geostationary, passive microwave, and synthetic aperture radar (SAR) satellite-based products, including GOES, METEOSAT, Himawari, MODIS, VIIRS, Metop, SSM/I, AMSR-2, and AMSU (NESDIS 2008;Siljamo and Hyvärinen 2011). IMS data are distributed in ASCII and GeoTiff formats at three different resolutions: 1 km, 4 km, and 24 km.
Daily IMS data in 1km distributed by the National Snow and Ice Data Center (NSIDC) were used to validate the results of the proposed algorithm. For the construction of collocation data, we resampled the resolution to 2 km, as shown in Figure 2d.

In situ data ASOS is used for ground observation in the Korea
Meteorological Administration (KMA) to understand the weather. 95 ASOSs are in operation, and they conduct observations at fixed times to acquire information on snow, clouds, visibility, ground temperature, evaporation, etc. The snow cover observed by ASOS is the actual depth of snow accumulated on the ground at the observation time, regardless of the period of accumulation.
We generated data by counting whether snow depth exists over time. In other words, if snow exists for 24 hours, a value of 24 is assigned to the corresponding observatory, and if snow exists for 12 hours, a value of 12 is assigned to the corresponding observatory. The values counted by hours were reserved according to the observatory, and the values were assigned in the image with 2 km pixels using the latitude and longitude of the observatory.

Accuracy assessment
For an accuracy assessment, we conducted qualitative and quantitative comparisons. For the qualitative comparison, we used the GK-2A natural color RGB and the snow cover products of VIIRS and IMS. In the case of the GK-2A natural color RGB, the production method provided by NMSC was utilized ( Table 2).
The natural-color RGB was designed to distinguish ice from water clouds and was constructed by applying histogram equalization to the visible channel (0:64 μm) and NIR channels (0:86 μm, 1:6 μm). In the natural color RGB, the snow-covered area and ice clouds are shown in cyan. For this reason, the pattern of an object needs to be considered if ice clouds are to be recognized.
For quantitative comparison, we used the confusion matrix indicated in Figure 3. We calculated the precision, recall, overall accuracy, and F1 score (Equations (1)-(4)) as validation scores, as follows: Precision represents how many snow pixels in the proposed algorithm correspond to the target snow pixels. Recall denotes whether the snow pixels were precisely captured without being dropped. Overall accuracy is an indicator of how accurately the proposed algorithm detects snow. The F1 score, which is calculated using the harmonic mean, denotes the difference between precision and recall. If the difference between these two values increases, the F1 score also increases. In addition, temporal validation was conducted to compare the SCE between the results of the proposed algorithm and the other snow cover maps. On a daily basis, we calculated correlation coefficients and compared trend differences between each snow cover product.

Methodology
The daily snow cover mapping algorithm proposed in this study is based on GK-2A multitemporal images. Figure 4 indicates a workflow of the entire process of the proposed algorithm, including image preprocessing. The proposed algorithm is applied to the GK-2A images, and a snow cover map at 10-minute intervals is first generated. Snow cover maps with 10 mins are accumulated, and the number of detected snowcovered pixels is calculated. After pixels with misclassification are removed, the daily snow cover map can be obtained.
where Ref is the TOA reflectance in a specific area with i, j pixel coordinates, and Rad denotes the radiance given by (a 1 �DN þ a 2 Þ. a 1 ; a 2 indicate the slope and offset, respectively, which are provided by NMSC. d is the distance (astronomical units) between the Earth and sun. SZA is the solar zenith angle, and ESUN is the mean solar exoatmospheric irradiance, which is calculated by integrating the solar spectral irradiance (SSIÞ and the spectral response function SRF ð Þ. In the channels with longwave radiance, DN can be calculated as radiance using the calibration table provided by NMSC. Specifically, the slope (a 1 ) and offset (a 2 ) are provided in the calibration table. Equations (7) and (8) are used for BT calculation.
where T e indicates the effective brightness temperature (k), k 1 and k 2 can be calculated with light speed, Planck's constant, and the Boltzmann constant, and the values are 1.1910E-16 and 1.4388E-02, respectively. v means wavenumber (cm À 1 ). In addition, the coefficients (c 0 , c 1 , c 2 ) for conversion to BT (T b ) are presented from the calibration table.

Discrimination between snow and clouds
In shortwave radiation (0:3-2 μm), fresh snow has a high albedo, from 75% to 95%. As time passes after the snow falls, the albedo decreases by 40%-60% for permanent snow, and the albedo of snow that becomes ice or soiled decreases (Geiger, Aron, and Todhunter 2010). The variables that affect snow albedo include crystal size, density, amount of dust, surface roughness, snow-water content, and the number of thaws (Geiger, Aron, and Todhunter 2010). Although the reflectance of snow and clouds is different enough to be distinguished, the difference is reduced due to the variables mentioned above. Therefore, the accuracy of the utilized cloud detection algorithm is important for finding the exact snow-covered pixel.
In the case of a geostationary satellite, images are captured under a wide range of illumination conditions (de Wildt et al. 2007). For this reason, there are changes in spectral information according to time changes. If a threshold-based cloud detection algorithm is used for geostationary satellite images, it can be difficult to identify pixels with cloud and snow. Therefore, we employed filtering techniques that comprised the angle-time variation and the discrepancy between the target and minimum TOA reflectance to detect the change in spectral information of the surface. We modified filtering techniques by Lee and Choi (2021) to be suitable for snow areas. The pixels removed using the discrimination between snow and cloud are indicated by red in Figure 8a.

Elimination of thick clouds using angle-time variation
Reflectance refers to the ratio of reflected radiance to incident solar flux, as indicated in Equation (5). If the effect of the SZA is corrected and the effect of atmospheric components is weak, reflectance under a clear sky is nearly unchanged. If the reflectance is expressed as an angle according to time change, the angle in a clear sky should be approximately 45 � . However, when clouds flow into a pixel, a considerable change occurs in the pixel. The reflectance of clouds is significantly higher than that of the surface, except for snow areas in the visible channel. For this reason, angle calculated by reflectance is also significantly changed. Therefore, Lee and Choi (2021) developed a method based on angle variation calculated by using TOA reflectance (0:8 μm).
In this study, NDSI is used as a supplementary variable, and the equation used by Lee and Choi (2021) is modified to diminish the misclassification. The NDSI can be calculated using Equation (9). In the visible and near-infrared wavelength, the reflectance of clouds and snow-covered areas is similarly high. Nonetheless, angle variation exists because the reflectance changes depending on the cloud thickness. In the case of SWIR (1:6 μm), because absorption by ice particles exists, different aspects of the wavelength shorter than SWIR can be found. Therefore, we supplemented the NDSI to increase cloud detection accuracy in snow-covered areas.

NDSI ¼
Ref 0 where Ref 0:51 and Ref 1:6 mean TOA reflectance at 0:51 μm and 1:6 μm from AMI. The area where SZA exceeds 80 � is masked so that only the daytime image is used. Thus, the presence of reflectance in images varied according to time changes in the area near SZA 80 � , and this leads to the incorrect calculation of angles. Therefore, the backward equation was additionally used in the original equation, as follows: where θ indicates the angle at observation time (t) and Δt denotes the time step, which we set to 10 minutes. D indicates TOA reflectance (0:8 μm) or NDSI. Figure 5 presents the variation in the TOA reflectance (0:8 μm) and the NDSI according to the observation time. In Figure 5a, the red and blue lines (dots) indicate changes in TOA reflectance and NDSI under clear-sky conditions, respectively. The orange and cyan lines (dots) show the fluctuating values when clouds are introduced into the pixels. In the case of Figure 5b, the graphs represent the angles calculated from Figure 5a using Equation (10). As shown in the red and blue lines in Figure 5b, the angle hardly changes in a clear sky. However, clouds flow into the pixel, and the angles of TOA reflectance and NDSI fluctuate (orange and cyan). In other words, the range of angles in the cloudy sky is larger than that in the clear sky. In addition, the inflection points of angles in the snow-covered area under cloudy conditions are dissimilar. In the case of NDSI, it reflects the characteristics of 1:6 μm that can distinguish ice crystals from water clouds. That is, in visible channels, ice particles at the top of convective clouds appear bright, whereas in the 1:6 μm, the reflectance is lowered due to absorption. Therefore, there are different aspects between the angles of the TOA reflectance and the NDSI. A supplementary explanation of these different aspects is added to the description of Figure 6. In the case of non-snow areas, negative NDSI values can be acquired, as shown in the bottom graph of Figure 5a. When a negative number is used in Equation (10), the negative angle can be obtained. If the NDSI is close to or slightly exceeds zero due to an influx of clouds, an angle is calculated around −100 or as close to 0 as the cyan line in the bottom graph of Figure 5b. Because the trend of angle variation is unspecific, the use of NDSI for areas without snow is not proper. In the case of NDSI, if the angle is less than 44 � or exceeds 46 � in the area where the NDSI value is over 0, it is considered to be clouds, and the pixel is removed. In the case of TOA reflectance, if the angle is calculated between 42 � and 48 � , it is considered a clear sky. We determined the range of thresholds through empirical testing.
The different aspects between angles calculated by TOA reflectance (0:8 μm) and NDSI can be found in Figure 6. The stratocumulus is on the west side of Figure 6b, and the cirrus is on the east. All the clouds are over the snow-covered area. In the case of the angle calculated by TOA reflectance, the area where the angle variation exists represents the location of the clouds. However, in general, the number of pixels that appear as extreme values (30 � or 60 � in Figure 6d) is limited. In the case of the angle calculated by NDSI in Figure 6f, the stratocumulus located in the west can be clearly captured. However, the range of angles in the areas where the cirrus clouds exist is not wide. Therefore, use in conjunction with angles calculated by TOA reflectance and NDSI is effective.

Semitransparent and stalled cloud mask using discrepancy in TOA reflectance
If clouds consistently flow into a specific pixel or if semitransparent clouds exist, the method using angletime variance is difficult to be utilized because the changes in TOA reflectance or NDSI are not considerable. In other words, an angle can be calculated within a predefined range of thresholds, although the clouds flow over the pixels. In order to solve the problem, Lee and Choi (2021) utilized the minimum TOA reflectance. Minimum TOA reflectance is a method for acquiring surface reflectance by calculating the minimum value among reflectance values observed in the same area during a specific period. For cloud masks with minimum TOA reflectance, a certain threshold that is pre-calculated between the target and minimum TOA reflectance is needed. Since cloudy pixels need to be fully separated from clear sky observation, Lee and Choi (2021) developed filtering technique modified by the cloth simulation filtering (CSF) method (Zhang et al. 2016). CSF was originally developed to determine the shape of terrain; it is a digital surface model (DSM) that operates by dropping a virtual cloth down on an inverted (upside-down) point cloud (Zhang et al. 2016). Constant reflectance in a clear sky must be maintained so that the pre-calculated threshold works properly. However, in some snow-covered areas, we found that the reflectance decreased as the SZA increased. Due to this reduction, some areas covered by snow or desert are removed after the filtering technique is applied.
In snow-covered areas, when the solar and satellite zenith angles are on the same line, the backscattering value appears to peak. As the solar zenith angle increases, forward scattering appears because of the anisotropic shape of the snow Dozier 2004a, 2004b). These changes are also seen in the desert because desert soils have relatively opaque vertical structures (Wang et al. 2005). Therefore, the variations in TOA reflectance are large in snow or desert areas where the mixed land cover proportion is low. For this reason, we improved the method for accurately detecting semitransparent and stalled clouds over snow-covered areas without masking the clearsky pixel by additionally using the NDSI and the average inverse 0:8 μm TOA reflectance as follows: where X i;j;t indicates the dClouds and inv 0:8 μm ð Þ denotes the inverse TOA reflectance of 0:8 μm, given by 1 À Ref ð Þ�10. inv min ð Þ is the inverse minimum TOA reflectance calculated by accumulating all images taken over a seven-day period.
Because the TOA reflectance of bright surfaces decreases when aerosol concentrations increase (Kim, Kim, and Yoon 2014), we calculated the minimum TOA reflectance using the solar zenith angle at 5 � intervals. In the case of snow-covered areas under clear sky, the inverse TOA reflectance and the displacement by gravity calculated with the original version of Equation (11) have values lower than the inverse minimum TOA reflectance, as shown in Figure 7a. Therefore, through the usage of NDSI and the average values of images, the dClouds have a higher value than the inverse minimum TOA reflectance. In other conditions, such as non-snow under cloudy sky, the dClouds have the same pattern as the value calculated in the original version as orange points in Figure 7. In Equation (11), NDSI is activated as a weight function that can be adjusted depending on the amount of snow in the pixel. If the NDSI is close to 0, the weight is reduced. Accordingly, the usage probability of the threshold can be increased so that clouds in the snow-covered area are classified. However, the range of NDSI is from −1 to 1, and the maximum value of the inverse data is 10.  Therefore, the average value is used to amplify the NDSI. The average value is calculated using all the pixels in images (row: n, column: m) acquired during a day (UTC 0 to k).
After dClouds (X i;j;t ) is calculated in images to which angle-time variation is applied, the discrepancy (discr i;j;t ) between the inverse minimum TOA reflectance and X i;j;t is utilized.
We determined the threshold of discr i;j;t in Equation (12) using the Otsu method. The Otsu method can discover the threshold at which the variance between classes is maximized after obtaining a histogram of an image and dividing it into classes (Otsu 1979). We calculated the average threshold for two months (1 December 2020 to 31 January 2021) and obtained a value of 5.7184. If the discr i;j;t exceeds this threshold, the pixel is determined to be clouds and is removed.

Residual cloud mask
The cirrus channel at 1:37 μm (from 1:37 μm to 1:38 μm) is located in strong water vapor absorption and is observed through radiance mainly reflected by the troposphere because upward radiance from the surface or the lower-middle atmosphere is blocked by water vapor. If a particular threshold in a non-snow area is utilized to detect cirrus clouds in snow-covered areas, snow pixels with low NDSI values can be removed. Therefore, in snow-covered areas, the threshold is slightly higher than that used for other land cover type (Frey et al. 2008). In this study, we calculated thresholds for each image applied with the cloud masks proposed in Section 3.2.1 and Section 3.2.2 using the Otsu method.
In addition, to eliminate cloudy pixels left over after the cloud masking steps above, we supplemented the BTD (3:7 μm-10:5 μm). The BTD is based on the absorption difference between these two wavelengths resulting from both water and ice cloud particles (Frey et al. 2008). Therefore, the BTD has been used in several studies related to cloud masks (Allen, Durkee, and Wash 1990;de Wildt et al. 2007;Frey et al. 2008;Zhang et al. 2020). In addition, the BTD has been shown to effectively improve the performance of cloud masks for snow mapping (Hall et al. 2001;Yeom, Han, and Lee 2009). Because 3:7 μm detects both the emitted radiation and the small amount of reflective radiation, we set the threshold by calculating the three-sigma value (plus the standard deviation) at every observation time, with consideration of the change in SZA (de Wildt et al. 2007;Derrien and Le Gléau 2005). A pixel with a value higher than the threshold was considered to be clouds and removed. The abovementioned masking process proceeded in order.

Daily snow cover map generation
Clouds are distinguished from snow and non-snow areas in Section 3.2. However, due to the characteristics of the optical sensor, if clouds persist in a specific area, detection of snow-covered areas is rarely accomplished. Therefore, the process of accumulating the snow-covered area was performed, and the snow pixels above a certain threshold were used to generate a daily snow cover map after the misclassified pixels were removed.

Removal of misclassified pixels
The channel of 13:3 μm was used to correct for false snow mapping in tropical forests because the NDSI value is sensitive to small changes in dark, dense forest areas and remaining clouds (Barton, Hall Dorothy, and Riggs George 2020). The presence of aerosols over the canopy and the water area under the canopy can cause misclassification (Barton, Hall Dorothy, and Riggs George 2020). Therefore, the misclassified pixels in the image were eliminated using the three-sigma rule.
In addition, the low TOA reflectance (0:8 μm) in the remaining NDSI pixels after performing the abovementioned processes was removed because low reflectance can be associated with cloud shadows, surface intrinsic properties, and topographic shadows (Hall et al. 2001). The threshold was obtained by applying the Otsu algorithm to values less than or equal to 0.3 below latitudes of 40 � , where these misclassifications mainly occurred. The 0.3 was determined to prevent the threshold calculated by the Otsu algorithm from rising and was obtained through empirical testing. If the number of extracted pixels was small, the Otsu algorithm was omitted.

Optimal stacking of snow cover map
After the masking process for clouds and misclassified pixels was accomplished, NDSI values greater than 0.1 were considered snow-covered areas, as shown in . We set the NDSI threshold as 0.1 is to improve the accuracy of snow detection in areas with a high proportion of forests. Theoretically, snow exists if NDSI is greater than 0. However, in the NDSI range from 0 to 0.1, uncertain snow detections or snow commission errors were common (Riggs, Hall, and Román 2016). Therefore, we defined an NDSI below 0.1 as a non-snow area.
To generate a daily snow cover map, we stacked the snow-covered pixels observed for three days before the target date to minimize the missing snowcovered area caused by clouds. When the pixels were stacked, a problem arose: pixels with clouds or false detection that microscopically existed in each NDSI image were gathered and formed into a chunk of pixel. For this reason, stacked pixels below three were removed. In other words, among the pixels that accumulated for three days, only pixels that stacked more than three were used as snow. The thresholds concerning stacking days and the number of accumulated pixels were determined through quantitative and qualitative analysis. As shown in Figure 8b, through the above process, it was possible to minimize the effect of clouds and discriminate clear snow-covered areas.

Qualitative comparison
For qualitative comparison, the GK-2A natural-color RGB and snow cover products, which are VIIRS and IMS, were utilized. To understand the visual strengths and weaknesses of the proposed algorithm, the four days from December 2020 to January 2021 were selected and compared.
On December 20, as shown in Figures 9 (a)-(d), the snow-covered pixels observed inland (40 � N À 50 � N, 80 � E À 140 � E) had similar trends between the result of the proposed algorithm and the VIIRS. In the north of Uvs Lake (50 � N, 92 � E) in Mongolia, there was a difference between the IMS and the optical sensorbased products. The continuous inflow of clouds in the area was confirmed in the GK-2A natural-color RGB. In addition, Japan and the Kamchatka Peninsula (55 � N, 159 � E) represented the difference from IMS due to the continuous inflow of clouds. In all these areas, because clouds continued to flow in for more than three days, the proposed algorithm could not detect the snow. In the vicinity of Ganzhou (26 � N, 115 � E), the snow pixels of VIIRS were judged as misclassifications when compared with GK-2A natural-color RGB and IMS. The misclassification is caused by discrimination between clouds and snow because VIIRS observes the area at a specific time.
On December 29, as shown in Figure 10d, the area near Nanjing (32 � N, 119 � E) was classified as snow in the IMS. However, the optical sensor-based snow cover map (Figures 10 (a)-(c)) could not detect the snow in the area due to the presence of stalled clouds. When the clouds cleared in this area on December 30, the proposed algorithm detected snow in the area (Figure 10e). In addition, when compared with the VIIRS (Figure 10g), the misclassification in southern China did not appear in the results of the proposed algorithm.
On January 23, as shown in Figures 11 (a)-(d), most of the snow-covered areas were detected when comparing the results of the proposed algorithm to the IMS (Figures 11 (a)-11(d)). In the case of VIIRS, there were areas that could not be detected due to the inflow of clouds, especially in northern Mongolia (47 � N, 102 � E). When IMS was compared with GK-2A natural color RGB and VIIRS (Figures 11  (b)-(d)), there were misclassified pixels in the vicinity of Torey Lake (50 � N, 115:5 � E). It was confirmed that the IMS was classified as a snow-covered area near the relevant area.

Validation with Landsat-8
There is a large difference in spatial resolution between the GK-2A (2 km) and Landsat-8 (30 m) snow cover maps. Because GK-2A acquires spectral information mixed from various land covers (soil, snow, vegetation, and impervious areas) within 2 km, the dissimilarity in snow-covered pixels with Landsat-8 unavoidably exists. Therefore, we found a starting point where high accuracy remained constant after the proportion of snow in the pixels was divided at 5% intervals. If 50% was set as a starting point, a pixel with a proportion of snow over 50% was classified as snow, while a pixel below 50% was classified as a non-snow area. The lowest O.A. was 0.78 at 0%, and the lowest F1 score was 0.70 at 95% as the median value. If 45% was used as the starting point, the highest validation scores, such as 0.92 for O.A. and 0.89 for the F1 score, could be obtained. However, the ratio at or above 30% was stable because the differences in the median values between 30% and 45% were 0.007 for O.A., 0.015 for precision, 0.022 for recall, and 0.003 for the F1 score. Therefore, a Landsat-8 snow cover map applied with a 30% ratio was utilized to calculate accuracy, as shown in Figure 12.
The outlier lower than the lower whisker in Figure 12 arose on 21 December 2020. The snowcovered area not detected by the proposed algorithm mainly occurred in the central inland region of China at 30 � N À 40 � N and 110 � E. The undetected snow-covered pixels, which mainly occurred in urban or mountainous areas, had low NDSI values, from 0 to 0.1. The range of low NDSI values was not considered snow in this study because of uncertain snow detection or snow commission errors, as mentioned above.

Validation with the IMS and VIIRS
IMS and VIIRS snow cover products were used to validate the proposed algorithm. The period of data used for validation was two months (1 December 2020-31 January 2021). We calculated validation scores using daily snow cover maps and expressed them as a boxplot, as shown in Figure 13.
In the case of VIIRS, the accuracy assessment was performed similarly to the procedure conducted in Section 4.2.1. When the proportion of snow in a pixel exceeded 25%, the best median F1 score of 0.7348 was obtained. The worst results were found at 95%, and the differences between 25% and 95% in median O.A., precision, recall, and F1 score were 0.023, 0.264, 0.128, and 0.123, respectively. Therefore, we used the snow cover map applied with a 25% ratio, which showed the best accuracy, as shown in the blue boxplots in Figure 13.
In the case of IMS, most GK-2A snow-covered pixels were included in the IMS pixels, with a median precision of 0.9743. However, the median recall value was 0.6668 due to some areas not being detected as clouds.
When the products of IMS and VIIRS were compared with the GK-2A result, the O.A. values were similar. However, the F1 score, which was compared with the IMS product, had a slightly high value. Given that misclassification exists in the VIIRS product, a comparison between the IMS and the proposed algorithm may yield more reasonable results. Therefore, since the median values of the O.A. and F1 scores were around 0.8 in comparison with the IMS product, as shown in Figure 13, the effectiveness of the proposed algorithm is confirmed.

Validation with in situ data
For validation, the snow data observed by ASOS in the Republic of Korea were used, as shown in Figure 1. The cases where snow was detected for more than 12 hours and where snow was detected only for 24 hours at the observatory were compared. In the case of more than 12 hours, the proposed algorithm has a precision close to 0.9 as shown in Table 3. This means that the proposed algorithm matched the pixels classified as snow at a rate close to 90%. However, the recall was 0.62, which is about 38% of non-detection, by the proposed algorithm. In the case of IMS, which is fused with various satellites, all validation scores presented higher values than the proposed algorithm. However, when only snow cover detected for 24 hours was used, the proposed algorithm presented higher validation scores than IMS, such as + 0.14 for precision, + 0.043 for the F1 score, and + 0.12 for O.A. In the case of VIIRS, all validation scores were lower than the proposed algorithm. As a result, the proposed algorithm is more effective in detecting snow cover than VIIRS. In the case of IMS, the detected snowcovered area may be overestimated because the precision of 0.56 over 24 hours is low. This accords with the results of qualitative comparison with IMS.
Because the topography of the Republic of Korea is not homogeneous, the spatial distribution of snowfall is different from region to region (Lee and Lee 2006), e.g. in the west coastal area of the Republic of Korea, when the Siberian High expands, snowfall occurs due to the influence of the ocean (Lee and Lee 2006). In the case of the inland area located in the southwest of  Quantitative validation results of snow cover maps generated by the proposed algorithm and Landsat-8 imagery. The period of data used for validation is two months (1 December 2020-31 January 2021). The validation scores were calculated using daily snow cover maps and expressed as boxplots. Figure 13. Quantitative validation results of the proposed algorithm with the IMS and VIIRS snow cover maps. The period of data used for validation is two months (1 December 2020-31 January 2021). The validation scores were calculated using daily snow cover maps and expressed as boxplots.
Republic of Korea, when the Siberian High is strengthened and its influence reaches inland, snowfall occurs (Lee and Lee 2006). In the case of the east coastal area of the Republic of Korea, snowfall is determined by the inflow of northeasterly winds. Therefore, when in situ data are compared with snow cover products with a 2 km spatial resolution, the accuracy tends to be low.

Temporal validation with the VIIRS and IMS
The SCE was calculated with daily snow cover maps based on the proposed algorithm (red), the IMS (green), and the VIIRS (yellow), as shown in Figure 14. As shown in the IMS line (green), the SCE decreased to 6; 211; 228 km 2 in mid-December 2020 and increased to 7; 970; 908 km 2 at the end of January 2021. Because the correlation between the SCEs of GK-2A and IMS was 0.9126, these products showed similar trends. The average difference between the GK-2A and IMS was 2; 244; 825 � 368; 549 km 2 , and a constant difference between the two datasets existed. The reason for this constant difference is the presence of undetected cloud pixels due to the limitation of the optical sensor, as mentioned in the qualitative comparison. In addition, the IMS contains overestimated snow-covered pixels in comparison with the VIIRS data in the area over 50 � N, as shown in Figures 10 and 11. In the case of VIIRS, the correlation with GK-2A is 0.4922, while the correlation with IMS is 0.3236. When the SCE decreased in mid-December 2020 and increased in the second half of January 2021, the VIIRS snow cover did not follow the trend of IMS. In other words, the VIIRS SCE shows a uniform trend without changes in slopes.

Comparison with methods in forest areas
In this study, an NDSI value greater than 0.1 was utilized to detect snow pixels in the image. To check the usefulness of the NDSI threshold determined in this study, it was compared with the method proposed by Klein, Hall, and Riggs (1998), which is used in the MODIS and VIIRS snow cover map algorithms. In forested areas, the forest canopy blocks the radiation emitted from snow-covered areas, and the mixed spectral data are thus acquired on sensors (Wang et al. 2018). For this reason, NDSI values are low in forested  areas because the spectral reflectance in the visible wavelength is much lower than that in the snow, and the spectral reflectance in the mid-infrared wavelength increases due to the influence of the soil and forest canopy (Klein, Hall, and Riggs 1998;Wang et al. 2018). In addition, varying grain sizes cause the greatest NDSI changes at the lowest canopy cover percentages, and increasing the illumination angle decreases the NDSI (Klein, Hall, and Riggs 1998). Therefore, an adjustment composed of the NDSI with the normalized difference vegetation index (NDVI) was suggested (Klein, Hall, and Riggs 1998). However, when the method proposed by Klein, Hall, and Riggs (1998) was used, a pixel with snow-covered areas was misclassified as an area without snow (Raleigh et al. 2013).
We compared the accuracy of the results generated by the proposed algorithm and the method of Klein, Hall, and Riggs (1998) in the same environment. In other words, the processes of cloud mask and misclassification removal that were developed in this study were utilized in the method of Klein, Hall, and Riggs (1998), and NDSI pixels with a range of 0.1 to 0.4 were applied with the linear regression equations as follows (Key et al. 2013): where a 1 and a 2 are the coefficients. Klein, Hall, and Riggs (1998) calculated the coefficients through the field boundary constructed by the GeoSAIL model and the Landsat TM observation. GeoSAIL is a relatively simple reflectance model for discontinuous canopies (Klein, Hall, and Riggs 1998). The model is used for the estimation of crown reflectance and areal proportions to generate the total stand reflectance. If the NDSI value from 0.1 to 0.4 falls in the boundary, it is considered snow. In this study, because the use of the GeoSAIL model is limited, the field boundary was generated using GK-2A images and the MODIS land use land cover (LULC) dataset. Among the LULC data, only the deciduous needleleaf tree areas in the plant functional type (PFTs) were used. We built the linear regression equation by constructing the field boundary. We set the two points, which are the low NDSI point when NDVI was 0.1 and an NDSI value of 0.4 when NDVI was 0.1, by referring to the figures of Klein, Hall, and Riggs (1998). After the NDSI and NDVI values Figure 15. Comparison of the accuracy between the results of the proposed algorithm and the results obtained by Klein, Hall, and Riggs (1998) method. The period of data used for validation is two months (1 December 2020-31 January 2021). The validation scores were calculated using daily snow cover maps and expressed as boxplots.
included in the two points were extracted, we applied the 1 sigma rule to the dataset and created a linear regression equation.
The data for two months from 31 December 2020 to 31 January 2021 were used for comparison between the proposed algorithm and the method of Klein, Hall, and Riggs (1998). As shown in Figure 15, the validation scores were calculated and expressed as boxplots. The differences in validation scores were as follows: 0.054 for O.A., −0.014 for precision, 0.146 for recall, and 0.08 for the F1 score. Although a slight increase in precision was found, the remaining three validation scores of Klein, Hall, and Riggs (1998) method were all lower than those of the proposed algorithm. In addition, the interquartile range and the lower whisker of the boxplot generated by Klein, Hall, and Riggs (1998) method are longer than the boxplot generated by the proposed algorithm. The number of pixels classified as snow-covered areas in Klein, Hall, and Riggs (1998) method was smaller than those of the proposed algorithm in the Far Eastern Federal District of Russia and near Harbin city in China. The pixels that are not included in the linear regression equation are likely snow-covered pixels. The result indicates the same as the opinion asserted by Raleigh et al. (2013). Therefore, we found that the proposed algorithm performed better than Klein, Hall, and Riggs (1998) method, as shown in Figure 15.

Strengths and weaknesses of the proposed algorithm
As a result of the accuracy assessment with Landsat-8-based snow cover maps, the median F1 score was 0.89. Additionally, the proposed algorithm can detect snow-covered pixels with these F1 scores when the proportion of snow in the pixel is at least 30%. The probability of false detection was lower than that of the VIIRS-based snow cover mapping algorithm. In addition, the trend of SCE has a correlation of 0.9 with the IMS product.
However, as mentioned above, the result of the proposed algorithm has a problem caused by clouds. We utilized the data accumulated over three days to minimize the effect of clouds. If the clouds in the area persist for three days, it becomes difficult to detect snow because the surface information cannot be obtained. Therefore, this limitation may be solved through a fusion process with a microwave sensorbased snow cover product.

Summary and conclusion
We proposed a daily snow cover mapping algorithm based on GK-2A multitemporal images. The steps of the algorithm involve (1) cloud masking if considerable changes in TOA reflectance and NDSI exist, (2) masking the remaining cloud pixels using the dClouds proposed in this study and by using the cirrus channel and BTD (3:7 μm-10:5 μm), (3) extracting NDSI over 0.1 and accumulating the NDSI images observed at 10-min intervals to create daily snow cover maps. We evaluated whether the proposed algorithm was well suited to discriminate snow and other land cover types and validated the results using Landsat-8, IMS, VIIRS, and in situ snow cover products. Although the snow cover mapping algorithm obtained satisfactory results, some limitations occurred due to the characteristics of the optical sensor when clouds stalled in a specific area for a long time. Therefore, in the future, we will try to develop a blending algorithm to solve the problems associated with undetected snow-covered pixels.