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.