Let it snow! Development of an operational snow cover product from Sentinel-2 and Landsat-8 data

=> "Winter is coming" ― George R.R. Martin, A Game of Thrones As Christmas holidays are approaching you might want to know if there is snow in your favorite spot of ski touring? A good knowledge of the snow cover variability is important - not only to plan your next week-end, but also because the snow is a key water resource in many regions, including here in south west France. The winter snowpack in the Pyrenees is like a natural water reservoir, releasing meltwater runoff during the spring, when the crops in the cultivated lowlands need to be irrigated. The snow melt is also extensively used for hydropower production. We are developing a snow cover product from Landsat-8 and Sentinel-2 data to provide the snow presence or absence at 20 m resolution every 5 days. The algorithm is simple because the snow surface is quite straightforward to detect from high resolution optical imagery. A challenge is typically to avoid the confusion between the snow cover and the clouds. Here the job was made easy since we take as input the level 2A product so we can use the cloud mask from the MACCS algorithm as a prior information. Another great advantage of the level 2A product is that it provides slope-corrected surface reflectances. Did I tell you that mountain regions are quite sloping? The slope correction enables to use the same detection thresholds in south-facing and north-facing slopes. Actually the main unresolved issue is the obstruction of the ground surface by the forest canopy, which can hinder the snow detection in evergreen forest areas. Unfortunately MACCS does not yet provide bottom of canopy reflectances... 

Overview of the algorithm

 We want a robust and efficient algorithm to process large mountain regions with a reasonable computation cost. Since the snow cover extent is largely controlled by the elevation we incorporated a digital elevation model in the input data to avoid false positive in the low elevation areas and false negative in the high elevation areas. A first pass of snow detection is performed with strict thresholds on the NDSI and red reflectance values ($N_1=0.4, R_1=0.2$) in order to retrieve a set of pixels for which we have a very high confidence that they are well classified as snow. A pixel is marked as snow covered if: $\textrm{snow}_1 = \textrm{cloud}_1 \textrm{ is false and } \textrm{NDSI} \gt N_1 \textrm{ and } \rho_\textrm{Red} > R_1$     (Pass 1) However, we do not take the original L2A cloud mask as input $\textrm{cloud}_1$ here because it is processed at a coarser resolution, very conservative and much useful information for snow cover mapping is lost. We allow the reclassification of some cloud pixels in snow or no-snow only if they have a rather low reflectance. We select only these potentially "dark clouds" because the NDSI test is robust to the snow/cloud confusion in this case. The cloud flag is removed for these pixels unless they were flagged as a cloud shadow. $\textrm{cloud}_1 = (\textrm{cloud}_{L2A} \textrm{ and } \rho_\textrm{Red}>R_c) \textrm{ or } \textrm{"cloud shadow" is true}$ Then we estimate the lowest elevation of the snow cover in the image ($Z_S$) using the SRTM digital elevation model.