SNAP vs. GDAL vs. Orfeo - reescalado de la imagen, corte del histograma

Question

Trabajo con una imagen Sentinel-2 .jp2 (banda roja, 10950 x 10950 píxeles). Mi objetivo es llegar al mismo resultado que SNAP con un script de Python. Ver el método y los parámetros de SNAP:

Así que esta es mi referencia (resultado con SNAP), quiero llegar a este resultado (representación en escala de grises QGIS, corte acumulativo - 2/98%):

Así que intenté replicarlo con GDAL:

import numpy as np
from osgeo import gdal, gdal_array

input = "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04.jp2"
output = "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_gdal.tif"

dataset = gdal.Open(input)
array = dataset.ReadAsArray()

percentile_025 = np.percentile(array, 2.5) # 349.0
percentile_975 = np.percentile(array, 97.5) # 3735.0

command = 'gdal_translate -scale ' + str(percentile_025) + ' ' + str(percentile_975)+ ' 0 255 -of GTiff -ot Byte' + ' ' + input + ' ' + output

os.system(command)

El resultado de GDAL no es el mismo, es un poco más brillante, las áreas blancas son más grandes. Los valores no son los mismos en el panel de capas (representación en escala de grises de QGIS, corte acumulativo - 2/98%):

El código Orfeo:

import otbApplication

Convert = otbApplication.Registry.CreateApplication("Convert")
Convert.SetParameterString("in", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04.jp2")
Convert.SetParameterString("out", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_hcut.tif")
Convert.SetParameterString("type","linear")
Convert.SetParameterString("hcp.high","2.5")
Convert.SetParameterString("hcp.low","2.5")
Convert.ExecuteAndWriteOutput()

Rescale = otbApplication.Registry.CreateApplication("Rescale")
Rescale.SetParameterString("in", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_hcut.tif")
Rescale.SetParameterString("out", "d:/bitbucket/cnn-lcm/T33TWM_A012703_20171127T100339_B04_orfeo.tif")
Rescale.SetParameterOutputImagePixelType("out", 1)
Rescale.SetParameterFloat("outmin", 0)
Rescale.SetParameterFloat("outmax", 255)
Rescale.ExecuteAndWriteOutput()

El resultado de Orfeo es muy similar al de GDAL (sólo hay 1-2 diferencias de valor en los píxeles). Y hay franjas grandes y problemáticas en el centro (representación en escala de grises de QGIS, corte acumulativo - 2/98%):

Así que finalmente mis preguntas:

¿Es posible eliminar las diferencias? ¿Es posible alcanzar exactamente el resultado de SNAP? ¿Y cómo?

Enlace de descarga de datos: http://sentinel-s2-l1c.s3.amazonaws.com/tiles/33/T/WM/2017/11/27/0/B04.jp2

Answer 1

He investigado esto y no he podido reproducir los resultados exactos que obtuviste de SNAP (específicamente no estoy seguro de por qué tu salida tiene un valor mínimo de 5). Sin embargo, pude confirmar que los resultados de SNAP y GDAL difieren significativamente. Dos ideas:

1. GDAL está implementando correctamente el clip del 95%. He examinado la trama original para identificar los valores mínimos/máximos que se han asignado a 0/255 (respectivamente) en la salida de GDAL. Resultaron ser 355 y 3729, que son diferentes de los umbrales de percentil 349/3735 que identificas arriba. Sin embargo, son precisamente los valores mínimos/máximos que redondo a 0/255 después de escalar con dichos umbrales percentiles ( ver aquí )
2. Puedes casi reproducir la salida de GDAL mediante SNAP si aumenta el parámetro "Precisión del histograma" al máximo (6) en la ventana "Estadísticas" antes de hacer la conversión. La salida SNAP/GDAL difiere ahora como máximo en 1, lo que se debe a las diferencias de redondeo.

No he investigado sobre Orfeo porque nunca lo he utilizado.