Continuity of MODIS and VIIRS Snow Cover Extent Data Products for Development of an Earth Science Data Record

An Earth Observing System global snow cover extent data products record at moderate spatial resolution (375–500 m) began in February 2000 with the Moderate-resolution Imaging Spectroradiometer (MODIS) instrument onboard the Terra satellite. The record continued with the Aqua MODIS in July 2002, the Suomi-National Polar Platform (S-NPP) Visible Infrared Imaging Radiometer Suite (VIIRS) in January 2012 and continues with the Joint Polar Satellite System-1 (JPSS-1) VIIRS, launched in November of 2017. The objective of this work is to develop a snow cover extent Earth Science Data Record (ESDR) using different satellites, sensors and algorithms. There are many issues to understand when data from different algorithms and sensors are used over a decade-scale time period to create a continuous dataset. Issues may also arise with sensor degradation and even differences in sensor band locations. In this paper we describe development of an ESDR derived from existing MODIS and VIIRS data products and demonstrate continuity among the products. The MODIS and VIIRS snow cover detection algorithms produce very similar daily snow cover maps, with 90–97% agreement in snow cover extent (SCE) in different landscapes. Differences in SCE between products ranged from 2–15% and are attributable to convolved factors of viewing geometry, pixel spread across a scan and time of observation. Compared at a common grid size of 1 km, there is a mean of 95% agreement in SCE and a difference range of 1–10% between the MODIS and VIIRS SCE maps. Mapping sensor observations to a coarser resolution grid reduces the effect of the factors convolved in the 500 m tile to tile comparisons. We conclude that the MODIS and VIIRS SCE data products are reliable constituents of a moderate-resolution ESDR.


Introduction
Snow cover extent (SCE) has been cited as an Essential Climate Variable in the Global Climate Observing System because of its tremendous importance in the Earth's energy balance and as a source of freshwater for more than 1 billion people worldwide [1]. Snow is also important globally for recreation and for municipal water-supply needs. Snow can cover up to 50% of the Northern Hemisphere land surface in winter, covering 1.9 to 45 million km 2 [2]. SCE is thus an important variable needed to address science questions and Earth science applications objectives, including those concerning societal needs [2].
In the Northern Hemisphere SCE is highly variable from year to year, thus decade-scale trends may be difficult to identify until a sufficient length of record is available. Knowledge of onset, duration and melt date of snow cover is critical to understand interannual and longer-term changes at regional to hemisphere scales. Changes in regional SCE are both an indicator and a result of climate change.
Continuity of the daily gridded snow cover extent products from the Moderate-resolution Imaging Spectroradiometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS) is presented here from the perspective of the products being used to investigate snow cover parameters of seasonal snow cover onset, duration and melt.
A 54-year climate-data record (CDR) of SCE in the Northern Hemisphere has been developed by the National Oceanic and Atmospheric Administration (NOAA) and Rutgers University [3]. This 25 km resolution CDR has been used to determine that spring snow cover is melting earlier, and that in many areas, the duration of the snow season is decreasing [4]. For example, Arctic spring snow cover extent is declining by approximately 18% per decade [5]. While global trends in SCE have been captured in the NOAA/Rutgers SCE CDR, the requirements for an SCE Essential Climate Variables product are that it be at a daily frequency and approximately 1-km resolution [1].
The capability now exists to develop a multi-decade, daily, moderate resolution data record of SCE using multiple spacecraft platforms and sensors. Data from the MODIS instruments on the Terra and Aqua satellites, combined with data from the VIIRS Suomi-National Polar Platform (S-NPP) and the Joint Polar Satellite System (JPSS) satellites may be used to create an Earth Science Data Record (ESDR) that will become a CDR when a sufficient number of years of data become available. The record begins in February 2000 and continues to the present and will be extended into the foreseeable future. The Terra and Aqua MODIS sensors are expected to continue nominal operations into 2022 when planned end of life operations will begin. However, depending on sensor and platform status and budget resources, data collection may continue through 2024. JPSS-1, now operational and designated as NOAA-20, will be followed by JPSS-2, JPSS-3 and JPSS-4, providing continuous coverage through 2038. These instruments will potentially provide a 40-year record of satellite snow cover observations.
A characteristic of a satellite-based CDR is that the record can be reprocessed as needed to take advantage of improvements in algorithms, software "fixes" to compensate for sensor degradation, orbit changes, differences between sensors and other factors. Each reprocessing, which creates a new data collection, must include provenance to document the changes.
The NASA VIIRS Land Science Team has worked to ensure that the VIIRS land and cryosphere data products have continuity with the corresponding MODIS land and cryosphere data products, because the MODIS products have been validated and used extensively by both the science and applications communities [6][7][8]. Numerous investigators have used the standard MODIS products to investigate snow cover parameters (e.g., [9][10][11][12][13][14][15][16]).
The objective of this work is to develop an ESDR using the Terra, Aqua and S-NPP daily SCE products, extending from 2000 to the present (and ultimately beyond). Development of an ESDR involves documentation of issues relating to use of different satellites, sensors and algorithms, to develop a continuous and well documented data record. Cross comparison methods were used to compare SCE among the MODIS and VIIRS daily SCE products. In this paper, we provide examples from the western United States of product continuity along with uncertainty information, both at native product resolution and with the products re-projected to a 1 km resolution grid. Documentation of data product continuity on the global scale, though the ultimate goal, is beyond the scope of the present paper.

MODIS and VIIRS Instruments and Satellites
The Terra, Aqua, S-NPP and NOAA-20 satellites are in near polar orbits, with Terra in a descending orbit with a morning equator crossing time and Aqua, S-NPP and NOAA-20 in ascending orbits with afternoon equator crossing times (Table 1). Characteristics of the sensors relevant to the snow cover detection algorithm based on the normalized-difference snow index (NDSI) are listed in Table 1. The NOAA-20 characteristics are included in Table 1 because of their relevance to continuing the SCE data product record beyond S-NPP. MODIS and VIIRS have similar visible and near infrared spectral bands but have different spatial resolutions. Both are whiskbroom instruments so the pixel spatial resolution grows larger away from nadir. The growth in pixel size is not controlled in MODIS but for VIIRS a unique on board processing technique is employed to limit the growth in pixel size so that pixels at the edges of scans are similar in size to those near nadir [17,18] while MODIS allows for scan overlap. That difference in sampling strategy is apparent in the Level 1B (L1B) products where MODIS retains the overlap pixels; the same features appear in successive lines at far edges of a swath where-as in VIIRS those overlap regions are trimmed on board and the L1B swath product has lines of bow tie trim data in the overlap regions. In L2G processing a gridding algorithm removes the overlap regions in MODIS products and the bow tie trim data in the VIIRS data products.

MODIS and VIIRS Snow Cover Products
The MODIS and VIIRS snow cover products are generated as a series from the instrument swath Level-2 (L2) to daily gridded and projected products at Level-3 (L3). The snow cover detection algorithm is applied at L2, then the L2 products are gridded to tiles of the sinusoidal projection and observations are composited to produce the L3 products. The L3 daily gridded snow cover products of Terra MOD10A1, Aqua MYD10A1 [19] and VIIRS S-NPP VNP10A1 [20] are used in this study. The MODIS products are provided at 500 m resolution in HDF-EOS2 format, and the VIIRS product is provided at 375 m resolution in HDF-EO5 format, and includes netCDF CF-1.6 conventions for attributes and mapping data to a projection. The snow cover extent data products are available to download from the National Snow and Ice Data Center (NSIDC) Distributed Active Archive Center (DAAC).
The data products are named using the Earth Science Data Type (ESDT) convention of ESDT.AYYYYDDD.HnnVnn.CCC.YYYYDDDHHMMSS.ext. The ESDT is MOD10A1, MYD10A1 or VNP10A1, and the date of acquisition is year (AYYYY) and day (DDD). Location of the tile on the sinusoidal projection is given by .Hnn the horizontal tile index and .Vnn the vertical tile index with the origin of the grid, (0, 0) at −180 • and 90 • ; the Collection denoting reprocessing is given as CCC. When referring to a Collection the convention of Cn, e.g. C6 is used. Product production date is YYYYDDDHHMMSS and product format is given by the .ext. The MODIS and S-NPP data products and user guides are available from the NSIDC DAAC: The snow cover data products are validated based on the Committee on Earth Observation Satellites validation stages 0-4 [21]. The MODIS snow cover products have been validated to Stage 2 [22] and the VIIRS snow cover product has been validated to Stage 1 [23]. Validation case studies summarized on the VIIRS Land website [23] report 97% agreement between MODIS and VIIRS SCE, while MODIS Aqua and VIIRS S-NPP products are found to have an average agreement of ±0.05-0.10 between NDSI values of collocated pixels.
The NOAA-20 snow cover product VJ110A1 is currently in science production/testing and planned for release in Collection 2, possibly in late 2020. Evaluation of the VJ110A1 snow cover data product is not included in this study because the product is currently not publicly available.
In addition to the L3 data products, new cloud-gap-filled MODIS and VIIRS data products are to be released in MODIS Collectin 6.1 (C6.1) and VIIRS Collection 2 (C2) data collections respectively. An evaluation of these products is presented in [24]. The daily L3 snow cover products are input to the cloud-gap-filled product algorithm so the continuity evaluation here is also relevant to the cloud-gap-filled products.

Snow Cover Detection Algorithm
A ratio or normalized difference ratio of bands in the visible (VIS) and the near-infrared (NIR) parts of the spectrum to detect snow has been used extensively to detect snow from aircraft and satellite instruments [25]. The NDSI is effective because snow has a high reflectance in the VIS and a low reflectance in the NIR. The NDSI is calculated from the difference in VIS and NIR top-of-atmosphere (TOA) reflectance data: where: ρ = top-of-atmosphere reflectance. The NDSI is a robust parameter for detection of snow cover which has been used with many VIS and NIR sensors, e.g. Landsat, MODIS and VIIRS, and implemented in different ways [26][27][28][29] with consistently good results. Snow cover is present when the NDSI is 0.0 ≤ NDSI ≤ 1.00. The snow cover algorithm limits NDSI to the 0.10 to 1.0 range due to uncertainty associated with low NDSI values [19,20]. For the snow cover data products the NDSI is scaled by 100 and stored in the 0-100 range.
The MODIS and VIIRS instruments have similar spectral bands that can be used to calculate the NDSI. For MODIS the NDSI is calculated using band B4 centered at 0.55 µm and band B6 centered at 1.66 µm. For VIIRS the NDSI is calculated using band I1 centered at 0.64 µm and band I3 centered at 1.61 µm. Though the VIS bands used for NDSI from the MODIS and VIIRS are different they are both in the VIS region where snow has characteristically-high reflectance as shown in Figure 1. The MODIS and VIIRS NIR bands have a similar central location and bandwidth. The relative spectral response functions for MODIS B4 and B6 [30] and VIIRS I1 and I3 [31] are shown in Figure 1. The snow reflectance curve in Figure 1 was created using a field spectral reflectance dataset collected with a mast mounted FieldSpec Pro JR spectroradiometer by ASD (Boulder, CO, USA) available from [32] and discussed in [33]. It is derived from an average of three dry snow samples collected at the same site on 17, 18 and 22 March 2010 when environmental conditions were similar [32].
Snow cover is highly reflective across all the visible wavelengths ( Figure 1). Thus the detection of snow in the visible part of the spectrum is relatively insensitive to which visible band is used, though it is preferable to use red and green wavelengths because of lesser atmospheric effects than in the blue wavelengths.
The snow cover detection algorithm is applied at the swath level (L2). The cloud mask product (described in Section 4.3) is used as input to the snow cover detection algorithm to mask clouds. The L2 snow cover products are then processed through a series of L2G gridding algorithms and intermediate products where the observations from the L2 swath products are gridded and projected on to the sinusoidal projection. In regions where swaths overlap, all of the observations covering a grid cell are stacked, and then an observation for the day is selected based on criteria related to solar zenith angle, location relative to nadir and spatial coverage in a grid cell. The observation of the day is stored in the L3 products. The algorithms and products are described in the MODIS [34] and VIIRS [20] user guides. The plots of relative spectral response for MODIS and VIIRS bands were created using data from [30,31] respectively. The spectral reflectance curve for snow was derived from an average of field spectral reflectance data collected at the same site with similar environmental conditions on 17, 18 and 22 March 2010 as described in the text and available from [32].
Observation data on sensor scan angle and viewing geometry are not passed to the L3 product, however time of the acquisition of the observation (the time of the swath acquisition) is passed to the L3 product. The L2G and L3 algorithms are the same for processing/selection of MODIS and VIIRS observations; the algorithms differ in gridding between the 500 m and 375 m resolution versions of the sinusoidal projection.

Study Region
Our study region is the western United States. Two tiles, outlined in black in Figure 2, of the MODIS and VIIRS daily products that include mountains and plains were used for the tile to tile comparisons. The region used for the 1 km gridded SCE comparisons is covered by six tiles of the products outlined by the orange rectangle in Figure 2. Map of tiles covering most of North America on the sinusoidal projection. Tiles h08v05 and h10v04, used for tile to tile comparisons are outlined in black, The SCE 1 km grid region used data from these six tiles h08v04, h09v04, h10v04, h08v05, h09v05, and h10v05 outlined in orange.
The MODIS to S-NPP VIIRS comparisons were done for 16 and 17 February 2012. In 2012 snow cover extent in western North American was below average in February [35]. The study region as seen on 16 February 2012 in Worldview Terra MODIS true color imagery is shown in Figure 3 with states and geographic features identified. Snow covers mountain ranges and some areas of the plains. These two dates were chosen because cloud cover was minimal and snow cover was near seasonal maximum. This image ( Figure 3) provides a geographic frame of reference for the results and discussion. Mountains are the focus of discussion regarding tile h08v05 covering parts of California, Nevada and Arizona and plains are the focus of discussion regarding tile h10v04 covering parts of eastern Montana, Wyoming, South Dakota and Nebraska.

Methodology
Cross comparison methods have been used to investigate continuity between MODIS and VIIRS land data products. The use of cross comparison methods on land data products generated from ratio based algorithms was reviewed to assess whether it is reasonable to apply similar methods to the snow cover products from MODIS and VIIRS. The Normalized Difference Vegetation Index (NDVI) and other ratio based vegetation indices are similar to the NDSI. Products derived from the NDVI and other indices were cross compared at a climate modeling grid (CMG) scale of~5 km spatial resolution to minimize effects that factors such as sensor characteristics and viewing geometry can have on cross comparisons of observations at sensor pixel resolution. Notable results from cross comparison studies of MODIS and VIIRS NDVI [36,37], enhanced vegetation products [38], vegetation indices [39] and the bidirectional reflectance distribution function and albedo products [40] studies show that despite convolved factors there is strong continuity between the products and that an ESDR created from MODIS and VIIRS can be used to monitor vegetation changes such as seasonal greenness.
For the present study, two cross-comparison analyses were conducted; the first compared two tiles of the products and the second compared a 1 km SCE grid covering the western USA constructed from six tiles of the products (Figure 2). The 1 km grid was created to compare SCE at a CMG scale. Both comparisons were done for 16 and 17 February 2012. Both Terra and Aqua are compared with the S-NPP though Terra is in a different orbit, and convolves several factors in the comparison.
To enable visual comparison of the SCE among the products the NDSI_Snow_Cover data layers were converted into un-projected colorized images to show the snow cover extent and other features such as clouds and lakes that are reported as flag_values in a data layer. To facilitate cross comparisons of cell-to-cell observations in the tiled products, the VNP10A1 product at 375 m resolution was resampled to the 500 m grid resolution of MOD10A1 and MYD10A1 using Geospatial Data Abstraction library (GDAL) functions. Additionally, the product file formats were converted to GeoTIFF in this processing step. In these cross comparisons the NDSI_Snow_Cover data layer from the products is used. Data layers were compared based on data array dimensions.
To enable comparisons of the products at a resolution that reduces the effects of sensor and orbit differences, we created a 1 km grid SCE product for the western United States. Products from the six daily tiles (Figure 2), were re-projected from the sinusoidal projection onto the SCE CMG 1 km grid in geographic projection using GDAL functions. Cell-to-cell comparisons of SCE were made for three pairs-Terra to Aqua, Aqua to S-NPP and Terra to S-NPP-using only cloud free observations in both products of a pair to map SCE.
For results presented in Table 2 and Table 3 the difference in SCE in paired comparisons is expressed as a percentage of SCE cell counts. A cell is considered snow covered if the NDSI value is in the 10-100 range and not snow covered if the NDSI = 0. The count of cells that are snow covered in both products is considered agreement and is used as the base to calculate the percentages of differences between the pair of products. The percent of difference is calculated as the count of SCE cells in only one product of the pair divided by the count of cells of SCE in both products. Cloud observations and flag_values were excluded from comparison.
NDSI relationships among the products are investigated with two-dimensional histograms (density plots) of paired observations and paired mean NDSI differences and standard deviations. The 2D density plots ( Figure 8, Figure 9, Figure 10 and Figure 11) reveal relationships in the data that may not be apparent in a scatter plot of the variables. For the 2D density plots the frequency of occurrence is given as a relative scale, from low to high to enable use of the same scale for all comparisons. Paired mean NDSI differences were calculated as the mean difference of pairs of observations for NDSI values of the sensor on the x-axis, with the difference on the y-axis being the mean sensor difference of (x sensor-y sensor), the solid line and +/− one standard deviation plotted as the dashed lines. The differences were binned using a bin size interval of two NDSI values, i.e., 10-11, 12-13, and so forth then the mean and standard deviation were calculated as shown in Figure 8, Figure 9, Figure 10 and Figure 11.
We used the NASA EOSDIS Worldview app (https://worldview.earthdata.nasa.gov/) for interactive browsing of MODIS and VIIRS imagery to assist in selecting mostly cloud free dates and comparing orbit tracks of satellites.

MODIS to VIIRS Comparison
To assess the continuity of SCE between products, cross comparisons of the products at native resolution and common 500 m resolution were made. Images of the NDSI_Snow_Cover data layer for MOD10A1, MYD10A1 and VNP10A1 for two days, 16 and 17 February 2012, two tiles covering different geographic regions were compared to assess data continuity among the products as shown in Location of snow cover is very similar in all three products (Figures 4 and 5) but the extent of snow cover varies around the peripheries of snow covered regions because of convolved factors such as acquisition time, viewing geometry, pixel size change across a scan, and projecting and gridding the observations. Large areas of cloud cover are similar among the three products on a given day but small areas of cloud and cloud cover at the periphery of SCE areas differ due to differences in the cloud mask algorithms for MODIS and VIIRS, time of acquisition and other convolved factors. The SCE for NDSI in the 10-100 range, is in reasonably good agreement among the images based on visual interpretation. The viewable SCE is consistent from day to day.
Cell-to-cell comparison of the MODIS and VIIRS NDSI_Snow_Cover using the S-NPP (VNP10A1) products resampled to 500 m resolution, projected onto the geographic grid and converted to GeoTIF format are shown in Figures 6 and 7. These comparisons were made for the h08v05 and h10v04 tiles as described above (Figures 4 and 5). In these paired comparisons NDSI_Snow_Cover values in the 10-100 range were interpreted as snow cover, so that if both products were mapped as snow cover they agree, and if not they disagree. The purpose of these SCE maps (Figures 6 and 7) is to clearly show where products agree and disagree on SCE under cloud free viewing conditions.    Results of cross comparisons of the SCE maps in Figures 6 and 7 quantifying differences in SCE between products are presented in Table 2. Cross-comparisons of the products as array images at native resolution (Figures 4 and 5) and as gridded products in which the sensor pixels observations have been projected and gridded (Figures 6 and 7), do not reveal the many factors that can affect pixel-to-pixel comparisons (e.g., [7,41,42]). The MODIS and VIIRS instruments each have different pixel resolutions at nadir and at the edges of a scan. For both sensors the area of a pixel increases towards the edges of a scan, i.e., the bow tie effect, which is processed with different scan sampling strategies as described in Section 2.1. The VIIRS pixel spread across a scan is about half of that which occurs in MODIS. The location of an observation within an instrument scan influences where the observation will be mapped in a grid cell. This can influence cell-to-cell comparisons. The sensors and platforms have many differences (Table 1) that affect observations such as: time of day, viewing geometry, and whether the orbit is ascending or descending. The snow cover detection algorithm uses TOA reflectance so atmospheric effects are convolved with other factors. Snow cover extent is very similar in all the products, with most of the differences appearing at the periphery of regions of common SCE (Figures 4-7) ( Table 2). Terra shows about 15% greater SCE compared to Aqua or S-NPP, and Aqua shows about 14% greater SCE compared to S-NPP in the mountainous regions on both days. These differences are expected due to differences in orbits, nadir tracks, viewing geometry, solar illumination, and pixel size growth across the terrain.
There is much less difference in SCE among the products over the plains landscape versus mountainous terrain on both days ( Figure 5) and in the calculated differences of the gridded products in Figure 7 and Figure 9 where differences can range from 1-8% (Table 2). This is due to the generally flat topography of the plains. Convolved factors appear to be less important over the plains than over the mountainous landscape. In the comparison maps of tile h10v04 ( Figure 7) the majority of the observed snow cover is seen as common in all three comparisons with differences in snow cover scattered around the periphery of the common snow. In all of the paired comparisons (Figure 7) differences in the location of snow cover, outside the common snow cover extent, are scattered over the landscape. In the mountainous regions the differences in snow cover are concentrated around the periphery of snow covered regions as seen in Figure 6. Such differences in SCE may be associated with changes in cloud cover or snow cover, or both, between data acquisition times.
In some situations, differences in SCE between products can be attributable to rapidly changing snow cover conditions. In the comparison on 16 February in Figures 5 and 7, there are differences in SCE in the plains in northeastern Montana. Those differences were investigated using the Worldview app for Terra and Aqua imagery and data products. (The S-NPP imagery and snow cover data product were not available in Worldview to assist in the investigation.) In that region on 16 February there is considerable difference in appearance and extent of snow cover and cloud cover between the morning overpass of Terra and the afternoon overpass of Aqua because some of the snow cover observed in the Terra imagery had melted by the afternoon. However, there is more snow cover in S-NPP which is more similar to Terra vs the Aqua SCE; this suggests that there may be other factors affecting the comparison with Aqua that need to be considered.
Comparisons of SCE between products within a day and between days are affected by location of a region in the swaths of the satellite. For example, the snow cover observed in the mountains of northwestern Arizona is in a different location in the Terra and Aqua swaths ( Figure 6). In Terra, the region is near nadir of the orbit track and in Aqua the region is on the far western edge of the orbit track on 16 February then on 17 February the region is on the far west of the Terra orbit swath and in the east of the Aqua swath thus confounding paired comparisons. (The S-NPP imagery and orbit tracks were not available in Worldview to help with this comparison.) Satellite orbit tracks vary from day to day with a near repeat cycle of every eight days, so a ground location may be found in different positions in a swath, i.e. east, center or west in the swath, from day to day.
Some of the differences in snow cover detection between sensors are likely caused by cloud/snow confusion because the three cloud mask products for Terra, Aqua and VIIRS are different for each sensor as described earlier. The primary concern is that clouds not detected as certain cloud could then be detected as snow in the L2 snow algorithm because the snow cover algorithms all mask cloud based on the certain cloud flag from a cloud mask product. Cloud/snow confusion that occurs in L2 is propagated into the L3 snow cover data products.
The SCE products are consistent for the observed common regions of snow cover; no large unexpected or unexplained differences in SCE are observed in any of the comparisons. The SCE maps are also consistent with visual interpretation of imagery on a regional scale.
To explore NDSI relationships between products, 2D histograms of paired NDSI observations and paired mean NDSI differences were calculated and are shown in Figures 8-11 for the tile to tile comparisons.   The 2D histograms and mean differences of NDSI in Figures 8-11 reveal relationships between products. The mean difference in NDSI values between paired product comparisons was determined by calculating the mean difference of pairs of observations for NDSI values of the sensor on the x-axis, with the difference on the y-axis being the mean sensor difference (x sensor-y sensor); the solid line and +/− one standard deviation line are plotted as dashed lines.
In general, for all the sensor comparisons there is good agreement between NDSI values over the range of 40-80 where the 2D histograms have the densest concentrations of data surrounding the 1:1 line in the paired cell-to-cell comparisons across both mountain (tile ho8v05) and plains (tile h10v04) landscapes on both dates (Figures 8-11). The NDSI histograms (not shown) also show the peak of distribution in the NDSI range of 40-80, with a tail of low counts in the 10-40 range and a sharp drop of NDSI values ≥ 90 reflecting the distribution of points in the scatter plots. There is more scatter in the NDSI relationships over the 10-40 range, which is the middle range of frequency of occurrence density of points in all the comparisons. The greatest differences between NDSI values occur far from the 1:1 line and occur in the lowest range of frequency of occurrence (Figures 8-11). The NDSI mean and difference plots (Figures 8-11) show relatively consistent differences in all the comparisons, across the full range of NDSI values though in some comparisons there is a spike in differences at NDSI values ≥ 90. Terra to Aqua NDSI comparison results are similar for both dates and landscapes. The 2D histograms (Figures 8-11) show a dense concentration of data in the 40-80 range. Scatter increases at NDSI < 40, and the amount of scatter seems related to the landscape with lesser amount of scatter on the plains (Figures 9 and 11) compared to the mountain landscape of tile h08v05 on both dates (Figures 8 and 10). The pattern of middle to low frequency of occurrence is similar in all the NDSI comparisons, indicating consistency among the sensors. In the NDSI comparisons the density pattern of highest frequency is generally elongated along and near the 1:1 line with slight variation in the pattern. The density pattern on 16 February for the mountain dominated tile h08v05 (Figure 8) has less elongation of the highest density region. In this mountain landscape the mean Terra-Aqua NDSI difference ranged from −10 at low NDSI values to +10 at high NDSI values with a spike at the highest values. Those results suggest that factors convolved in the product comparisons had a relatively consistent effect among the products.
The convolved factors in comparing Terra and Aqua NDSI values did not appear to have a significant effect on the observations in these comparisons albeit they may be a factor on 16 February. The MYD10A1 (Aqua) tile h08v05 on 16 February is composed of two swaths of L2 data mapped into the tile, with each swath covering about half the tile in a diagonal line north to south thus contributing to greater differences in acquisition time and viewing geometry convolved in the comparisons. On 17 February the tile was mostly covered by a single Terra orbit track and Aqua and S-NPP single orbit coverage and orbit tracks close together and near nadir coverage. The location of an orbit track over a region has an effect not only on the pixel size, as discussed earlier, but also on the geolocation for mapping a pixel observation to the surface and projection to the grid. A pixel may thus be larger in mountainous terrain. The spikes in mean NDSI differences at about NDSI = 90 on 16 and 17 February for tile h08v05 in Figures 8 and 10 respectively, for all three comparisons have a very low population of observations. Thus, further investigation of them was not pursued for this paired comparison.
Aqua to S-NPP NDSI comparisons show high frequency of occurrence in the NDSI 40 to 90 range (Figures 8-11) on both dates and tiles with S-NPP NDSI values being greater than MYD10A1 in the NDSI mean difference. The NDSI mean differences were in the −10 to 0 range with the difference trending from −10 at low NDSI values to 0 at about NDSI = 30 then steady to about NDSI = 90 where there was a spike in differences with MYD10A1 values being much greater than S-NPP. The greatest differences were found for tile h08v05 (primarily mountainous landscape (c.f. Figure 3)) on 16 February where the NDSI differences ranged from −20 to 0 across the range of NDSI 10-90 ( Figure 8). On that date Aqua has two swaths from different orbits mapped into the tile, with each swath covering about half the tile in a diagonal line from north to south thus greater differences in acquisition time and viewing geometry from east to west convolved in the comparisons to the single orbit acquisition of S-NPP. In all of the comparisons there was a spike in the differences at Aqua values of NDSI ≥ 90 (Figures 8-11) where Aqua was greater than S-NPP. There are very few observations in that NDSI range. There were approximately 1700 Aqua cells with an NDSI value ≥ 85 which is about 0.6 % of the total of SCE cells in the tile, which we consider inconsequential.
Terra to S-NPP comparisons show a similar spread in NDSI in the 2D histograms (Figures 8-11) however, the mean differences are in the range of NDSI −10 to 0 with S-NPP tending to have slightly greater values. The greatest spread in NDSI differences is seen on 16 February for tile h08v05 covering a mostly mountain landscape. Over the plains as seen in tile h10v04, the NDSI comparisons have a lesser amount of spread on 17 February compared to 16 February, suggesting that day to day differences in NDSI values may be influenced by differing viewing conditions and changes in snow cover.
The spike in mean NDSI difference at Terra NDSI ≥ 85 in tile h08v05 on both days is the result of a small number of NDSI observations that are ≥85 in the valid range of NDSI values, consisting of about 0.4% on 16 February and 0.02% on 17 February. This is reasonable considering the many factors convolved in the comparison.
The important result of these cross comparisons done at the 500 m tile level resolution is that Terra, Aqua and S-NPP snow cover products all detected very similar SCE across the landscape. Differences in NDSI values in cross comparison of cells have similar patterns of frequency of occurrence as seen in Figures 8-11. The amount of difference, tends to be small and constant or to have a gradual trend across the range of NDSI values. Exceptions may occur at NDSI values ≥ 90. The NDSI values are expected to vary between sensors, within a day and day to day due to the many convolved factors that affect the NDSI. The similarity of results for all the gridded SCE comparisons (Figures 6 and 7) and the lack of any observed abrupt changes or unexpected discontinuities indicate excellent continuity of SCE between the products.

Comparison of 1 km Gridded SCE
The convolved factors in satellite cross comparison studies can be minimized by resampling or compositing the data to a coarser spatial resolution when comparing data products (c.f. [37,39]). To enable comparisons of the products at a resolution that reduces the effects of sensor and orbit differences, we created a 1 km grid SCE product for the western US. Six tiles of the data products ( Figure 2) were re-projected from their native sinusoidal projection to the geographic projection at 1 km resolution. The 1 km resolution SCE region is bounded by 35 • -50 • N, 105 • -125 • W ( Figure 12). The study area is 2,487,610 km 2 in area. The six tiles cover a slightly larger area (bounded by 30 • -50 • N, 90 • -130 • W) than the 1 km SCE study region (c.f. Figure 3).
Cell-to-cell comparisons of SCE were made for three pairs-Terra to Aqua, Aqua to S-NPP and Terra to S-NPP using only cloud free observations in both products of a pair. The SCE maps of the western U.S. for 16 and 17 February 2012 are shown in Figure 12 and the statistics of the comparisons are reported in Table 3.  Table 3. Paired SCE CMG comparison statistics for the western United States. Clear only observations in both products were compared. The percentages are based on the count of SCE cells in one product and snow free land (SCE = 0), and in the other as a percentage of snow cover common to both products. The SCE regions are consistent among these three paired 1 km gridded SCE comparisons with the common SCE shown in white in Figure 12. The differences in SCE between the sensors are presented as a percentage of SCE as shown in Table 3. Differences in SCE between days are due primarily to changes in cloud cover from day to day. Same day differences in SCE can be caused by melting snow cover, or by differing cloud conditions or viewing geometry of swath acquisitions. The Terra SCE tends to be 5-10% greater than either Aqua or S-NPP in some regions as depicted in blue in the Terra comparisons in Figure 12. The causes of those differences are differences in orbit tracks, viewing of a region near nadir or near the edge of a swath, and cloud masking. As an example, the region where more Terra snow cover (blue) is seen on 16 February in northwestern California (Figure 12), is the result of a combination of the fact that the region was imaged at the far eastern edge of a swath and some thin cloud resulting in confusion with snow in the morning. When imaged by Aqua in the afternoon the sky was clear and the same region was in the central part of the swath. Much of the thin cloud in that region was not flagged as certain cloud in Terra and that cloud was erroneously detected as snow in the L2 algorithm and product which then appears in the L3 SCE product. The differences between Terra and S-NPP can be attributed to similar factors.

Paired
The Aqua to S-NPP paired comparisons ( Figure 12) are very similar in extent of snow cover with only~5-6% difference ( Table 3). The snow cover mapped by Aqua only (blue) occurs at the periphery of common snow regions and is attributable to the large spread in pixel size as compared to the limited growth of VIIRS pixels across a scan. The SCE mapped by S-NPP only (red) also occurs at the periphery of common snow cover, but the cells are more scattered on the landscape and are fewer in number.
Results at the 1 km grid scale show greater similarities in SCE and lesser amounts of single sensor snow cover outside of the common SCE regions than were observed in the 500 m tile comparisons. This is demonstrated by the lower amounts of SCE only for a product in Table 3 versus Table 2. In the tile to tile comparisons the SCE only differences ranged from 1-16% and in the 1 km grid comparisons the SCE only differences ranged from 2-10%. This supports the assertion that mapping the sensor observations to a coarser resolution grid reduces the effect of the factors convolved in the 500 m tile to tile comparisons.

Cloud Masking
Though the topic of snow/cloud confusion is not specifically addressed in this paper, we present a brief discussion to highlight potential cloud/snow confusion problems, plans for the cloud mask sensor specific products and the cloud continuity product.
The snow cover detection algorithm uses the MODIS Terra or Aqua cloud mask [43][44][45][46][47] or the VIIRS cloud mask [48] product to mask clouds in a scene. We have found these cloud mask products to have generally good MODIS to VIIRS data continuity for the purpose of cloud masking in the snow cover algorithms. However a special concern exists with MODIS Aqua in which a majority of band 6 (1.6 µm) detectors are non-functional. Band 6 is the NIR band used in the Terra cloud mask and snow-cover algorithms. The Aqua cloud mask algorithm was modified to use NIR band 7 (2.13 µm) in place of Aqua's non-functional band 6 and this renders cloud detection in the Terra and Aqua cloud masking products somewhat different. See [49] for more details on cloud masking of the MODIS and VIIRS SCE products. As a result we observe increased amounts of cloud cover in some situations in Aqua as compared to Terra. The Aqua snow cover algorithm uses the Quantitative Image Restoration (QIR) algorithm [50] to restore band 6 data for use in the Aqua MODIS algorithm.
The cloud mask algorithm developers use a continuity algorithm that is the same for both MODIS and VIIRS to generate a new MODIS-VIIRS cloud mask (MVCM) product which is based on the MODIS cloud mask algorithm and uses only the bands that are common to both sensors [48]. To adapt the algorithm to use fewer bands (there are fewer VIIRS bands than MODIS bands) the MVCM algorithm includes some new spectral tests and different test thresholds that are not used in the MODIS algorithm. Validation of the MVCM has demonstrated cloud detection agreement at the pixel scale to be in the range of ±2% as compared to the MODIS Aqua cloud detection and Cloud-Aerosol Lidar with Orthogonal Polarization instrument cloud detection [48]. Based on that limited validation research the MVCM appears to provide consistent cloud masking between MODIS and VIIRS, though more extensive validation studies are needed to firmly establish confidence in the MVCM.
While evaluation of the MVCM is continuing, the cloud mask products generated in the MODIS Adaptive Processing System and Land Science Investigator-led Processing System are used in the snow cover algorithm to produce products in the MODIS C6.1 and VIIRS C2 data collections, respectively. In the snow products we have observed overall consistency in cloud masking albeit with minor differences. If it is decided that the MVCM should be used in the future, then a reprocessing of MODIS and VIIRS data collections would be indicated.

Conclusions
In this paper we show data product continuity among the NASA Terra MODIS, Aqua MODIS and S-NPP VIIRS SCE standard data products. Continuity is essential to develop an ESDR of SCE at moderate spatial resolution.
The MODIS and VIIRS snow cover detection algorithms are very similar in performance, producing very similar daily snow cover maps with 90-97% agreement in SCE when compared across different landscapes at a grid cell resolution of 500 m. Differences in SCE between the products ranged from 2-15% and are attributable to convolved factors of viewing geometry, pixel spread across a scan and time of observation. When compared at the SCE CMG resolution of 1 km, to mitigate the effect of convolved factors, there was a mean of 95% agreement in SCE and a difference range of 1-10% between the MODIS and VIIRS SCE maps. No significant random or systematic differences in SCE among the products were observed. This supports the assertion that mapping sensor observations to a coarser resolution CMG reduces the effect of the factors convolved in comparing higher resolution products from different sensors.
Cell to cell comparisons of NDSI values were generally close to the 1:1 line, though in some situations the concentration of data was large and spread out. NDSI values from the different sensors were most similar in the NDSI range of 40-80, which is also the range where the distributions of observations are greatest. Differences in NDSI values at <20 and >90, at the tails of the NDSI distributions, increased in some situations but represented a relatively small number of observations. The mean NDSI differences between observations (Figures 8-11) tended to be near zero except near the tails of the NDSI range where some scenes had upward or downward swings. No significant random or systematic differences in NDSI values among sensors were observed. This analysis demonstrates that there is very good continuity among the MODIS Terra, MODIS Aqua and VIIRS S-NPP SCE data products both within a day and from day to day, making the products suitable for use in development of an ESDR covering the MODIS and VIIRS eras. When at least three decades of the data record become available the data product record can be reevaluated for suitability as a climate data record.