Sentinel hub Statistical API- costum calculations

Apps / Services Python SDK
1

Hello,

I am reading the statistical API and I was wondering if there is any way to add my own calcualtions, for example. if I want to calculate the precentage of pixels in polygon that their NDVI value is higher than 0.8, or between two values…
Is there a way to do this with the statistical api?

Best
Reut

2

Hi Reut,

this is possible. If you calculate in your custom script a “mask” for NDVI larger than 0.8,
like

let NDVI_ABOVE_08 = (NDVI > 0.80) ? 1 : 0;

and add it to the output bands, then StatAPI will aggregate values of this band over all pixels. mean value of NDVI_ABOVE_08 in this case will represent the fraction of pixels within the request geometry with NDVI value above 0.8.

Best,
Anze

1 Like
3

thank you for the fast respond! if I want to get precentage between two values, how will that be written?

4

Two options:

  • add output for larger than and lower than
let NDVI_UNDER_04 = (NDVI < 0.40) ? 1 : 0;
let NDVI_UNDER_08 = (NDVI < 0.80) ? 1 : 0;
  • or join the conditions in a single output
let NDVI_IN_BETWEEN = (NDVI > 0.40 && NDVI < 0.8) ? 1 : 0;

In the first case you would then subtract fractions of pixels under 0.8 and 0.4.

1 Like
5

I have tried to add this to my evalscript but get errors:

ndvi_evalscript= """
//VERSION=3

function setup() {
  return {
    input: [
      {
        bands: [
          "B04",
          "B08",
          "dataMask"
        ]
      }
    ],
    output: [
      {
        id: "ndvi",
        bands: 1
      },
      {
        id: "dataMask",
        bands: 1
      }
    ]
  }
}

function evaluatePixel(samples) {
    return {
      let NDVI=index(samples.B08, samples.B04)
      let NDVI_UNDER_04 = (NDVI < 0.40) ? 1 : 0;
      ndvi: [NDVI,NDVI_UNDER_04],
      dataMask: [samples.dataMask]
    };
}
"""


#run statistics like in the documentation
yearly_time_interval = '2020-01-01', '2020-12-31'


aggregation = SentinelHubStatistical.aggregation(
    evalscript=ndvi_evalscript,
    time_interval=yearly_time_interval,
    aggregation_interval='P1D',
    resolution=(10, 10)
)

input_data = SentinelHubStatistical.input_data(
    DataCollection.SENTINEL2_L2A
)

histogram_calculations = {
    "ndvi": {
        "histograms": {
            "default": {
                "nBins": 20,
                "lowEdge": -1.0,
                "highEdge": 1.0
            }
        }
    }
}

ndvi_requests = []

for geo_shape in polygons_gdf.geometry.values:
    request = SentinelHubStatistical(
        aggregation=aggregation,
        input_data=[input_data],
        geometry=Geometry(geo_shape, crs=CRS(polygons_gdf.crs)),
        calculations=histogram_calculations,
        config=config
    )
    ndvi_requests.append(request)

%%time

download_requests = [ndvi_request.download_list[0] for ndvi_request in ndvi_requests]

client = SentinelHubStatisticalDownloadClient(config=config)

ndvi_stats = client.download(download_requests)

len(ndvi_stats)

That part runs but when I print the ndvi_stats can see error of bad request/execution error:

[{'data': [{'interval': {'from': '2020-01-04T00:00:00Z',
     'to': '2020-01-05T00:00:00Z'},
    'error': {'type': 'EXECUTION_ERROR'}},
   {'interval': {'from': '2020-01-09T00:00:00Z', 'to': '2020-01-10T00:00:00Z'},
    'error': {'type': 'EXECUTION_ERROR'}},
   {'interval': {'from': '2020-01-14T00:00:00Z', 'to': '2020-01-15T00:00:00Z'},
    'error': {'type': 'EXECUTION_ERROR'}},
   {'interval': {'from': '2020-01-19T00:00:00Z', 'to': '2020-01-20T00:00:00Z'},
    'error': {'type': 'EXECUTION_ERROR'}},
   {'interval': {'from': '2020-01-24T00:00:00Z', 'to': '2020-01-25T00:00:00Z'},
    'error': {'type': 'BAD_REQUEST'}},
.....

When I run this withou the NDVI_UNDER etc. it works just like in the tutorial.

Do you see my mistake?
thanks a lot :slight_smile:

6

The mistake is in the setup part: you defined “ndvi” output to contain single band, but your evalscript is returning 2.

Additionally, the let parts of the evaluatePixel should be outside the return statement.

1 Like
7

Thank you, that worked but it seems like the query returns 1 or 0 but not the percentage itself

8

How large is your area? If there is only a single Sentinel-2 pixel, that would be expected.

9

it is more than one pixel for sure

10

maybe is something about the data type?

11

One thing to check is the geometry and it’s CRS. If you are making requests for geometries in WGS84, then the resolution you specify in aggregation is in degrees! Your geometries are probably smaller than 10deg x 10deg, meaning the result is for 1 pixel.

1 Like
12

Thank you Batic, that was the issue :slight_smile:

1 Like