Sentinel-2 L2A multi-temporal stack

Apps / Services Python SDK
1

Dear SH forum,

I am currently working with Sentinel-2 L2A data to create a multi-temporal stack based on a box and a time interval, with specific features. I briefly report the requirements the stack should have:

  • 224 x 224 pixels for the final patch
  • 20 meters spatial resolution
  • Masking of clouds/cirrus/cloud shadow/saturated pixels/pixels whose value exceeds 1 (i.e. set the related pixels’ value to 0)
  • Spectral bands of interest B1, B2, B3, B4, B8A, B11, B12
  • Output format: 1 GeoTiff file for each sensing date and for each spectral band (e.g. <sensing_date1>_B1.tif, <sensing_date1>_B2.tif, …, <sensing_date2>_B1.tif, <sensing_date2>_B2.tif, etc…)

I started testing the sentinelhub python library for the purpose, exploiting the WcsRequest. Everything appears to work fine but I didn’t manage to solve a couple of points:

  • The masking of pixels whose value exceeds 1 (i.e. when a pixel has a value greater than 1, set it to 0)
  • The generation of a single GeoTiff for each sensing date and spectral band (by default, I obtain a multi-band GeoTiff for each sensing date containing all the requested spectral bands)

For the first case, I suppose it is necessary to set a condition within the evaluatePixel function but, since I am not familiar with JS, I am having some troubles getting this done.
For the second case, I tried to find a solution in the documentation but I was not able to find anything that could help me.

Is there a solution for these two points?

I add below the code of the WCSRequest along with the EvalScript V3 used

wcs_request = WcsRequest(data_collection=DataCollection.SENTINEL2_L2A,
                               image_format=MimeType.TIFF,
                               data_folder = data_folder_path,
                               layer='AIDEO-LAYER',
                               bbox=ip_box,
                               time=(start_str, end_str),
                               resx='{}m'.format(spatial_res),
                               resy='{}m'.format(spatial_res),
                               config=config_aideo)
wcs_request.save_data()

//VERSION=3

function evaluatePixel(sample) {
  	if ([0, 1, 2, 3, 8, 9, 10].includes(sample.SCL) ){
      return [0, 0, 0, 0, 0, 0, 0]
    } else{
      	return [sample.B01, sample.B02, sample.B03, sample.B04, sample.B8A, sample.B11,sample.B12];
    }
}

function setup() {
  return {
    input: [{
      bands: [
        "B01",
        "B02",
        "B03",
        "B04",
        "B8A",
        "B11",
        "B12",
        "SCL"
      ]
    }],
    output: {
      bands: 7,
      sampleType:"FLOAT32"
    }
  }
}
2

Hi Pietro,

By explaining your issues so clearly it makes it so much easier for us to help you: thank you! I’m going to answer your questions in the reverse order because it will be easier to explain in steps.

1. Save 1 Geotiff / band

The generation of a single GeoTiff for each sensing date and spectral band (by default, I obtain a multi-band GeoTiff for each sensing date containing all the requested spectral bands)

If you are not tied to using OGC services in Python, you can easily return a single GeoTiff per band using SentinelHubRequest. If you want a file per band using the WCSRequest, then I think you need to specify 1 layer per band in your instance.

With WCSRequest

In this case I would create as many layers as bands that you want returned, each with the same evalscript (e.g for Band 1, Band 2 etc…). Then in Python I would loop over the layers that you want to retrieve and run the WCSRequest.

With SentinelHubRequest

To export the results for each spectral band to a GeoTiff, you can make the most of the multiple output objects. You would then specify 1 output / band. Your setup function in the Evalscript would then look like this (I’ve reduced the number of bands to keep it short):

function setup() {
  return {
    input: [{
      bands: [
        "B01",
        "B02",
        "B03"
      ]
    }],
    output: [{id: "B1",
      bands: 1,
      sampleType: "FLOAT32"
    },
    {
      id: "B2",
      bands: 1,
      sampleType: "FLOAT32"
     },
     {
      id: "B3",
      bands: 1,
      sampleType: "FLOAT32"
     }
     ]
  }

Then you set the id of your product in the evalscript:

function evaluatePixel(sample) {
  	<do some things>
    return {
      B1: [sample.B01],
      B2: [sample.B02],
      B3: [sample.B03]
    }
}

Note that you are calling the identifiers from the output in the setup function (e.g. B1, B2…).

Then you set the responses as in the example here: https://sentinelhub-py.readthedocs.io/en/latest/examples/processing_api_request.html#Example-6-:-Multi-response-request-type

By the way, is there a specific reason that you return FLOAT32 values (such as the returned bands being your final product)? If not, a trick to save on data volume is to return your (band * 10000) in the evalscript and set the output sampleType to UINT16. Then just divide the returned tiffs by 10000 later.

2. Mask pixels > 1

The masking of pixels whose value exceeds 1 (i.e. when a pixel has a value greater than 1, set it to 0)

Now there is most probably a more elegant solution to do this, but I would write the following:

function evaluatePixel(sample) {
    // Initialise an array to contain results
    var output_bands = [];
  	if ([0, 1, 2, 3, 8, 9, 10].includes(sample.SCL) ){
      output_bands = [0, 0, 0]
    } else{
      	var all_bands =  [sample.B01, sample.B02, sample.B03];
        // Set values to zero if > 1 in the array
        output_bands = all_bands.map(function(item) { return item > 1 ? 0 : item; });
   
     return {
          B1: [output_bands[0]],
          B2: [output_bands[1]],
          B3: [output_bands[2]]
    }
  }

If you are using 1 layer / band with WcsRequest you can adjust the evalscript for each band:

function evaluatePixel(sample) {
  	if ([0, 1, 2, 3, 8, 9, 10].includes(sample.SCL) ){
        var band_result = 0
    } else{
        if (sample.B01 < 1){
      	    var band_result = sample.B01
        } else {
            var band_result = 0
        }
    }
        
     return [band_result]
}

Hope this helps,

Maxim