Define original image ID for each pixel of Quality Mosaic in Google Earth Engine

Question

I used a set of LS8 TOA images to produce a quality mosaic based on NDVI values. I wonder if there is a way do figure out a source / origin image ID for each pixel of the result mosaic. I guess it should be some kind of a thematic raster of the same resolution.

Source code I used to produce mosaic is attached below.

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');
// Get the least cloudy image in 2015.
var image = ee.Image(
  l8.filterBounds(ee.Geometry.Point(-122.2929, 37.8206))
    .filterDate('2015-01-01', '2015-12-31')
    .sort('CLOUD_COVER')
    .first()
);
Map.setCenter(-122.2929, 37.8206, 10)
// Compute the Normalized Difference Vegetation Index (NDVI).
var nir = image.select('B5');
var red = image.select('B4');
var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};
var withNDVI = l8.map(addNDVI);
// Make a "greenest" pixel composite.
var greenest = withNDVI.qualityMosaic('NDVI');
// Display the result.
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};
Map.addLayer(greenest, visParams, 'Greenest pixel composite');

Answer 1

I had the same problem when I created the 'Best Available Pixel' code (https://github.com/fitoprincipe/geebap), and I solved this way:

1. Add a band in which the value of every pixel is the number of days since the Epoch (1970-01-01T00:00:00Z)
2. Add a band with an encoded number for each collection and a property exposing the relation. For example, Landsat 5 TOA: 1, L7TOA: 2, etc.

This way you will have the date and the collection in each pixel, and therefor be able to 'track' the pixel down to its root. It should work fine for Landsat and Sentinel. I am not sure about MODIS because there is one image per day (should work as well).

Of course this is just my approach and there could be better ones. May be bits codification. For you, the code may be something like:

var tools = require('users/fitoprincipe/geetools:tools')

var get_days = function(date) {
  var m = date.millis()
  return m.divide(1000).divide(3600).toInt()
}

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
         .filterBounds(ee.Geometry.Point(-122.2929, 37.8206))
         .filterDate('2015-01-01', '2015-12-31');

var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

// add a 'date' band: number of days since epoch
var addDate = function(img) {
  var d = ee.Date(img.date())
  var days = get_days(d)
  var days_img = ee.Image.constant(days).rename('date').toInt32()
  return img.addBands(days_img)
}

var withNDVI = l8.map(addNDVI).map(addDate);
// Make a "greenest" pixel composite.
var greenest = withNDVI.qualityMosaic('NDVI');
// Display the result.
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};
Map.addLayer(greenest, visParams, 'Greenest pixel composite');

// Visualize the 'date' image
Map.addLayer(greenest.select('date'), {min:get_days(ee.Date('2015-01-01')).getInfo(),
                                       max:get_days(ee.Date('2015-12-31')).getInfo()}, 'date image')

// print out the Image Id from where the pixel came
// At this point you could add the original image to the Map
var getId = function(obj) {
  var lat = obj.lat
  var lon = obj.lon
  var point = ee.Geometry.Point([lon, lat])
  var date = tools.get_value(greenest, point, 30).get('date')
  date = ee.Date(ee.Number(date).multiply(1000).multiply(3600))
  var img = ee.Image(l8.filterDate(date, date.advance(1, 'day')).first())
  print('Img id in point ['+lon+","+lat+"] is", img.id())
}

Map.onClick(getId)

When you click on the Map, the Image ID will be printed on the console (may not work if you are in the inspector mode, switch to the console)