Monthly Snow Report Script 

//VERSION=3
/*
Source: @nkarasiak / www.karasiak.net

Monthly snow report.
Tired of waiting the perfect image with no cloud to show the snow cover ? This monthly snow report is here for you.
This script will find where the snow is persistent within the last 31 days (from chosen date).
In order to well represent the land-cover, the script will store each pertinent date in a list and will represents the median value.

This sentinel-2 script voluntarily color in green the image to better differentiate the snow from the other land-cover.

Script requires multi-temporal processing.
*/

// Put 3 to have synthesis of the last 3 months
var numberOfMonthsToUse = 1;
// Thresold to consider pixel as snow
var NDSIthresold = 0.20;
// In order to dismiss snow from water
var redThresold = 0.2;
// In order to dismiss clouds
var blueThresold = 0.18;

function setup() {
  return {
    input: [{
      bands: [
          "B02",
          "B03",
          "B04",
          "B11"
      ]
    }],
    output: { bands: 3 },
    mosaicking: "ORBIT"
  }
}


function NDSI(sample) {
    return ((sample.B03 - sample.B11) / (0.01 + sample.B03 + sample.B11));
}

function median(values) {
    // from https://stackoverflow.com/questions/45309447/calculating-median-javascript
    if (values.length === 0) return 0;

    values.sort(function(a, b) {
        return a - b;
    });

    var half = Math.floor(values.length / 2);

    if (values.length % 2)
        return values[half];

    return (values[half - 1] + values[half]) / 2.0;
}

function evaluatePixel(samples, scenes) {
    // for snow scene
    let snowyCount = 0;
    let snowB02 = [];
    let snowB03 = [];
    let snowB04 = [];

    // for unsnow scene
    let B02 = [];
    let B03 = [];
    let B04 = [];

    // to manage image length between tiles
    realSampleLength = 0;

    for (i = 0; i < samples.length; i++) {
        // in order to avoid black pixel (the ones between tiles)

        if ((samples[i].B02 > 0) || (samples[i].B03 > 0)) {
            realSampleLength++;

            // found snow
            if ((NDSI(samples[i]) > NDSIthresold) & (samples[i].B04 > redThresold)) {
                snowyCount++;
                snowB02.push(samples[i].B02);snowB03.push(samples[i].B03);snowB04.push(samples[i].B04);
            }
            // maybe it is snowy but clouds...
            // so count as snow but not memorizing pixel values.
            else if ((samples[i].B02 > blueThresold) & (samples[i].B04 > redThresold)) {
                snowyCount++;
            }
            // not snow
            else if (samples[i].B02 < blueThresold) {
                B02.push(samples[i].B02);B03.push(samples[i].B03);B04.push(samples[i].B04);
            }
        }
    }

    if ((snowB02 == 0) & (B02.length == 0)) {
        // no valid data available : dark green color
        colorMap = [.05,.2,.05];
    } else if ((snowyCount == realSampleLength) || (B02.length<1)) {
        // snowColorMap
        colorMap = [1.1 * median(snowB04), 1.3 * median(snowB03), 1.1 * median(snowB02)];
    } else if (B02.length >= 1) {
        // defaultColorMap
        colorMap = [1.5 * median(B04), 2.5 * median(B03), 1.5 * median(B02)];
    }

    return colorMap;
}

function preProcessScenes (collections) {
    var scenes = collections.scenes.orbits;
    scenes = scenes.sort((s1, s2) => new Date(s2.dateFrom) - new Date(s1.dateFrom));
    var latest_scene = new Date(scenes[0].dateFrom)
    collections.scenes.orbits = collections.scenes.orbits.filter(function (scene) {
        return new Date(scene.dateFrom) >= latest_scene - numberOfMonthsToUse * 31 * 24 * 3600 * 1000
    })
    return collections
}

Evaluate and Visualize

General description of the script

Tired of waiting the perfect image with no cloud to show the snow cover? This monthly snow report script is for you.

This code will find where the snow is persistent for the last 30 days (from the chosen date). In order to well represent the land-cover, the script will store each pertinent date in a list and will represent the median value.

As the aim of the script is to represent the snow cover, all the other land-cover are saturated in green in order to easily see the snow.

The limitations are essentially:

  • a few white rooftops
  • when no uncloudy pixel is available in the previous 30 days

Thus, it is possible to make the script working with 90 days in order to have a trimonthly synthesis by changing settings to: “numberOfMonthsToUse = 3”.

Author of the script

Nicolas Karasiak

Description of representative images

Persistent snow cover in February 2019 in Corsica

Persistent snow cover in February 2019 in Corsica

Persistent snow cover in Pyrénées in March 2019

Pyrénées persistent snow cover in March 2019

Mer de Glace / Alpes persistent snow cover in the previous 30 days of 24 april 2019.

Chamonix persistent snow cover in the previous 30 days of 24 April 2019.

See the supplementary material for more examples.

Credits