Adapting a custom script from process api to statistical api

EO Products / Scripts
1

Hi,

I’m having trouble in using one of the scripts on the custom-scripts repository, when I try to adapt it to the statistical api call. I’m trying it out in the requests builder, and it works well while using the process api. I just can’ really pinpoint how exactly I should adapt the script in order to fit the statistical api request. The script in use is [Vegetation Productivity Indicator | Sentinel-Hub custom scripts](https://Vegetation Productivity Indicator)

My main doubt is on how to edit the evaluatePixel return, and also the time range, which in the Process API is part of the date filter, but in statistical api is part of aggregation (together with aggregation interval).

Thank you in advance!

2

StatisticalAPI requires dataMask output to be added, so change the setup part of the script to this:

  return {
    input: ["B04", "B08"],
    output: [
      { id: "vpi", bands: 1},
      { id: "dataMask", bands: 1}
    ],
    mosaicking: "ORBIT"
  };

and then change the return inside the evaluatePixel function to return an object with two properties instead of an array, like this:

  return {
    dataMask: [hist_ndvi.length > 0 ? 1 : 0],
    vpi: [percentileOfScore(hist_ndvi, observed) / 100]
  }

As for the time intervals - statistical API will always run evaluatePixel with data from distinct (non-overlapping) intervals, and for this indicator you need long-term history, so the best you can do is compute the indicator for a single long interval in one request. To do this, choose the start and end dates, and set a really long aggregation interval e.g. 100 years (P100Y).

An example CURL request (that you can copy to Request Builder and click Parse): 
curl -X POST https://services.sentinel-hub.com/api/v1/statistics 
 -H 'Content-Type: application/json'
 -d '{
  "input": {
    "bounds": {
      "bbox": [
        12.44,
        41.87,
        12.46,
        41.89
      ]
    },
    "data": [
      {
        "dataFilter": {},
        "type": "sentinel-2-l1c"
      }
    ]
  },
  "aggregation": {
    "timeRange": {
      "from": "2016-03-01T00:00:00Z",
      "to": "2022-01-06T23:59:59Z"
    },
    "aggregationInterval": {
      "of": "P100Y",
      "lastIntervalBehavior": "SHORTEN"
    },
    "width": "20",
    "height": "20",
    "evalscript": "      //VERSION=3\nfunction setup() {\n  return {\n    input: [\"B04\", \"B08\"],\n    output: [\n      {id: \"vpi\", bands: 1},\n      {id: \"dataMask\", bands: 1}\n    ],\n    mosaicking: \"ORBIT\"\n  };\n}\n\nconst msInDay = 24 * 60 * 60 * 1000;\nconst msInYear = 365.25 * msInDay;\nconst toleranceDays = 10;\nconst toleranceMs = toleranceDays * msInDay;\n\nvar metadata = undefined;\n\nfunction filterScenes(scenes, inputMetadata) {\n  scenes = scenes.sort((s1, s2) => s2.date - s1.date);\n  const observed = scenes[0].date;\n  var newScenes = [scenes[0]];\n  for (var historical = observed - msInYear; historical >= inputMetadata.from - toleranceMs; historical -= msInYear) {\n    newScenes.push(findClosest(scenes, historical));\n  }\n  newScenes = newScenes.filter(scene => scene != null);\n  metadata = {\n    observed: observed.toISOString(),\n    historical: newScenes.slice(1).map(scene => scene.date.toISOString())\n  }\n  return newScenes;\n}\n\nfunction findClosest(scenes, date) {\n  var closestDt = toleranceMs + 1, closestScene = null;\n  for (var i = 0; i < scenes.length; i++) {\n    const dt = Math.abs(scenes[i].date - date);\n    if (dt < closestDt) {\n      closestDt = dt;\n      closestScene = scenes[i];\n    }\n  }\n  return closestScene;\n}\n\nfunction percentileOfScore (data, value)  {\n  // Calculate the percentile rank of a value relative to a list of values.\n\n    if (!data.length) {return [0];}\n\n    data.sort();\n\n    let lowerCount = 0;\n    let sameCount = 0;\n\n    for (let i = 0; i < data.length; i++) {\n      if (data[i] < value) {\n        lowerCount++;\n      } else if (data[i] === value) {\n        sameCount++;\n      } else {\n        break;\n      }\n    }\n\n    return (lowerCount + 0.5 * sameCount) / data.length * 100;\n}\n\nfunction updateOutputMetadata(scenes, inputMetadata, outputMetadata) {\n  outputMetadata.userData = metadata;\n}\n\nfunction evaluatePixel(samples, scenes) {\n  const observed = index(samples[0].B08, samples[0].B04);\n  var hist_ndvi = [];\n  for (var i = 1; i < samples.length; i++) {\n    const ndvi = index(samples[i].B08, samples[i].B04);\n    hist_ndvi.push(ndvi);\n  }\n\n  return {\n    dataMask: [hist_ndvi.length > 0 ? 1 : 0],\n    vpi: [percentileOfScore(hist_ndvi, observed) / 100]\n  }\n}"
  },
  "calculations": {
    "default": {}
  }
}'
3

Thanks a lot for the reply, that worked! Can the statsitics api also take in data fusion? I tried to adapt the script s1+s2 for calculating ndvi, following the same logic as above, but only got NaN results or errors :

curl -X POST https://services.sentinel-hub.com/api/v1/statistics 
-H 'Content-Type: application/json' 
-d '{
 "input": {
   "bounds": {
     "bbox": [
       0.576168,
       5.889986,
       0.592098,
       5.900488
     ]
   },
   "data": [
     {
       "dataFilter": {},
       "processing": {
         "orthorectify": "true",
         "backCoeff": "SIGMA0_ELLIPSOID"
       },
       "id": "s1",
       "type": "sentinel-1-grd"
     },
     {
       "dataFilter": {},
       "id": "l2a",
       "type": "sentinel-2-l2a"
     }
   ]
 },
 "aggregation": {
   "timeRange": {
     "from": "2019-04-05T00:00:00Z",
     "to": "2019-04-29T23:59:59Z"
   },
   "aggregationInterval": {
     "of": "P5D",
     "lastIntervalBehavior": "SHORTEN"
   },
   "resx": "10",
   "resy": "10",
   "evalscript": "//VERSION=3\nfunction setup() {\n  return {\n    input: [{\n        datasource: \"s1\",\n        bands: [\"VV\", \"VH\", \"dataMask\"]\n      },\n      {\n        datasource: \"l2a\",\n        bands: [\"B02\", \"B03\", \"B08\", \"B04\", \"SCL\", \"dataMask\"]\n      }\n    ],\n    output: [\n      {\n        id: \"data\",\n        bands: 1\n    },\n      {\n        id: \"dataMask\",\n        bands: 1\n      }\n    ]\n  }\n}\n\nfunction toDb(linear) {\n  // Convert the linear backscatter to DB (Filgueiras et al. (2019), eq. 3)\n  return 10 * Math.LN10 * linear\n}\n\nfunction calc_s1_ndvi(sigmaVV, sigmaVH) {\n  // Convert sigma0 to Decibels\n  let vh_Db = toDb(sigmaVH)\n  let vv_Db = toDb(sigmaVV)\n  // Calculate NRPB (Filgueiras et al. (2019), eq. 4)\n  let NRPB = (vh_Db - vv_Db) / (vh_Db + vv_Db)\n  // Calculate NDVI_nc with approach A3 (Filgueiras et al. (2019), eq. 14)\n  let NDVInc = 2.572 - 0.05047 * vh_Db + 0.176 * vv_Db + 3.422 * NRPB\n  return NDVInc\n}\n\nfunction evaluatePixel(samples) {\n  var s1 = samples.s1[0]\n  var s2 = samples.l2a[0]\n\n  // Calculate S2 NDVI\n    let ndvi = (samples.B08 - samples.B04)/(samples.B08 + samples.B04)\n    \n  // Calculate S1 NDVI\n  let s1_ndvi = calc_s1_ndvi(s1.VV, s1.VH)\n  // Use the S2-L2A classification to identify clouds\n  if ([7, 8, 9, 10].includes(s2.SCL)) {\n    // If clouds are present use S1 NDVI\n    return {\n      data: [s1_ndvi],\n      dataMask: [samples.dataMask]\n    }\n  } else {\n    // Otherwise use s2 NDVI\n    return {\n        data: [ndvi],\n        // Exclude nodata pixels, pixels where ndvi is not defined and water pixels from statistics:\n        dataMask: [samples.dataMask]\n    }\n  }\n}"
 },
 "calculations": {
   "default": {}
 }
}'

Any suggestions on if this might work?

Thanks again!

4

Hi @miha do you think what I wrote is possible? thanks!

5

It should work. I don’t know the details of what you want to achieve, but I noticed that the line let ndvi = (samples.B08 - samples.B04)/(samples.B08 + samples.B04) looks like it should use s2 instead of samples.

6

I had another look at your request, there were a couple of issues:

  • samples.dataMask is undefined in fusion requests; you can use s1.dataMask or s2.dataMask instead

  • the request uses EPSG:4326 and resolution of 10 degrees (resX and resY are always in CRS units, which in case of geographic CRS means degrees). This seems to have caused some errors on the server. I changed the aggregation to use width: 256, height: 196.66 - this not only prevents the error, but also calculates statistics on ~45k pixels within the area-of-interest; with too coarse resolution (10 degrees) you would only evaluate one pixel per interval

curl -X POST https://services.sentinel-hub.com/api/v1/statistics 
 -H 'Content-Type: application/json'
 -d '{
  "input": {
    "bounds": {
      "bbox": [
        0.576168,
        5.889986,
        0.592098,
        5.900488
      ]
    },
    "data": [
      {
        "dataFilter": {},
        "processing": {
          "orthorectify": "true",
          "backCoeff": "SIGMA0_ELLIPSOID"
        },
        "id": "s1",
        "type": "sentinel-1-grd"
      },
      {
        "dataFilter": {},
        "id": "l2a",
        "type": "sentinel-2-l2a"
      }
    ]
  },
  "aggregation": {
    "timeRange": {
      "from": "2019-04-05T00:00:00Z",
      "to": "2019-04-29T23:59:59Z"
    },
    "aggregationInterval": {
      "of": "P5D",
      "lastIntervalBehavior": "SHORTEN"
    },
    "width": 256,
    "height": 169.666,
    "evalscript": "//VERSION=3\nfunction setup() {\n  return {\n    input: [{\n        datasource: \"s1\",\n        bands: [\"VV\", \"VH\", \"dataMask\"]\n      },\n      {\n        datasource: \"l2a\",\n        bands: [\"B02\", \"B03\", \"B08\", \"B04\", \"SCL\", \"dataMask\"]\n      }\n    ],\n    output: [\n      {\n        id: \"data\",\n        bands: 1\n    },\n      {\n        id: \"dataMask\",\n        bands: 1\n      }\n    ]\n  };\n}\n\nfunction toDb(linear) {\n  // Convert the linear backscatter to DB (Filgueiras et al. (2019), eq. 3)\n  return 10 * Math.LN10 * linear;\n}\n\nfunction calc_s1_ndvi(sigmaVV, sigmaVH) {\n  // Convert sigma0 to Decibels\n  const vh_Db = toDb(sigmaVH);\n  const vv_Db = toDb(sigmaVV);\n  // Calculate NRPB (Filgueiras et al. (2019), eq. 4)\n  const NRPB = (vh_Db - vv_Db) / (vh_Db + vv_Db);\n  // Calculate NDVI_nc with approach A3 (Filgueiras et al. (2019), eq. 14)\n  const NDVInc = 2.572 - 0.05047 * vh_Db + 0.176 * vv_Db + 3.422 * NRPB;\n  return NDVInc;\n}\n\nfunction evaluatePixel(samples) {\n  const s1 = samples.s1[0];\n  const s2 = samples.l2a[0];\n\n  // Use the S2-L2A classification to identify clouds\n  if ([7, 8, 9, 10].includes(s2.SCL)) {\n    // If clouds are present use S1 NDVI\n    const s1_ndvi = calc_s1_ndvi(s1.VV, s1.VH);\n    return {\n      data: [s1_ndvi],\n      dataMask: [s1.dataMask]\n    };\n  } else {\n    const ndvi = (s2.B08 - s2.B04)/(s2.B08 + s2.B04);\n    // Otherwise use s2 NDVI\n    return {\n        data: [ndvi],\n        // Exclude nodata pixels, pixels where ndvi is not defined and water pixels from statistics:\n        dataMask: [s2.dataMask]\n    };\n  }\n}"
  },
  "calculations": {
    "default": {}
  }
}'

This request now does the following:

  • for each 5-days interval (there are 5 such intervals in the time range), create an S1 mosaic and L2A mosaic (using default mosaicking order, which is mostRecent) for the requested area-of-interest and downsampled to match the requested width & height
  • for every pixel in the mosaic, calculate result using evaluatePixel (if the latest L2A image is cloudy, use S1 for NDVI approximation; else use L2A NDVI)
  • for all per-pixel results in the interval (256 * 169 = 43264 for each interval), aggregate statistics