Count the number of days per year with snow - Google Earth Engine

Question

I want to generate a raster for each year, where the pixel value refers to the number of days per year where there is some snow on the ground. In this example I want to create this raster for Germany and the year range is [2015,2017], which means I will get 2 images. I use NDSI for this. The code I wrote is below; but I feel like the number of days is underestimated. At some areas where I expect to have > 90 days of snow, I have about 60. I would then like to get the average number of days across Germany, i.e. for each year, take the number of days with snow from each pixel and average it.

Any solution?

var startDate = ee.Date.fromYMD(2020, 1, 1)
var endDate = ee.Date.fromYMD(2021, 1, 1) // Exclusive
var dataset = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');
var aoi = dataset.filter(ee.Filter.eq('country_na', 'Germany'));

var collection = ee.ImageCollection('MODIS/006/MOD10A1')
  .filterBounds(aoi)
print(aoi)

var numberOfYears = endDate.difference(startDate, 'years').floor()
print(numberOfYears)
var daysSnowCollection = ee.ImageCollection(
  ee.List.sequence(0, numberOfYears.subtract(1))
    .map(daysSnow)
)

print(daysSnowCollection)
Map.addLayer(daysSnowCollection.first(), {min: 0, max: 100, palette: 'red,yellow,blue,white'}, 'First image')

//var output = daysSnowCollection.reduce(ee.Reducer.mean())
//print(output)

function daysSnow(yearOffset) {
  var date = startDate.advance(yearOffset, 'years')
  return collection
    .select('NDSI_Snow_Cover')
    .filterDate(date, date.advance(1, 'years'))
    .map(function (image) {
      var snow = image.gte(30)
      return image.updateMask(snow)
    })
    .count()
    .unmask(0)
    .clip(aoi)
    .set('date', date.format('yyyy-MM-dd'))
}

Answer 1

While you have an image for every day, some pixels are masked. Some of those masked pixels might actually be snow days. This could be one of the reasons for your underestimation - you're simply counting the number of unmasked snow days.

An alternative to this is to calculate the fraction of snow pixels and multiply with the number of days of the year. By the look of it, snow is still underestimated, but a little bit less so.

  var cover = collection
    .filterDate(date, date.advance(1, 'years'))
    .map(function(image) {
      var mask = image.select('NDSI_Snow_Cover_Basic_QA').not() // Only take best quality pixels
      return image.select('NDSI_Snow_Cover').gte(30)
        .updateMask(mask)
    })
    .reduce(ee.Reducer.sum().combine(ee.Reducer.count(), null, true))
  var fractionOfDaysWithSnow = cover.select('NDSI_Snow_Cover_sum')
    .divide(cover.select('NDSI_Snow_Cover_count'))
  var daysInYear = date.advance(1, 'years').difference(date, 'days')
  return fractionOfDaysWithSnow.multiply(daysInYear)
    .round().int16()
    .set('date', date.format('yyyy-MM-dd'))
    .clip(aoi)

To get the average for the whole region, use ee.Image.reduceRegion().

var aggregateDays = daysSnowCollection.map(function (image) {
  var mean = image.reduceRegion({
    reducer: ee.Reducer.mean(), 
    geometry: aoi, 
    scale: 500,
    maxPixels: 1e13
  }).values().get(0)
  return ee.Feature(null, {mean: mean})
})

https://code.earthengine.google.com/6c71a6ecbc062e4c96398d82b8d3c00a