Quarterly median or mean composites of Sentinel 1 data

I’m trying to get quarterly (3-month window) median composites of the VV/VH sentinel 1 imagery.

I’ve set up the following script, but it doesn’t seem to work (it returns all the images rather than the median, and the entire array is zeros).

The script has mean/median functions, and I’m testing them both out to see which is a better compositing method.

Curious if anyone has input into either/and:

  • How to set up a median composite script
  • How to run the script for each quarter within a year (currently I’m planning to just call the function 4 times with different start and end dates
//VERSION=3 (auto-converted from 1) + dataMask to control background color + scaling to UINT16 range

function mean(values) {
    var total = 0, i;
    for (i = 0; i < values; i += 1) {
        total += values[i];
    }
    return total / values.length;
}

function median(values) {
    if (values.length === 0) return 65535;
    if (values.length === 1) return values[0];

    values.sort(function(a, b) {
        return a - b;
    });

    var half = Math.floor(values.length / 2);
    return values[half];
}

function validate(samples) {
  if (samples.dataMask == 0) {
    return false;
  } else {
    return true;
  }
}


function evaluatePixel(samples) {
    var factor = 65535
    var vv = []
    var vh = []
    var a = 0;

    for (var i = 0; i < samples.length; i++) {
      var sample = samples[i];
      if (sample.VV > 0 && sample.VH > 0) {
        var isValid = validate(sample);

        if (isValid) {
          vv[a] = sample.VV * factor;
          vh[a] = sample.VH * factor;
          a = a + 1;
        }
      }
    }
    return [median(vv), median(vh)];
}

function setup() {
  return {
    input: [{
      bands: [
        "VV",
        "VH",
        "dataMask",
             ]
    }],
    output: {
      bands: 2,
      sampleType:"UINT16",
      mosaicking: "ORBIT"
    }
  }
}

Dear John,

How to set up a median composite script

With a few modifications your script works:

  • You have specified the mosaicking parameter in the wrong place: it should be outside of output in your setup function (see documentation)

  • When using mosaicking: "ORBIT", samples becomes an array that you need to loop through (as mentionned here)

Here is a working version if that helps:

function median(values) {
    if (values.length === 0) return 0;
    if (values.length === 1) return values[0];

    values.sort(function(a, b) {
        return a - b;
    });

    var half = Math.floor(values.length / 2);
    return values[half];
}


function evaluatePixel(samples) {
    // Initialise arrays
    var VV_samples = [];
    var VH_samples = [];
    
    // Loop through orbits and add data
    for (let i=0; i<samples.length; i++){
      // Ignore noData
      if (samples[i].dataMask != 0){
        VV_samples.push(samples[i].VV);
        VH_samples.push(samples[i].VH);
       }
    }  
    
    const factor = 65535;
    
    if (VV_samples.length == 0){
      var VV_median = [factor];
    } else {
      var VV_median = median(VV_samples) * factor;
    }
    if (VH_samples.length == 0){
      var VH_median = [factor];
    } else{
      var VH_median = median(VH_samples) * factor;
    }
    
    return [VV_median, VH_median];
}

function setup() {
  return {
    input: [{
      bands: [
        "VV",
        "VH",
        "dataMask",
             ]
    }],
    output: {
      bands: 2,
      sampleType:"UINT16"
    },
    mosaicking: "ORBIT"
  }
}

For your second question:

How to run the script for each quarter within a year (currently I’m planning to just call the function 4 times with different start and end dates

You can run the script for a whole year and return the quarterly medians by setting up the Evalscript. I have written an example for VV below, you would just need to adapt it to do VH as well. Disclaimer: there are probably more efficient ways to code this, I am just proposing a way that works for me. If any one has improvements, they are most welcome (for example you could optimise it to just have to set the year etc…)!

//VERSION=3
function setup() {
  return {
    input: [{
      bands: [
        "VV",
        "dataMask",
             ]
    }],
    // You can add multiple outputs here (e.g. VV & VH)
    output: {
      id: "VV",
      bands: 4,
      sampleType:"UINT16",
    },
    mosaicking: "ORBIT"
  };
}

function median(values) {
    if (values.length === 0) return 0;
    if (values.length === 1) return values[0];
    values.sort(function(a, b) {
        return a - b;
    });
    var half = Math.floor(values.length / 2);
    return values[half];
}

function check_quarter(indate, q1, q2, q3, q4){
  // Check if date is in one of the intervals
  if (indate >= q1[0] && indate <= q1[1]){
    return 1;
  } else if (indate >= q2[0] && indate <= q2[1]){
    return 2;
  } else if (indate >= q3[0] && indate <= q3[1]){
    return 3;
  } else if (indate >= q4[0] && indate <= q4[1]){
    return 4;
  } else {
    throw new Error("Date out of range!");
  }
}

function evaluatePixel(samples, scenes) {
  // Set date intervals
  var q1 = [new Date("2020-01-01"), new Date("2020-03-31")];
  var q2 = [new Date("2020-04-01"), new Date("2020-06-30")];
  var q3 = [new Date("2020-07-01"), new Date("2020-09-30")];
  var q4 = [new Date("2020-10-01"), new Date("2020-12-31")];
  
  // Initialise object with arrays
  var VV_samples = {q1: [], q2: [], q3:[], q4:[]};

  // Loop through orbits and add data
  for (var i = 0; i < samples.length; i++) {
    // Get quarter
    quarter = check_quarter(new Date(scenes[i].date), q1, q2, q3, q4)
    // Assign to correct array
    if (samples[i].dataMask !=0){
      if (quarter == 1){
        VV_samples.q1.push(samples[i].VV);
      } else if (quarter == 2){
        VV_samples.q2.push(samples[i].VV);
      } else if (quarter == 3){
        VV_samples.q3.push(samples[i].VV);
      } else if (quarter == 4){
        VV_samples.q4.push(samples[i].VV);
      }
    }
  }

  const factor = 65535;

  // Make an object to contain median values
  var VV_median = {};
  
  // Compute median for each quarter
  for (const key in VV_samples) {
    if (VV_samples.hasOwnProperty(key)) {
      if (VV_samples[key].length == 0){
        VV_median[key] = [factor];
      } else {
        VV_median[key] = median(VV_samples[key]) * factor;
      }
    }
  }

  return {
    VV: Object.values(VV_median)
  }
}

Of course in your payload you will have to set the date correctly not to get an error. Below is an example:

curl -X POST https://services.sentinel-hub.com/api/v1/process 
 -H 'Content-Type: application/json' 
 -H 'Authorization: Bearer <my-token>' 
 -d '{
  "input": {
    "bounds": {
      "bbox": [
        12.452079,
        41.860948,
        12.462509,
        41.868378
      ]
    },
    "data": [
      {
        "dataFilter": {
          "timeRange": {
            "from": "2020-01-01T00:00:00Z",
            "to": "2020-12-31T23:59:59Z"
          }
        },
        "type": "sentinel-1-grd"
      }
    ]
  },
  "output": {
    "width": 86.47203647638206,
    "height": 82.71038165936085,
    "responses": [
      {
        "identifier": "VV",
        "format": {
          "type": "image/tiff"
        }
      }
    ]
  },
  "evalscript": "//VERSION=3\nfunction setup() {\n  return {\n    input: [{\n      bands: [\n        \"VV\",\n        \"dataMask\",\n             ]\n    }],\n    // You can add multiple outputs here (e.g. VV & VH)\n    output: {\n      id: \"VV\",\n      bands: 4,\n      sampleType:\"UINT16\",\n    },\n    mosaicking: \"ORBIT\"\n  };\n}\n\nfunction median(values) {\n    if (values.length === 0) return 0;\n    if (values.length === 1) return values[0];\n    values.sort(function(a, b) {\n        return a - b;\n    });\n    var half = Math.floor(values.length / 2);\n    return values[half];\n}\n\nfunction check_quarter(indate, q1, q2, q3, q4){\n  // Check if date is in one of the intervals\n  if (indate >= q1[0] && indate <= q1[1]){\n    return 1;\n  } else if (indate >= q2[0] && indate <= q2[1]){\n    return 2;\n  } else if (indate >= q3[0] && indate <= q3[1]){\n    return 3;\n  } else if (indate >= q4[0] && indate <= q4[1]){\n    return 4;\n  } else {\n    throw new Error(\"Date out of range!\");\n  }\n}\n\nfunction evaluatePixel(samples, scenes) {\n  // Set date intervals\n  var q1 = [new Date(\"2020-01-01\"), new Date(\"2020-03-31\")];\n  var q2 = [new Date(\"2020-04-01\"), new Date(\"2020-06-30\")];\n  var q3 = [new Date(\"2020-07-01\"), new Date(\"2020-09-30\")];\n  var q4 = [new Date(\"2020-10-01\"), new Date(\"2020-12-31\")];\n  \n  // Initialise object with arrays\n  var VV_samples = {q1: [], q2: [], q3:[], q4:[]};\n\n  // Loop through orbits and add data\n  for (var i = 0; i < samples.length; i++) {\n    // Get quarter\n    quarter = check_quarter(new Date(scenes[i].date), q1, q2, q3, q4)\n    // Assign to correct array\n    if (samples[i].dataMask !=0){\n      if (quarter == 1){\n        VV_samples.q1.push(samples[i].VV);\n      } else if (quarter == 2){\n        VV_samples.q2.push(samples[i].VV);\n      } else if (quarter == 3){\n        VV_samples.q3.push(samples[i].VV);\n      } else if (quarter == 4){\n        VV_samples.q4.push(samples[i].VV);\n      }\n    }\n  }\n\n  const factor = 65535;\n\n  // Make an object to contain median values\n  var VV_median = {};\n  \n  // Compute median for each quarter\n  for (const key in VV_samples) {\n    if (VV_samples.hasOwnProperty(key)) {\n      if (VV_samples[key].length == 0){\n        VV_median[key] = [factor];\n      } else {\n        VV_median[key] = median(VV_samples[key]) * factor;\n      }\n    }\n  }\n\n  return {\n    VV: Object.values(VV_median)\n  }\n}"
}'

Let us know if you have any other questions!

Thanks so much Maxim! This works inasmuch as the values returned are no longer zero, but it still returns every date and not the median composite, at least on my end. Do you know how to set it up to return just the median? I can’t see what is still wrong

I have it set up like the following, which returns for Q1 (jan - march) for the following BBox:

Bbox: [17.77777822222222,-18.999999777777777,17.833333777777778,-18.94444422222222]
The following dates will be downloaded: [5, 17, 29, 41, 53, 64, 76, 88]
The S1 shape is: (8, 309, 292, 2) # This should be (309, 292, 2) not (8, ...)
year = 2020
dates_q1 = (f'{str(year)}-01-01' , f'{str(year)}-04-01')
dates_q2 = (f'{str(year)}-04-01' , f'{str(year)}-07-01')
dates_q3 = (f'{str(year)}-07-01' , f'{str(year)}-10-01')
dates_q4 = (f'{str(year)}-10-01' , f'{str(year)}-12-31')

for date in [dates_q1, dates_q2, dates_q3, dates_q4]:

    image_request = WcsRequest(
        layer=layer, bbox=box,
        time=date,
        image_format = MimeType.TIFF_d16,
        data_source=source, maxcc=1.0,
        resx='20m', resy='20m',
        instance_id=api_key,
        custom_url_params = {constants.CustomUrlParam.DOWNSAMPLING: 'NEAREST',
                            constants.CustomUrlParam.UPSAMPLING: 'NEAREST'},
        time_difference=datetime.timedelta(hours=72),
    )

    s1_dates_dict = [x for x in image_request.get_dates()]
    s1_dates = extract_dates(s1_dates_dict, year)
    print(f"The following dates will be downloaded: {s1_dates}")

    if len(image_request.download_list) >= 2:
        
            s1 = np.array(image_request.get_data())
            print(f"The S1 shape is: {s1.shape}")
            if not isinstance(s1.flat[0], np.floating):
                assert np.max(s1) > 1
                s1 = np.float32(s1) / 65535.
            assert np.max(s1) <= 1
            s1_all.append(s1)
s1 = np.stack(s1)
print(s1.shape)

Is there a particular reason for you to use WCS? If not I would suggest using the Process API for your call. I mimicked your approach:

evalscript = """
//VERSION3
function median(values) {
    if (values.length === 0) return 0;
    if (values.length === 1) return values[0];

    values.sort(function(a, b) {
        return a - b;
    });

    var half = Math.floor(values.length / 2);
    return values[half];
}


function evaluatePixel(samples) {
    // Initialise arrays
    var VV_samples = [];
    var VH_samples = [];
    
    // Loop through orbits and add data
    for (let i=0; i<samples.length; i++){
      // Ignore noData
      if (samples[i].dataMask != 0){
        VV_samples.push(samples[i].VV);
        VH_samples.push(samples[i].VH);
       }
    }  
    
    const factor = 65535;
    
    if (VV_samples.length == 0){
      var VV_median = [factor];
    } else {
      var VV_median = median(VV_samples) * factor;
    }
    if (VH_samples.length == 0){
      var VH_median = [factor];
    } else{
      var VH_median = median(VH_samples) * factor;
    }
    
    return [VV_median, VH_median];
}

function setup() {
  return {
    input: [{
      bands: [
        "VV",
        "VH",
        "dataMask",
             ]
    }],
    output: {
      bands: 2,
      sampleType:"UINT16"
    },
    mosaicking: "ORBIT"
  }
}
"""
# Set coords in WGS84
coords = [17.77777822222222,-18.999999777777777,17.833333777777778,-18.94444422222222]

# Set resolution
resolution = 10

# Make bbox and size
bbox = BBox(bbox=coords, crs=CRS.WGS84)
size = bbox_to_dimensions(bbox, resolution=resolution)


year = 2020
dates_q1 = (f'{str(year)}-01-01' , f'{str(year)}-04-01')
dates_q2 = (f'{str(year)}-04-01' , f'{str(year)}-07-01')
dates_q3 = (f'{str(year)}-07-01' , f'{str(year)}-10-01')
dates_q4 = (f'{str(year)}-10-01' , f'{str(year)}-12-31')

for date in [dates_q1, dates_q2, dates_q3, dates_q4]:
    
    request = SentinelHubRequest(
        evalscript=evalscript,
        input_data=[
            SentinelHubRequest.input_data(
                data_collection=DataCollection.SENTINEL1_IW,          
                time_interval=date,          
            ),
        ],
        responses=[
        SentinelHubRequest.output_response('default', MimeType.TIFF)
        ],
        bbox=bbox,
        size=size,
        config=config
    )

    response = request.get_data()
    print(f"Shape of returned images for {date} = {response[0].shape[0:2]}")
    print(f"Number of bands = {response[0].shape[-1]}")

As a side note, make sure to consider using ortorectify/gamma0/RTC/speckle options in order to get ARD…

@gmilcinski is there a way to integrate that into the SentinelHubRequest framework? I couldn’t figure out how to. The Process API works for this (thanks Maxim!) but now I don’t know how to specify terrain correction, etc.

Are you referring to our sentinelhub-py package?
@maleksandrov - do you perhaps know?

One just needs to pust another tag in the request, i.e.:

        "processing": {
            "backCoeff": "GAMMA0_ELLIPSOID",
            "orthorectify": "true"
        }

But I do not know, how to…

Hi @john.brandt,

You can provide these parameters to SentinelHubRequest.input_data method under parameter other_args. Altogether, your request can look something like this:

request = SentinelHubRequest(
    input_data=[
        SentinelHubRequest.input_data(
            ...,
            other_args={
                "processing": {
                    "backCoeff": "GAMMA0_ELLIPSOID",
                    "orthorectify": "true"
                }
            }
        ),
    ...
)