Google Earth Engine - get true color (RGB) of Landsat 8 SR

Question

I try to export RGB images from ee.ImageCollection('LANDSAT/LC08/C01/T1_SR'). As explained in the documentation, I select my region of interest and the three bands related to RGB.

To transform my ImageCollection to a single image, I apply the median reduction on that collection. When it is done I compute the percentiles of the image to have the minimum value and the maximum value for each band. When it is done I put the values in the VisParams.

Is my methodology is correct to have the true color of the SR image ?

In tutorials given the value is hard-coded (generally min: 0 and max: 3000). the other possibility is stretch method but this method only takes into account one value (which is the max of min quantile and max and max quantile).

As result the color image sightly change. Each method differ and the color as well and there is no information that discuss the true color of the image. I think defining the value for each band is the more realistic approach but I am not an expert in that field.

Note: The second method that use stretch is not really feasible in my case since the final code is written in python, thus 'tricks' that can be used on the browser editor does not always work with Python API.

Javascript code:

 var geometry = ee.Geometry.Polygon([[39.05789266 , 13.59051553],
                                     [39.11335033 , 13.59051553],
                                     [39.11335033 , 13.64477783],
                                     [39.05789266 , 13.64477783],
                                     [39.05789266 , 13.59051553]]);

/**
 * Function to mask clouds based on the pixel_qa band of Landsat 8 SR data.
 * @param {ee.Image} image input Landsat 8 SR image
 * @return {ee.Image} cloudmasked Landsat 8 image
 */
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterBounds(geometry)
                  .filterDate('2016-01-01', '2016-12-31')
                  .map(maskL8sr)
                  .select(['B4', 'B3','B2']);

var image = dataset.median()
var percentiles = image.reduceRegion(ee.Reducer.percentile([0, 100], ['min', 'max']), geometry,30).getInfo();
print(percentiles)
var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: [percentiles['B4_min'], percentiles['B3_min'], percentiles['B2_min']],
  max: [percentiles['B4_max'], percentiles['B3_max'], percentiles['B2_max']],
  gamma: 1
};
Map.centerObject(geometry);
Map.addLayer(image, visParams);

Python code:

import ee

ee.Initialize()
SATELLITE_SR = 'LANDSAT/LC08/C01/T1_SR'
RGB = ['B4', 'B3', 'B2']
geometry = ee.Geometry.Polygon([[39.05789266, 13.59051553],
                                [39.11335033, 13.59051553],
                                [39.11335033, 13.64477783],
                                [39.05789266, 13.64477783],
                                [39.05789266, 13.59051553]])


def mask_l8_sr(image):
    # Bits 3 and 5 are cloud shadow and cloud, respectively.
    cloud_shadow_bit_mask = (1 << 3)
    clouds_bit_mask = (1 << 5)
    # Get the pixel QA band.
    qa = image.select('pixel_qa')
    # Both flags should be set to zero, indicating clear conditions.
    mask = qa.bitwiseAnd(cloud_shadow_bit_mask).eq(0) and (qa.bitwiseAnd(clouds_bit_mask).eq(0))
    return image.updateMask(mask)


dataset = ee.ImageCollection(SATELLITE_SR).filterBounds(geometry).filterDatefilterDate('2016-01-01', '2016-12-31').map(
    mask_l8_sr).select(RGB)
image = dataset.reduce('median')
PERCENTILE_SCALE = 30  # Resolution in meters to compute the percentile at
percentiles = image.reduceRegion(ee.Reducer.percentile([0, 100], ['min', 'max']),
                                 geometry, PERCENTILE_SCALE).getInfo()
# Extracting the results is annoying because EE prepends the channel name
minVal = [val for key, val in percentiles.items() if 'min' in key]
# splitVal = next(val for key, val in percentiles.items() if 'split' in key)
maxVal = [val for key, val in percentiles.items() if 'max' in key]
print(minVal)
print(maxVal)
reduction = image.visualize(bands=RGB,
                            min=list(reversed(minVal)), # reverse since bands are given in the other way (b2,b3,4b)
                            max=list(reversed(maxVal)),
                            gamma=1)

# export reduction to drive

Answer 1

There is an error for python code:

dataset = ee.ImageCollection(SATELLITE_SR) \
    .filterBounds(geometry) \
    .filterDatefilterDate('2016-01-01', '2016-12-31') \
    .map(mask_l8_sr) \
    .select(RGB)

AttributeError: 'ImageCollection' object has no attribute 'filterDatefilterDate'

it should be:

dataset = ee.ImageCollection(SATELLITE_SR) \
    .filterBounds(geometry) \
    .filterDate('2016-01-01', '2016-12-31') \
    .map(mask_l8_sr) \
    .select(RGB)