Cloudless NDVI mosaic (-1 to 1 values)

Hi,

It is possible to create a cloudless mosaic of monthly average NDVI using the evalscript. First you need to use preProcessScenes to sort the acquisitions (or select leastRecent to have all acquisitions in the chronological order) and updateOutput to reset the number of output bands (which is the same as the number of months) within your time range.

Next, in the evaluatePixel you can compute the NDVI for all acquisitions and then average them by month. Lastly, you can output the monthly average NDVI into separate bands in the output TIFF. Below is an example script that creates a cloudless mosaic of monthly average NDVI.

//VERSION=3

function setup() {
  return {
    input: ["B04", "B08", "CLM", "dataMask"],
    output: { bands: 1, sampleType: "FLOAT32" },
    mosaicking: "ORBIT"
  };
}

// Sort acquisitions
function preProcessScenes (collections) {
  collections.scenes.orbits = collections.scenes.orbits.sort(
    (s1, s2) => new Date(s1.dateFrom) - new Date(s2.dateFrom)
  );
  return collections
}

// define the number of bands needed for the time range
function updateOutput(outputs, collection) {
  const scenes = collection.scenes
  if (scenes.length === 0) {
    n_bands = 1;
  } else {
    const first_year = scenes[0].date.getFullYear()
    const first_month = scenes[0].date.getMonth() + 1
    const last_year = scenes[scenes.length - 1].date.getFullYear()
    const last_month = scenes[scenes.length - 1].date.getMonth() + 1
    if (first_year === last_year) {
      n_bands = last_month - first_month + 1;
    }
    if (first_year < last_year) {
      n_bands = 
        (12 - first_month + 1) + 
        (last_year - (first_year + 1)) * 12 +
        (last_month)
    }
  }
  Object.values(outputs).forEach((output) => {
    output.bands = n_bands;
  });
}


function evaluatePixel(samples, scenes) {
  // create a map for month and the indeices of acquisitions in the month
  let month_to_idx = {}
  scenes.forEach((scene, idx) => {
    year = scene.date.getFullYear()
    month = scene.date.getMonth() + 1
    key = [year, month].join("-")
    if (key in month_to_idx) {
      month_to_idx[key].push(idx)
    } else {
      month_to_idx[key] = [idx]
    }
  })

  let n_months = Object.keys(month_to_idx).length;
  let avg_ndvi = new Array(n_months).fill(NaN);

  // loop through months in the map and compute compute the average of ndvi
  Object.entries(month_to_idx).forEach(([, idxs], month) => {
    ndvi_monthly_count = 0
    ndvi_monthly_avg = 0;
    monthly_ndvi = []
    idxs.forEach((idx) => {
      if (samples[idx].CLM === 0 && samples[idx].dataMask === 1) {
        ndvi = (samples[idx].B08 - samples[idx].B04) / (samples[idx].B08 + samples[idx].B04);
        if (isFinite(ndvi)) {
          ndvi_monthly_count ++;
          monthly_ndvi.push(ndvi);
          ndvi_monthly_avg = (ndvi_monthly_avg * (ndvi_monthly_count - 1) + ndvi) / ndvi_monthly_count;
          avg_ndvi[month] = ndvi_monthly_avg;
        }

      }
    });
  });
  return avg_ndvi;
}

Thank you for your reply chung! I’m just slightly confused. Is this all done in the EO browser GUI? or does this need to be done through the API? which part is done where?

I get this error when I put it into eo browser : ‘Dataset with id: 0 not found’

The example multi-temporal script is mainly written for requesting NDVI monthly average within a specific time range using Process API or Batch Processing API. The output will be a multi-band TIFF and each band represents the average NDVI of a month, e.g., band 1 is the average NDVI in November and band 2 is the average NDVI in December if the time range is set from 2022-11-01 to 2022-12-31.

You could try the following request with Requests Builder. Simply copy the following code and paste it to the Request preview window then parse it. You should get a 2-band TIFF which has the average NDVI in November and December of the AOI.

curl -X POST https://services.sentinel-hub.com/api/v1/process \
 -H 'Content-Type: application/json' \
 -H 'Authorization: Bearer ' \
 -d '{
  "input": {
    "bounds": {
      "bbox": [
        4524858.253997,
        2088521.994524,
        4532692.194797,
        2093982.407454
      ],
      "properties": {
        "crs": "http://www.opengis.net/def/crs/EPSG/0/3035"
      }
    },
    "data": [
      {
        "dataFilter": {
          "timeRange": {
            "from": "2022-11-01T00:00:00Z",
            "to": "2022-12-31T23:59:59Z"
          }
        },
        "type": "sentinel-2-l2a"
      }
    ]
  },
  "output": {
    "resx": 10,
    "resy": 10,
    "responses": [
      {
        "identifier": "default",
        "format": {
          "type": "image/tiff"
        }
      }
    ]
  },
  "evalscript": "//VERSION=3\n\nfunction setup() {\n  return {\n    input: [\"B04\", \"B08\", \"CLM\", \"dataMask\"],\n    output: { bands: 1, sampleType: \"FLOAT32\" },\n    mosaicking: \"ORBIT\"\n  };\n}\n\n// Sort acquisitions\nfunction preProcessScenes (collections) {\n  collections.scenes.orbits = collections.scenes.orbits.sort(\n    (s1, s2) => new Date(s1.dateFrom) - new Date(s2.dateFrom)\n  );\n  return collections\n}\n\n// define the number of bands needed for the time range\nfunction updateOutput(outputs, collection) {\n  const scenes = collection.scenes\n  if (scenes.length === 0) {\n    n_bands = 1;\n  } else {\n    const first_year = scenes[0].date.getFullYear()\n    const first_month = scenes[0].date.getMonth() + 1\n    const last_year = scenes[scenes.length - 1].date.getFullYear()\n    const last_month = scenes[scenes.length - 1].date.getMonth() + 1\n    if (first_year === last_year) {\n      n_bands = last_month - first_month + 1;\n    }\n    if (first_year < last_year) {\n      n_bands = \n        (12 - first_month + 1) + \n        (last_year - (first_year + 1)) * 12 +\n        (last_month)\n    }\n  }\n  Object.values(outputs).forEach((output) => {\n    output.bands = n_bands;\n  });\n}\n\n\nfunction evaluatePixel(samples, scenes) {\n  // create a map for month and the indeices of acquisitions in the month\n  let month_to_idx = {}\n  scenes.forEach((scene, idx) => {\n    year = scene.date.getFullYear()\n    month = scene.date.getMonth() + 1\n    key = [year, month].join(\"-\")\n    if (key in month_to_idx) {\n      month_to_idx[key].push(idx)\n    } else {\n      month_to_idx[key] = [idx]\n    }\n  })\n\n  let n_months = Object.keys(month_to_idx).length;\n  let avg_ndvi = new Array(n_months).fill(NaN);\n\n  // loop through months in the map and compute compute the average of ndvi\n  Object.entries(month_to_idx).forEach(([, idxs], month) => {\n    ndvi_monthly_count = 0\n    ndvi_monthly_avg = 0;\n    monthly_ndvi = []\n    idxs.forEach((idx) => {\n      if (samples[idx].CLM === 0 && samples[idx].dataMask === 1) {\n        ndvi = (samples[idx].B08 - samples[idx].B04) / (samples[idx].B08 + samples[idx].B04);\n        if (isFinite(ndvi)) {\n          ndvi_monthly_count ++;\n          monthly_ndvi.push(ndvi);\n          ndvi_monthly_avg = (ndvi_monthly_avg * (ndvi_monthly_count - 1) + ndvi) / ndvi_monthly_count;\n          avg_ndvi[month] = ndvi_monthly_avg;\n        }\n\n      }\n    });\n  });\n  return avg_ndvi;\n}"
}'

The purpose of EO Browser is more focused on visualisation, so the script needs to be modified a bit. Note that with EO Browser this script only works when the time range is set to one month, e.g., from 2022-12-01 to 2022-12-31. Please take a look at the visualisation of NDVI monthly average in December 2022.

This is brilliant chung thank you! i’ll give try this