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.
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