Change in brightness for cloudless mosaic and raw band values generated from a range of two months

I am trying to create a cloudless mosaic using the following script as adapted from the custom script repo with minor modification

//VERSION=3

// adapted from here:
// https://github.com/sentinel-hub/custom-scripts/blob/master/sentinel-2/cloudless_mosaic/L2A-first_quartile_4bands.js

function setup() {
  return {
    input: [{
      bands: [
        "B04", // red
        "B03", // green
        "B02", // blue
        "B08", // nir
        "B11", // swir
        "SCL" // pixel classification
      ],
      units: "DN"
    }],
    output: {
      bands: 5,
      sampleType: SampleType.UINT16
    },
    mosaicking: "ORBIT"
  };
}

function updateOutputMetadata(scenes,inputMetadata, outputMetadata) {
    outputMetadata.userData = { "scenes":  scenes.orbits }
}

function filterScenes(availableScenes, inputMetadata) {
  return availableScenes.filter(function(scene) {
    var m = scene.date.getMonth();
    var y = scene.date.getFullYear();
    var years = [2020];  #year 
    var months = [6, 7]; // 0 indexed for months
    return months.includes(m) && years.includes(y);
  });
}

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

// function for pulling first quartile of values
function getFirstQuartile(sortedValues) {
  var index = Math.floor(sortedValues.length / 4);
  return sortedValues[index];
}

// function for pulling median (second quartile) of values
function getMedian(sortedValues) {
  var index = Math.floor(sortedValues.length / 2);
  return sortedValues[index];
}


function validate(samples) {
  var scl = samples.SCL;

  if (scl === 3) { // SC_CLOUD_SHADOW
    return false;
  } else if (scl === 9) { // SC_CLOUD_HIGH_PROBA
    return false;
  } else if (scl === 8) { // SC_CLOUD_MEDIUM_PROBA
    return false;
  } else if (scl === 7) { // SC_CLOUD_LOW_PROBA
    // return false;
  } else if (scl === 10) { // SC_THIN_CIRRUS
    return false;
  } else if (scl === 11) { // SC_SNOW_ICE
    return false;
  } else if (scl === 1) { // SC_SATURATED_DEFECTIVE
    return false;
  } else if (scl === 2) { // SC_DARK_FEATURE_SHADOW
    // return false;
  }
  return true;
}

function evaluatePixel(samples, scenes) {
  var clo_b02 = [];
  var clo_b03 = [];
  var clo_b04 = [];
  var clo_b08 = [];
  var clo_b11 = [];
  var clo_b02_invalid = [];
  var clo_b03_invalid = [];
  var clo_b04_invalid = [];
  var clo_b08_invalid = [];
  var clo_b11_invalid = [];
  var a = 0;
  var a_invalid = 0;

  for (var i = 0; i < samples.length; i++) {
    var sample = samples[i];
    if (sample.B02 > 0 && sample.B03 > 0 && sample.B04 > 0 && sample.B08 > 0 && sample.B11 > 0) {
      var isValid = validate(sample);

      if (isValid) {
        clo_b02[a] = sample.B02;
        clo_b03[a] = sample.B03;
        clo_b04[a] = sample.B04;
        clo_b08[a] = sample.B08;
        clo_b11[a] = sample.B11;
        a = a + 1;
      } else {
        clo_b02_invalid[a_invalid] = sample.B02;
        clo_b03_invalid[a_invalid] = sample.B03;
        clo_b04_invalid[a_invalid] = sample.B04;
        clo_b08_invalid[a_invalid] = sample.B08;
        clo_b11_invalid[a_invalid] = sample.B11;
        a_invalid = a_invalid + 1;
      }
    }
  }

  var rValue;
  var gValue;
  var bValue;
  var nValue;
  var sValue;
  if (a > 0) {
    rValue = getValue(clo_b04);
    gValue = getValue(clo_b03);
    bValue = getValue(clo_b02);
    nValue = getValue(clo_b08);
    sValue = getValue(clo_b11);
  } else if (a_invalid > 0) {
    rValue = getValue(clo_b04_invalid);
    gValue = getValue(clo_b03_invalid);
    bValue = getValue(clo_b02_invalid);
    nValue = getValue(clo_b08_invalid);
    sValue = getValue(clo_b11_invalid);
  } else {
    rValue = 0;
    gValue = 0;
    bValue = 0;
    nValue = 0;
    sValue = 0;
  }
  return [rValue, gValue, bValue, nValue, sValue];
}

It works pretty well, however, I get an image for the range of two months and the image (RGB) that appears to be bright and enhanced, although I did not increase the brightness or use visualiser. What I wanted is to get raw values (median of the pixels ranging across the months of interest).

And I ran the simple script to get raw values as I use to:

//VERSION=3
    function setup() {
        return {
            input: [{
                bands: ["B02","B03","B04","B08"],
                units: "DN"
            }],
            output: {
                bands: 4,
                sampleType: "INT16"
            }
        };
    }
    function evaluatePixel(sample) {
        return [sample.B02,
                sample.B03,
                sample.B04,
                sample.B08];
    }

Note: The input data is same for both the scripts

"data": [
      {
        "dataFilter": {
          "timeRange": {
            "from": "2020-06-01T00:00:00Z",
            "to": "2020-07-26T23:59:59Z"
          },
          "maxCloudCoverage": 50,
          "mosaickingOrder": "leastRecent"
        },
        "type": "sentinel-2-l2a"
      }
    ]
  },

The RGB image I construct from the raw values is less bright than the former. I am not sure how I would explain this difference to myself. I was thinking maybe it was because of the mosaicking type.

Hi,

In the second script that you have provided (the simple script to get raw values), are you still expecting to return the median values from this script? Currently, it only uses simple mosaicking so will return the leastRecent acquisition pixels to you.

Therefore, it is not unexpected that the images look different because the returns are very different. One is a single acquisition and the other is a mosaic of several acquisitions.

Happy to help you get to the bottom of this though :slight_smile:

Thanks for the response @william.ray ! Alright, I get it. Thanks for the clarification. However, I am not able to execute the first script in Python, however, it works with Request builder. I am not sure why

bbox = BBox(bbox=[12.457598, 41.922684, 12.480912, 41.939893], crs=CRS.WGS84)

request = SentinelHubRequest(
    evalscript=evalscript,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A,          
            time_interval=('2022-07-28', '2022-08-28'),          
        ),
    ],
    responses=[
        SentinelHubRequest.output_response('default', MimeType.TIFF),
    ],
    bbox=bbox,
    resolution=(10,10),
    config=config
)

response = request.get_data()

The above script throws error if I specify resolution=(10,10) although the area is just aorun 4 km square. The error is as follows:

DownloadFailedException: Failed to download from:
https://services.sentinel-hub.com/api/v1/process
with HTTPError:
400 Client Error: Bad Request for url: https://services.sentinel-hub.com/api/v1/process
Server response: "{"error":{"status":400,"reason":"Bad Request","message":"Your request of 1929.15 meters per pixel exceeds the limit 1500.00 meters per pixel of the collection S2L2A. Please revise the resolution (or corresponding width/height) to make sure it is in supported range.","code":"RENDERER_EXCEPTION"}}"

I assume This limit is only for the process api and not for the request builder?

Hi Aswin,

Actually, the error you have here is related to the projection and resolution that you are using. You are using WGS84 in your request and therefore you need to use the correct units when specifying the resolution.

This will be in degrees which is the native units of WGS84. For approximately 10m resolution, I would recommend using 0.0001 degrees. Alternatively, you can request the data in the appropriate UTM zone and use metres as the resolution units.

Thanks for the response @william.ray ! I’ll try changing the CRS then

Hi @william.ray, I am getting a blank image when I run the following evalscript (please see the first evalscript in the question) in python :slight_smile:

bbox = BBox(bbox=[
  12.457598,
  41.922684,
  12.480912,
  41.939893
], crs=CRS.WGS84)

request = SentinelHubRequest(
    data_folder='/Users/amanohar/bid4best_internship/tree_species_classification/data/',
    evalscript=evalscript,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A,          
            time_interval=('2020-06-01', '2020-07-26'),          
            other_args={"dataFilter": {"maxCloudCoverage": 50,"mosaickingOrder": "leastRecent"}}
        ),
    ],
    responses=[
        SentinelHubRequest.output_response('default', MimeType.TIFF),
        SentinelHubRequest.output_response("userdata", MimeType.JSON),
    ],
    bbox=bbox,
    resolution=(0.0001,0.0001),
    config=config,
    #results_dir='/Users/amanohar/bid4best_internship/tree_species_classification/data/'
)

response = request.get_data(save_data=True)

However, the evalscript works in request builder and I get a nice image.

Please can you share the full Python code with us and we can help debug the code for you Aswin. Thanks.

Here is the python script:

import geopandas as gpd
import numpy as np
#upgrade shapely to 1.8.2
import datetime
import os

import matplotlib.pyplot as plt
import numpy as np

from math import radians, cos, sin, asin, sqrt
from sentinelhub import (
    CRS,
    BBox,
    DataCollection,
    DownloadRequest,
    MimeType,
    SentinelHubDownloadClient,
    SentinelHubRequest,
    bbox_to_dimensions,
)

from sentinelhub import SHConfig

config = SHConfig()

# Set up configuration tool

config.sh_client_id = sh_client
config.sh_client_secret = sh_secret

if not config.sh_client_id or not config.sh_client_secret:
    print("Warning! To use Process API, please provide the credentials (OAuth client ID and client secret).")

evalscript = """
//VERSION=3

// copied from here:
// https://github.com/sentinel-hub/custom-scripts/blob/master/sentinel-2/cloudless_mosaic/L2A-first_quartile_4bands.js

function setup() {
  return {
    input: [{
      bands: [
        "B04", // red
        "B03", // green
        "B02", // blue
        "B08", // nir
        "B11", // swir
        "SCL" // pixel classification
      ],
      units: "DN"
    }],
    output: {
      bands: 5,
      sampleType: SampleType.UINT16
    },
    mosaicking: "ORBIT"
  };
}

function updateOutputMetadata(scenes,inputMetadata, outputMetadata) {
    outputMetadata.userData = { "scenes":  scenes.orbits }
}

// acceptable images are ones collected in months of 2018, 2019, 2020 ...
function filterScenes(availableScenes, inputMetadata) {
  return availableScenes.filter(function(scene) {
    var m = scene.date.getMonth();
    var y = scene.date.getFullYear();
    var years = [2020];
    var months = [6, 7]; // 0 indexed, 5 = June, 6 = July, 7 = August
    return months.includes(m) && years.includes(y);
  });
}

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

// function for pulling first quartile of values
function getFirstQuartile(sortedValues) {
  var index = Math.floor(sortedValues.length / 4);
  return sortedValues[index];
}

// function for pulling median (second quartile) of values
function getMedian(sortedValues) {
  var index = Math.floor(sortedValues.length / 2);
  return sortedValues[index];
}

// adapted from https://forum.sentinel-hub.com/t/l2a-scene-classification-for-sentinel-2/51

function validate(samples) {
  var scl = samples.SCL;

  if (scl === 3) { // SC_CLOUD_SHADOW
    return false;
  } else if (scl === 9) { // SC_CLOUD_HIGH_PROBA
    return false;
  } else if (scl === 8) { // SC_CLOUD_MEDIUM_PROBA
    return false;
  } else if (scl === 7) { // SC_CLOUD_LOW_PROBA
    // return false;
  } else if (scl === 10) { // SC_THIN_CIRRUS
    return false;
  } else if (scl === 11) { // SC_SNOW_ICE
    return false;
  } else if (scl === 1) { // SC_SATURATED_DEFECTIVE
    return false;
  } else if (scl === 2) { // SC_DARK_FEATURE_SHADOW
    // return false;
  }
  return true;
}

function evaluatePixel(samples, scenes) {
  var clo_b02 = [];
  var clo_b03 = [];
  var clo_b04 = [];
  var clo_b08 = [];
  var clo_b11 = [];
  var clo_b02_invalid = [];
  var clo_b03_invalid = [];
  var clo_b04_invalid = [];
  var clo_b08_invalid = [];
  var clo_b11_invalid = [];
  var a = 0;
  var a_invalid = 0;

  for (var i = 0; i < samples.length; i++) {
    var sample = samples[i];
    if (sample.B02 > 0 && sample.B03 > 0 && sample.B04 > 0 && sample.B08 > 0 && sample.B11 > 0) {
      var isValid = validate(sample);

      if (isValid) {
        clo_b02[a] = sample.B02;
        clo_b03[a] = sample.B03;
        clo_b04[a] = sample.B04;
        clo_b08[a] = sample.B08;
        clo_b11[a] = sample.B11;
        a = a + 1;
      } else {
        clo_b02_invalid[a_invalid] = sample.B02;
        clo_b03_invalid[a_invalid] = sample.B03;
        clo_b04_invalid[a_invalid] = sample.B04;
        clo_b08_invalid[a_invalid] = sample.B08;
        clo_b11_invalid[a_invalid] = sample.B11;
        a_invalid = a_invalid + 1;
      }
    }
  }

  var rValue;
  var gValue;
  var bValue;
  var nValue;
  var sValue;
  if (a > 0) {
    rValue = getValue(clo_b04);
    gValue = getValue(clo_b03);
    bValue = getValue(clo_b02);
    nValue = getValue(clo_b08);
    sValue = getValue(clo_b11);
  } else if (a_invalid > 0) {
    rValue = getValue(clo_b04_invalid);
    gValue = getValue(clo_b03_invalid);
    bValue = getValue(clo_b02_invalid);
    nValue = getValue(clo_b08_invalid);
    sValue = getValue(clo_b11_invalid);
  } else {
    rValue = 0;
    gValue = 0;
    bValue = 0;
    nValue = 0;
    sValue = 0;
  }
  return [rValue, gValue, bValue, nValue, sValue];
}
"""
bbox = BBox(bbox=[
  12.457598,
  41.922684,
  12.480912,
  41.939893
], crs=CRS.WGS84)

request = SentinelHubRequest(
    data_folder='/Users/amanohar/bid4best_internship/tree_species_classification/data/',
    evalscript=evalscript,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A,          
            time_interval=('2020-06-01', '2020-07-26'),          
            other_args={"dataFilter": {"maxCloudCoverage": 50,"mosaickingOrder": "leastRecent"}}
        ),
    ],
    responses=[
        SentinelHubRequest.output_response('default', MimeType.TIFF),
        SentinelHubRequest.output_response("userdata", MimeType.JSON),
    ],
    bbox=bbox,
    resolution=(0.0001,0.0001),
    config=config,
    #results_dir='/Users/amanohar/bid4best_internship/tree_species_classification/data/'
)

response = request.get_data(save_data=True)

Thanks for the quick response @william.ray and for the help.

Hi,

It seems that your filterScenes function was causing the error in the return. Reinserting the original function makes the script work properly again. Are you able to explain what you were trying to do differently?

In addition, was the evalscript in Request Builder exactly the same or different? Might be worth double checking this.

The below script will work for you:

//VERSION=3

// copied from here:
// https://github.com/sentinel-hub/custom-scripts/blob/master/sentinel-2/cloudless_mosaic/L2A-first_quartile_4bands.js

function setup() {
    return {
      input: [{
        bands: [
          "B04", // red
          "B03", // green
          "B02", // blue
          "B08", // nir
          "B11", // swir
          "SCL" // pixel classification
        ],
        units: "DN"
      }],
      output: {
        bands: 5, sampleType:"UINT16"
      },
      mosaicking: "ORBIT"
    };
  }
  
  function updateOutputMetadata(scenes,inputMetadata, outputMetadata) {
      outputMetadata.userData = { "scenes":  scenes.orbits }
  }
  
  // acceptable images are ones collected in months of 2018, 2019, 2020 ...
  function filterScenes (scenes, inputMetadata) {
    return scenes.filter(function (scene) {
      return scene.date.getTime()>=(inputMetadata.to.getTime()-12*31*24*3600*1000);
    });
  }
  
  function getValue(values) {
    values.sort(function (a, b) {
      return a - b;
    });
    return getMedian(values);
  }
  
  // function for pulling first quartile of values
  function getFirstQuartile(sortedValues) {
    var index = Math.floor(sortedValues.length / 4);
    return sortedValues[index];
  }
  
  // function for pulling median (second quartile) of values
  function getMedian(sortedValues) {
    var index = Math.floor(sortedValues.length / 2);
    return sortedValues[index];
  }
  
  // adapted from https://forum.sentinel-hub.com/t/l2a-scene-classification-for-sentinel-2/51
  
  function validate(samples) {
    var scl = samples.SCL;
  
    if (scl === 3) { // SC_CLOUD_SHADOW
      return false;
    } else if (scl === 9) { // SC_CLOUD_HIGH_PROBA
      return false;
    } else if (scl === 8) { // SC_CLOUD_MEDIUM_PROBA
      return false;
    } else if (scl === 7) { // SC_CLOUD_LOW_PROBA
      // return false;
    } else if (scl === 10) { // SC_THIN_CIRRUS
      return false;
    } else if (scl === 11) { // SC_SNOW_ICE
      return false;
    } else if (scl === 1) { // SC_SATURATED_DEFECTIVE
      return false;
    } else if (scl === 2) { // SC_DARK_FEATURE_SHADOW
      // return false;
    }
    return true;
  }
  
  function evaluatePixel(samples, scenes) {
    var clo_b02 = [];
    var clo_b03 = [];
    var clo_b04 = [];
    var clo_b08 = [];
    var clo_b11 = [];
    var clo_b02_invalid = [];
    var clo_b03_invalid = [];
    var clo_b04_invalid = [];
    var clo_b08_invalid = [];
    var clo_b11_invalid = [];
    var a = 0;
    var a_invalid = 0;
  
    for (var i = 0; i < samples.length; i++) {
      var sample = samples[i];
      if (sample.B02 > 0 && sample.B03 > 0 && sample.B04 > 0 && sample.B08 > 0 && sample.B11 > 0) {
        var isValid = validate(sample);
  
        if (isValid) {
          clo_b02[a] = sample.B02;
          clo_b03[a] = sample.B03;
          clo_b04[a] = sample.B04;
          clo_b08[a] = sample.B08;
          clo_b11[a] = sample.B11;
          a = a + 1;
        } else {
          clo_b02_invalid[a_invalid] = sample.B02;
          clo_b03_invalid[a_invalid] = sample.B03;
          clo_b04_invalid[a_invalid] = sample.B04;
          clo_b08_invalid[a_invalid] = sample.B08;
          clo_b11_invalid[a_invalid] = sample.B11;
          a_invalid = a_invalid + 1;
        }
      }
    }
  
    var rValue;
    var gValue;
    var bValue;
    var nValue;
    var sValue;
    if (a > 0) {
      rValue = getValue(clo_b04);
      gValue = getValue(clo_b03);
      bValue = getValue(clo_b02);
      nValue = getValue(clo_b08);
      sValue = getValue(clo_b11);
    } else if (a_invalid > 0) {
      rValue = getValue(clo_b04_invalid);
      gValue = getValue(clo_b03_invalid);
      bValue = getValue(clo_b02_invalid);
      nValue = getValue(clo_b08_invalid);
      sValue = getValue(clo_b11_invalid);
    } else {
      rValue = 0;
      gValue = 0;
      bValue = 0;
      nValue = 0;
      sValue = 0;
    }
    return [rValue, gValue, bValue, nValue, sValue];
  }

Thanks for the response @william.ray! I was trying to specify the months and years within the script to not bother changing the python function to search for a certain range of months. Thought of making it standard to get only a two months mosaic. But this should work as well I think, however, my trail period is over and I cannot test this until I sometime.

I have another related question, do have a standard method to make RGB images without losing much on manipulating contrasts and brightness, which is inconsistent for different images.

I can refer you to this forum thread so I would watch this space!

Ah, okay, It is my thread. I will watch out then