28

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

8
  • 5
    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.
    – poetyi
    Jan 4 '18 at 16:03
  • This gis.stackexchange.com/questions/172721/… talks about statistics as well.
    – user30184
    Jan 4 '18 at 17:21
  • 2
    any nodata inside your raster? percentile_025 = np.percentile(array[array!=nodata], 2.5), etc might give a different value? Nov 12 '18 at 4:21
  • 1
    Even though you use similar visualisation parameters in QGIS (95% cum. cut), the images you processed (GDAL and Orfeo) are scaled from 0 to 255, and the SNAP image is scaled from 5 to 255. Could you calculate the difference between the SNAP and other images (with raster calculator), together with a histogram of that difference image? This way we can take a look at the actual image values, which are IMHO easier to interpret than the visualization
    – danscr
    Jan 8 '20 at 11:42
2

I looked into this and was not able to reproduce the exact results you got from SNAP (specifically I'm not sure why your output has a minimum value of 5). However, I was able to confirm that the results from SNAP and GDAL differ significantly. Two insights:

  1. GDAL is implementing the 95% clip correctly. I examined the original raster to identify the minimum/maximum values that were mapped to 0/255 (respectively) in the GDAL output. These turned out to be 355 and 3729, which are different from the 349/3735 percentile thresholds you identify above. However, they are precisely the minimum/maximum values that round to 0/255 after scaling with said percentile thresholds (see here)
  2. You can almost reproduce the GDAL output using SNAP if you increase the 'Histogram accuracy' parameter to the maximum (6) in the 'Statistics' window before doing the conversion. The SNAP/GDAL output now differ by at most 1, which is due to differences in rounding.

I did not look into Orfeo as I have never used it.

2
  • Seems an interesting answer, if it reach a like 5 up votes (added one), I will accept it. Since it seems there is no obvious answer for the question. After some feedback on the post, I guess I should update the it: - percentile_025 = np.percentile(array[array!=nodata], 2.5) - Calculate raster differences; Donno what can be in the background of the GDAL percentile differences (349, 355..)
    – poetyi
    Oct 16 '20 at 8:30
  • 349 vs 355 is simply due to the fact the output datatype is integer and therefore the scaled values must be rounded. Look a the final formula in the post I linked to above. Since we are scaling from [349.3735] to [0,255] we plug in min = 349, max = 3735, a = 0, and b = 255. Then, f(349) = 0, as expected. f(350) =0.075... , but that winds up getting rounded to 0. f(355) = 0.451... is the highest value of x that will get rounded to zero ... after that f(356) = 0.521... rounds up to 1. Same thing happens at the high-end of the range (3735 vs 3729). Hope that clarifies
    – Nick
    Oct 16 '20 at 16:45

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.