What do values from 100 to 200 mean in MOD10A1 NDSI snow cover layer?

Question

I am working with MODIS snow cover products (MOD10A1) and am unable to understand some of the values that are returned. I am trying to get % snow cover from the NDSI (normalised difference snow index) snow cover layer. The MODIS user manual states that the NDSI snow cover layer has values from 0 - 100, representing % snow cover in each pixel, and eight values between 200 and 255 that represent all other possible features/masks (cloud, missing data etc). In processing the images I am finding values between 100 and 200 and cannot find any reference to such values in the MODIS documentation.

I downloaded the MOD10A1 product as .hd files from the NSIDC.org site. I work in R, but have not been able to work with the .hd files in R, so I converted the NDSI snow cover layer to .tif files using the HEG converter program that is recommended on the MODIS NASA website. I imported the .tif files into RStudio using the raster package and used the getValues and unique functions to find the value in each pixel. The returned values are anything from 0 to 255, including values in the range 100-200.

Does anyone know what these values mean? Have they come with the product or is there an error in the file conversion? Thank you for your help.

EDIT: thanks for the suggestion. One of the exact file names is 'MOD10A1.A2015364.h25v06.006.2016182181418.hdf' and a link to the file https://drive.google.com/file/d/1HeEpIL15EC_PSBWsuGT4FJMZOPr4_oND/view?usp=sharing

I have tried using the rast function in the terra package and get the same result.

Answer 1

You can look at this with terra.

library(terra)
f <- "MOD10A1.A2015364.h25v06.006.2016182181418.hdf"

You can do

x <- rast(f)
names(x)
#[1] "MOD_Grid_Snow_500m:NDSI_Snow_Cover"                    "MOD_Grid_Snow_500m:NDSI_Snow_Cover_Basic_QA"          
#[3] "MOD_Grid_Snow_500m:NDSI_Snow_Cover_Algorithm_Flags_QA" "MOD_Grid_Snow_500m:NDSI"                              
#[5] "MOD_Grid_Snow_500m:Snow_Albedo_Daily_Tile"             "MOD_Grid_Snow_500m:orbit_pnt"                         
#[7] "MOD_Grid_Snow_500m:granule_pnt"                       
 

And use the first layer

ndsisc <- x[[1]]
ndsisc
#class       : SpatRaster 
#dimensions  : 2400, 2400, 1  (nrow, ncol, nlyr)
#resolution  : 463.3127, 463.3127  (x, y)
#extent      : 7783654, 8895604, 2223901, 3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source      : MOD10A1.A2015364.h25v06.006.2016182181418.hdf:MOD_Grid_Snow_500m:NDSI_Snow_Cover 
#names       : MOD_Grid_Snow_500m:NDSI_Snow_Cover 

Or acknowledge the sub-dataset structure of the file

s <- sds(f)
s
#class       : SpatDataSet 
#subdatasets : 7 
#dimensions  : 2400, 2400 (nrow, ncol)
#nlyr        : 1, 1, 1, 1, 1, 1, 1 
#resolution  : 463.3127, 463.3127  (x, y)
#extent      : 7783654, 8895604, 2223901, 3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source(s)   : MOD10A1.A2015364.h25v06.006.2016182181418.hdf 
#names       : NDSI_Snow_Cover, NDSI_Snow_Cover_Basic_QA, NDSI_Snow_Cover_Algorithm_Flags_QA, NDSI, Snow_Albedo_Daily_Tile, orbit_pnt, granule_pnt 

names(s)
#[1] "NDSI_Snow_Cover"                    "NDSI_Snow_Cover_Basic_QA"           "NDSI_Snow_Cover_Algorithm_Flags_QA"
#[4] "NDSI"                               "Snow_Albedo_Daily_Tile"             "orbit_pnt"                         
[7] "granule_pnt"  

And use the first sub-dataset

ndsisc <- s[1]
ndsisc
#class       : SpatRaster 
#dimensions  : 2400, 2400, 1  (nrow, ncol, nlyr)
#resolution  : 463.3127, 463.3127  (x, y)
#extent      : 7783654, 8895604, 2223901, 3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source      : MOD10A1.A2015364.h25v06.006.2016182181418.hdf:MOD_Grid_Snow_500m:NDSI_Snow_Cover 
#names       : NDSI snow cover from best observation of the day 
 

Now check the values with freq (or unique)

fq <- freq(ndsisc)
tail(fq)
#      layer value  count
#[84,]     1    92     28
#[85,]     1    93     15
#[86,]     1    94      3
#[87,]     1   201  33636
#[88,]     1   237  37178
#[89,]     1   250 807785

There are values below 100 and above 200, but not in between, as you assert.

A next step could be

nd <- clamp(ndsisc, 0, 100, values=FALSE)
nd
#class       : SpatRaster 
#dimensions  : 2400, 2400, 1  (nrow, ncol, nlyr)
#resolution  : 463.3127, 463.3127  (x, y)
#extent      : 7783654, 8895604, 2223901, 3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source      : memory 
#names       : NDSI snow cover from best observation of the day 
#min values  :                                                0 
#max values  :                                               94  

Perhaps you where accidentally using the "Snow_Albedo_Daily_Tile" subdataset?

albedo <- s[5]
albedo
#class       : SpatRaster 
#dimensions  : 2400, 2400, 1  (nrow, ncol, nlyr)
#resolution  : 463.3127, 463.3127  (x, y)
#extent      : 7783654, 8895604, 2223901, 3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source      : MOD10A1.A2015364.h25v06.006.2016182181418.hdf:MOD_Grid_Snow_500m:Snow_Albedo_Daily_Tile 
#names       : Snow albedo of the corresponding snow cover observation 

That one has values between 100 and 200. But is that unexpected?

as.vector(unique(albedo))
#[1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30
#[32]  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61
#[63]  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  88  91  92  93  94  95
#[94] 100 101 125 137 150 251 253