Skip to main content

I used the following modified evalscript to download statistics for NDVI and NDMI for an AOI. here is the evalscript:


ndvi_new_evalscript = """

//VERSION=3

function setup() {
return {
input: [{
bands: ["B04","B05","B08", "B8A","B11", "B12", "CLM", "CLP", "dataMask"],
units: "DN"
}],
output: [
{
id: "bands",
bands: ["B04"],
sampleType: "UINT16"
},
{
id: "masks",
bands: ["CLM"],
sampleType: "UINT16"
},
{
id: "indices",
bands: ["NDVI", "NDVI_RE1", "NDMI","CLP"],
sampleType: "UINT16"
},
{
id: "dataMask",
bands: 1
}]
}
}

function evaluatePixel(samples) {
// Normalised Difference Vegetation Index and variation
let NDVI = (samples.B08 - samples.B04) / (samples.B08 + samples.B04);
let NDVI_RE1 = index(samples.B08, samples.B05);
let NDMI = (samples.B08 - samples.B11) / (samples.B08 + samples.B11);

// Bare Soil Index
// let NBSI = index((samples.B11 + samples.B04), (samples.B08 + samples.B02));

// cloud probability normalized to interval [0, 1]
let CLP = samples.CLP / 255.0;

// masking cloudy pixels
let combinedMask = samples.dataMask
if (samples.CLM > 0) {
combinedMask = 0;
}

// masking water pixels
let waterMask = samples.dataMask
if (samples.SCL == 6) {
waterMask = 0;
}



const f = 5000;
return {
bands: [samples.B04],
masks: [samples.CLM],
indices: [toUINT(NDVI, f), toUINT(NDVI_RE1, f), toUINT(CLP, f), toUINT(NDMI, f) ],
dataMask: [combinedMask*waterMask]
};
}

function toUINT(product, constant){
// Clamp the output to [-1, 10] and convert it to a UNIT16
// value that can be converted back to float later.
if (product < -1) {
product = -1;
} else if (product > 10) {
product = 10;
}
return Math.round(product * constant) + constant;
}
"""

The python call was like this:


aggregation = SentinelHubStatistical.aggregation(
evalscript=ndvi_new_evalscript, time_interval=yearly_time_interval, aggregation_interval="P1D", resolution=(10, 10)
)

#calculations = {"default": {"statistics": {"default": {"percentiles": {"k": [5, 50, 95]}}}}}
histogram_calculations = {"ndvi": {"histograms": {"default": {"nBins": 20, "lowEdge": -1.0, "highEdge": 1.0}}}}

features_requests = []
for geo_shape in AOI_polygon.geometry.values:
request = SentinelHubStatistical(
aggregation=aggregation,
input_data=[SentinelHubStatistical.input_data(DataCollection.SENTINEL2_L2A)],
geometry=Geometry(geo_shape, crs=CRS(AOI_polygon.crs)),
calculations=histogram_calculations,
config=config
)

features_requests.append(request)

I got an negative NDVI value and I wanted to check if it is the same Sentinel-2 image that is returned with just getting the least cloud cover image from that particular date or interval (2018-6-29 to 2018-6-30). Looks like the image from 2018-6-29 to 2018-6-30 returns a median 0.3 NDMI value as opposed to -1 that was returned by the statistical API.


It can also be that I am not using the exact evalscript and the masks that were used in the statistical API. Do you know if I can do a one-to-one comparison as further checks.

Hi,

To use the image acquisition with the least cloud cover in your time interval you will need to explicitly define this as an argument in your Statistical API request; this will be part of the input_data section of your request payload. By default the API will obtain the most recent image, so this could be the reason why you are returning different NDMI values.

input_data=[
SentinelHubStatistical.input_data(
DataCollection.SENTINEL2_L2A,
other_args={"dataFilter": {"mosaickingOrder": "leastCC"}},
),

If I have misinterpreted your question, then let us know 👍


Thanks for the answers. Is there a way that I can also return RGB images for which the mean NDVI was calculated from the statistical API


Hi, you can do this using a separate Process API request. Statistical API doesn’t return any pixels, only the underlying statistics.


Yes, I am wondering how I can adapt the same script to produce images since it has data masks and cloud masks.


As I said, you need to use Process API to produce images; your current script uses Statistical API so you cannot produce images.


If you are looking for a True Color evalscript that uses a cloud mask, then I can point you to this example.


//VERSION=3
function setup() {
return {
input: ["B02", "B03", "B04", "CLM"],
output: { bands: 3 }
}
}

function evaluatePixel(sample) {
if (sample.CLM == 1) {
return [0.75 + sample.B04, sample.B03, sample.B02]
}
return [3.5*sample.B04, 3.5*sample.B03, 3.5*sample.B02];
}

Reply