0

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'))
}

1 Answer 1

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

3
  • Thanks for your reply! It sounds perfect and indeed, less underestimated. However, I'm not sure I understand how all the code works. With this method, we always take one image per day, right? Then, for each of these images, we keep only the pixels of good quality and with an NDSI > 30. This is the 'fractionOfDaysWithSnow' variable which I am not sure I understand what it does. Since just before, we sum it up, it means that I should have, for each pixel, the number of days in the year when there was snow, right? So I don't understand why multiply by the number of days.
    – Elodie
    Feb 27, 2023 at 9:00
  • Masked pixels will not be included in the count or the sum, but I'm not masking non-snowy pixels. So you have a images with 0 for no snow, 1 for snow. The count gives the total number of valid observations, the sum gives the number of snowy observations. The code doesn't do a sum then a count, it's doing both in parallel. Dividing the count with the sum gives you the fraction of snowy days. Assuming the fraction is the same for any masked pixels, you get the number of snowy days by multiplying with the number of days in the year. Feb 28, 2023 at 6:45
  • Great, thank you very much for this detailed explanation and your help! I understand better :) and it helps me a lot.
    – Elodie
    Mar 1, 2023 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.