22

I work with a Sentinel-2 .jp2 image (red band, 10950 x 10950 pixels). My aim is to reach the same result what SNAP does with a Python script. See the SNAP method and parameters:

SNAP settings

So this is my reference (result with SNAP), I want to reach this result (QGIS grayscale representation, cumulative cut - 2/98%):

Reference in SNAP

So I tried to replicate it with 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)

The GDAL result is not the same, its a bit brighter, the white areas are bigger. The values are not the same on the layers panel (QGIS grayscale representation, cumulative cut - 2/98%):

Band in GDAL

The Orfeo code:

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()

The Orfeo result is very similar to GDAL (only 1-2 value differences in pixels). And there are big, problematic strips in the middle (QGIS grayscale representation, cumulative cut - 2/98%):

Image with Orfeo

So finally my questions:

Is it possible to eliminate the differences? Is it possible to reach exactly the result of SNAP? And how?

Download link to data: http://sentinel-s2-l1c.s3.amazonaws.com/tiles/33/T/WM/2017/11/27/0/B04.jp2

  • 3
    If the result from the SNAP method "between 95% clipped histogram" and from your "numpy percentile 2.5 - 97.5" method are different then it probably means that SNAP is doing it in another way, maybe based on average and standard deviation. Perhaps it is possible to find it from here github.com/senbox-org/s2tbx. – user30184 Jan 4 '18 at 15:37
  • Yes, I though that too, but I am not sure about it, since I didn't find the corresponding code. – pnz Jan 4 '18 at 16:03
  • This gis.stackexchange.com/questions/172721/… talks about statistics as well. – user30184 Jan 4 '18 at 17:21
  • I would move this question to forum.step.esa.int to find out what's happening exactly when doing your process in SNAP – GCGM Jan 10 '18 at 8:27
  • any nodata inside your raster? percentile_025 = np.percentile(array[array!=nodata], 2.5), etc might give a different value? – user1269942 Nov 12 '18 at 4:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.