Cloudless Sentinel-2 image sourced from a single orbit

I have been working on obtaining cloudless mosaic images but have been struggling with strong orbit-wise artifacts in some areas (see my post).

As an alternative to seamless mosaics I am interested in obtaining orbit-wise cloudless images. The cloud-filling results in my tests so far look really good but the stripes introduced by the mosaicking process are going to complicate downstream applications of the imagery.

Is there a way to specify a request to only pull images from a particular orbit rather than sourcing images from multiple orbits?

For example, in the image below I would only want to pull raw images from the center orbit (relative orbit 54) and avoid pulling images from relative orbits 97 and 11.

Is the Sentinel product ID available in the scenes object? If it is then I could parse the relative orbit from the product ID:

This scene is in relative orbit 31:
S2A_MSIL1C_20170105T013442_N0204_R031_T53NMJ_20170105T013443

I have been combing through the docs and it looks like orbitId is supposed to be carried in the metadata for scenes.

I have tried to add an orbit filter in my filterScenes function, but that is not yielding any scenes at all in my tests. Is the orbitId value the relative or absolute orbit?

\\ function to convert absolute orbit to relative orbit
function getRelativeOrbit(o) {
  return o % 143 + 3;
}

function filterScenes(availableScenes, inputMetadata) {
  return availableScenes.filter(function(scene) {
    var m = scene.date.getMonth();
    var y = scene.date.getFullYear();
    var o = scene.orbitId
    var years = [2018, 2019, 2020];
    var months = [5, 6, 7]; // 0 indexed, 5 = June, 6 = July, 7 = August
    var orbit = 54; // target relative orbit
    return months.includes(m) && years.includes(y) && orbit == o;
  });
}

Is there a better way to debug the filterScenes behavior? My process has been to make adjustments then order a low resolution preview image but when I get an empty image back I can tell that the filtering is not working as I intended.

The metadata returned following this example does not contain orbitId, so it seems as if my filter was doomed!

{\"date\":\"2020-08-28T00:00:00.000Z\",\"idx\":0,\"bandBuffers\":[{},{},{},{},{},{}]}, ...

Is there some way to pull more metadata from the scenes?

Would you be so kind and put in here an example of the whole request, so that we can debug it?

Also, can you try setting

 mosaicking: "TILE"

(instead of ORBIT)

You have ventured in quite advanced use of our API (which is good thing!), so it might be that some things are more buggy than what would have been needed…

Here is the request (copied out of Request Builder). The process does appear to execute without any errors but the image that is returned from the low-resolution preview is all no-data, so it seems like all collections are getting filtered out by the query.

curl -X POST https://services.sentinel-hub.com/api/v1/batch/process/ 
 -H 'Content-Type: application/json' 
 -H 'Authorization: Bearer [authorization string]' 
 -d '{
  "processRequest": {
    "input": {
      "bounds": {
        "bbox": [
          -75.602668,
          43.755225,
          -73.494418,
          44.059959
        ],
        "properties": {
          "crs": "http://www.opengis.net/def/crs/EPSG/0/4326"
        }
      },
      "data": [
        {
          "type": "S2L2A",
          "dataFilter": {
            "timeRange": {
              "from": "2018-06-01T00:00:00Z",
              "to": "2020-09-15T23:59:59Z"
            },
            "maxCloudCoverage": 50
          }
        }
      ]
    },
    "output": {
      "responses": [
        {
          "identifier": "default",
          "format": {
            "type": "image/tiff"
          }
        }
      ]
    },
    "evalscript": "//VERSION=3\n\n// copied from here:\n// https://github.com/sentinel-hub/custom-scripts/blob/master/sentinel-2/cloudless_mosaic/L2A-first_quartile_4bands.js\n\nfunction setup() {\n  return {\n    input: [{\n      bands: [\n        \"B04\", // red\n        \"B03\", // green\n        \"B02\", // blue\n        \"B08\", // nir\n        \"B11\", // swir\n        \"SCL\" // pixel classification\n      ],\n      units: \"DN\"\n    }],\n    output: {\n      bands: 5,\n      sampleType: SampleType.UINT16\n    },\n    mosaicking: \"ORBIT\"\n  };\n}\n\n// function to convert absolute orbit to relative orbit\nfunction getRelativeOrbit(o) {\n  return o % 143 + 3;\n}\n\n// acceptable images are ones collected in summer months of 2018, 2019, 2020\nfunction filterScenes(availableScenes, inputMetadata) {\n  return availableScenes.filter(function(scene) {\n    var m = scene.date.getMonth();\n    var y = scene.date.getFullYear();\n    var o = getRelativeOrbit(scene.orbitId);\n    var years = [2018, 2019, 2020];\n    var months = [5, 6, 7]; // 0 indexed, 5 = June, 6 = July, 7 = August\n    var orbits = [54];\n    return months.includes(m) && years.includes(y) && orbits.includes(o);\n  });\n}\n\nfunction getValue(values) {\n  values.sort(function (a, b) {\n    return a - b;\n  });\n  return getMedian(values);\n}\n\n// function for pulling first quartile of values\nfunction getFirstQuartile(sortedValues) {\n  var index = Math.floor(sortedValues.length / 4);\n  return sortedValues[index];\n}\n\n// function for pulling median (second quartile) of values\nfunction getMedian(sortedValues) {\n  var index = Math.floor(sortedValues.length / 2);\n  return sortedValues[index];\n}\n\n\nfunction validate(samples) {\n  var scl = samples.SCL;\n\n  if (scl === 3) { // SC_CLOUD_SHADOW\n    return false;\n  } else if (scl === 9) { // SC_CLOUD_HIGH_PROBA\n    return false;\n  } else if (scl === 8) { // SC_CLOUD_MEDIUM_PROBA\n    return false;\n  } else if (scl === 7) { // SC_CLOUD_LOW_PROBA\n    // return false;\n  } else if (scl === 10) { // SC_THIN_CIRRUS\n    return false;\n  } else if (scl === 11) { // SC_SNOW_ICE\n    return false;\n  } else if (scl === 1) { // SC_SATURATED_DEFECTIVE\n    return false;\n  } else if (scl === 2) { // SC_DARK_FEATURE_SHADOW\n    // return false;\n  }\n  return true;\n}\n\nfunction evaluatePixel(samples, scenes) {\n  var clo_b02 = [];\n  var clo_b03 = [];\n  var clo_b04 = [];\n  var clo_b08 = [];\n  var clo_b11 = [];\n  var clo_b02_invalid = [];\n  var clo_b03_invalid = [];\n  var clo_b04_invalid = [];\n  var clo_b08_invalid = [];\n  var clo_b11_invalid = [];\n  var a = 0;\n  var a_invalid = 0;\n\n  for (var i = 0; i < samples.length; i++) {\n    var sample = samples[i];\n    if (sample.B02 > 0 && sample.B03 > 0 && sample.B04 > 0 && sample.B08 > 0 && sample.B11 > 0) {\n      var isValid = validate(sample);\n\n      if (isValid) {\n        clo_b02[a] = sample.B02;\n        clo_b03[a] = sample.B03;\n        clo_b04[a] = sample.B04;\n        clo_b08[a] = sample.B08;\n        clo_b11[a] = sample.B11;\n        a = a + 1;\n      } else {\n        clo_b02_invalid[a_invalid] = sample.B02;\n        clo_b03_invalid[a_invalid] = sample.B03;\n        clo_b04_invalid[a_invalid] = sample.B04;\n        clo_b08_invalid[a_invalid] = sample.B08;\n        clo_b11_invalid[a_invalid] = sample.B11;\n        a_invalid = a_invalid + 1;\n      }\n    }\n  }\n\n  var rValue;\n  var gValue;\n  var bValue;\n  var nValue;\n  var sValue;\n  if (a > 0) {\n    rValue = getValue(clo_b04);\n    gValue = getValue(clo_b03);\n    bValue = getValue(clo_b02);\n    nValue = getValue(clo_b08);\n    sValue = getValue(clo_b11);\n  } else if (a_invalid > 0) {\n    rValue = getValue(clo_b04_invalid);\n    gValue = getValue(clo_b03_invalid);\n    bValue = getValue(clo_b02_invalid);\n    nValue = getValue(clo_b08_invalid);\n    sValue = getValue(clo_b11_invalid);\n  } else {\n    rValue = 0;\n    gValue = 0;\n    bValue = 0;\n    nValue = 0;\n    sValue = 0;\n  }\n  return [rValue, gValue, bValue, nValue, sValue];\n}\n"
  },
  "tilingGrid": {
    "id": 0,
    "resolution": 20
  },
  "bucketName": "silviaterra-sentinel-hub"
}'

I tried setting mosacking: "TILE" but that returned this error:

Failed to evaluate script! evalscript.js:35: TypeError: Cannot read property 'getMonth' of undefined var m = scene.date.getMonth(); ^ TypeError: Cannot read property 'getMonth' of undefined at evalscript.js:35:23 at Array.filter (native) at filterScenes (evalscript.js:34:26)

I assumed that the orbitId values in the scene metadata are absolute orbits (not relative orbits). @gmilcinski can you confirm which type are supposed to be stored in the metadata?

Hi @hrodman,

with filterScenes you can not atm filter out the scenes based on the orbit id. So, you will have to pull all the data into your evaluatePixel function but there are original tile ids available in scenes object, e.g. try scenes.tiles[0].tileOriginalId (mosaicking TILE), from where you can parse absolute orbit id (check Granule naming convention here).

I hope that helps. If not, let us know.

Two other comments:

  1. You can also write to userdata.json file from within the evaluatePixel function, which makes it easy to check what properties of scenes are available, e.g. this is slightly adjusted example from the documentation:
    curl -X POST \
      https://services.sentinel-hub.com/api/v1/process \
      -H 'Authorization: Bearer <your access token>' \
      -H 'accept: application/tar' \
      -F 'request={
        "input": {
            "bounds": {
                "bbox": [
                    13.822174072265625,
                    45.85080395917834,
                    14.55963134765625,
                    46.29191774991382
                ]
            },
            "data": [
                {
                    "type": "S2L2A",
                    "dataFilter": {
                        "timeRange": {
                            "from": "2018-12-27T00:00:00Z",
                            "to": "2018-12-27T23:59:59Z"
                        }
                    }
                }
            ]
        },
        "output": {
            "width": 200,
            "height": 100,
            "responses": [
                {
                    "identifier": "default",
                    "format": {
                        "type": "image/png"
                    }
                },
                {
                    "identifier": "userdata",
                    "format": {
                        "type": "application/json"
                    }
                }
            ]
        }
    }' \
      -F 'evalscript=//VERSION=3
    function setup() {
      return {
        input: ["B02", "B03", "B04"],
        mosaicking: Mosaicking.TILE,
        output: { id:"default", bands: 3}
      }
    }

    function evaluatePixel(samples, scenes, inputMetadata, customData, outputMetadata)) {
      outputMetadata.userData = { "metadata":  scenes.tiles }
      return [ 2.5*samples[0].B04, 2.5*samples[0].B03, 2.5*samples[0].B02 ]
    }'
  1. We have discussed your request/question and we realize that it would be useful to have an option to filter by orbit id even before the processing reaches evaluatePixel function. We will look into that and I will post any updates about it here.

Best, Anja

Thank you for the recommendation, @avrecko. I wound up implementing a process that uses the WebFeatureService function from the python API to query for all dates for a given bounding box then parses the scene IDs (granule naming convention) to filter down to the desired relative orbits. I can add the list of exact dates that I want to include to my evalscript so only images from the desired relative orbits make it past filterScenes.

Here are the python functions that I wrote to get the job done:

from sentinelhub import BBox, bbox_to_dimensions, CRS, DataCollection, \
    Geometry, MimeType, SentinelHubBatch, SentinelHubRequest, SHConfig, \
    WebFeatureService

# absolute orbits translate to relative orbits differently for 2A and 2B
# I arrived at this solution after some trial and error...
def absolute_to_relative_orbit(x, sat):
    assert sat in ['2A', '2B']
    if sat == '2A':
        adj = -140
    if sat == '2B':
        adj = -26

    return (x + adj) % 143

def get_dates_by_orbit(bbox, start_date, end_date, max_cc, target_orbit, config):
    assert target_orbit is not None, "relative_orbit must be specified"

    # convert target_orbit to list if just a single orbit
    if type(target_orbit) is int:
        target_orbit = [target_orbit]

    # define time window
    search_time_interval = (f'{start_date}T00:00:00', f'{end_date}T23:59:59')

    # query scenes
    wfs_iterator = WebFeatureService(
        bbox,
        search_time_interval,
        data_collection=DataCollection.SENTINEL2_L2A,
        maxcc=max_cc,
        config=config
    )

    # filter down to dates from specified orbit(s)
    dates = []
    for tile_info in wfs_iterator:
        # raw product ID
        id = tile_info['properties']['id']

        # acquisition date
        date = tile_info['properties']['date']

        # absolute orbit is buried in ID after _A string
        absolute_orbit = int(re.search("(?<=_A)[\d]+(?=_)", id).group())

        # satellite (either 2A or 2B) is defined at beginning of id
        sat = id[1:3]

        # convert to relative orbit
        relative_orbit = absolute_to_relative_orbit(absolute_orbit, sat)

        if relative_orbit not in target_orbit:
            continue

        # add date if not already added to list
        if date not in dates:
            dates.append(date)

    assert len(dates) > 0, \
        f'No dates available for this bounding box and relative orbit {target_orbit}'

    return dates

def filter_dates(dates, months, years):
    # convert date strings to date objects
    dates = [dt.datetime.strptime(date, '%Y-%m-%d').date() for date in dates]

    # filter down to supplied months/years
    filtered = [str(date) for date in dates if date.month in months and date.year in years]

    assert len(filtered) > 0, \
        'None of supplied dates satisfy desired months/years'

    return filtered

The logic in the absolute to relative orbit function seems fragile but it has worked in my tests so far.

Assuming you have a Geometry object called geom, a config object with credentials, and a evalscript string that is waiting for a string of dates to pass to filterScenes, you can implement the orbit-wise filter like this:

dates = get_dates_by_orbit(
    geom.bbox,
    start_date='2019-06-01',
    end_date='2020-09-01',
    max_cc=0.5,
    target_orbit = [126, 83, 40], # can be one or multiple orbits
    config=config)

# filter down to summer only
dates_filt = filter_dates(dates, months=[6, 7, 8], years=[2019, 2020])

# make list of dates into string to pass to evalscript
date_string = ', '.join([f'"{date}"' for date in dates_filt])
eval_formatted = evalscript % date_string

eval_formatted has the explicit list of dates in the filterScenes function and can be passed to SentinelHubRequest

My evalscript is stored as a python string, the filterScenes function looks like this:

'''
function filterScenes(availableScenes, inputMetadata) {
  var allowedDates = [%s]; // format with python
  return availableScenes.filter(function (scene) {
    var sceneDateStr = scene.date.toISOString().split("T")[0]; //converting date and time to string and rounding to day precision
    return allowedDates.includes(sceneDateStr);
  });
}
'''
1 Like