所以这是我的参考（使用SNAP的结果），我想达到此结果（QGIS灰度表示，累积切割-2 / 98％）：
所以我尝试使用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)
GDAL结果是不一样的，它的亮度更高，白色区域更大。图层面板上的值不同（QGIS灰度表示，累积切割-2/98％）：
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()
Orfeo结果与GDAL非常相似（像素差异只有1-2个值）。中间有大的有问题的条带（QGIS灰度表示，累积剪切-2/98％）：
所以最后我的问题：
是否可以消除差异？是否有可能完全达到SNAP的结果？以及如何操作？
下载数据链接：
http://sentinel-s2-l1c.s3.amazonaws.com/tiles/33/T/WM/2017/11/27/ 0 / B04.jp2
#1 楼
我调查了一下，但无法重现从SNAP获得的准确结果（特别是我不确定为什么输出的最小值为5）。但是，我能够确认SNAP和GDAL的结果存在显着差异。有两种见解：GDAL正在正确实施95％剪辑。我检查了原始栅格，以确定在GDAL输出中分别映射到0/255的最小/最大值。结果是355和3729，这与您在上面确定的349/3735百分位数阈值不同。但是，它们恰好是在按所述百分位数阈值进行缩放后舍入为0/255的最小值/最大值（请参见此处）
如果将“直方图精度”参数增加到，则几乎可以使用SNAP再现GDAL输出。进行转换之前，请在“统计信息”窗口中输入最大值（6）。 SNAP / GDAL输出现在最多相差1，这是由于舍入的差异所致。
我从未使用过Orfeo，因为我从未使用过它。
评论
似乎有一个有趣的答案，如果它获得5票赞成（加1票）的支持，我会接受。由于似乎没有明显的问题答案。在对帖子发表一些反馈之后，我想我应该对其进行更新：-percentile_025 = np.percentile（array [array！= nodata]，2.5）-计算栅格差异；不要在GDAL百分位数差异的背景下出现什么情况（349，355 ..）
– pnz
20-10-16在8:30
349 vs 355仅仅是因为输出数据类型是整数，因此缩放后的值必须四舍五入。在我上面链接的帖子中查找最终公式。由于我们从[349.3735]缩放到[0,255]，因此我们插入min = 349，max = 3735，a = 0和b =255。然后，f（349）= 0，如预期的那样。 f（350）= 0.075 ...，但是最终将其舍入为0。f（355）= 0.451 ...是x的最大值，它将被舍入为零...之后f（356）= 0.521 ...向上舍入为1。在范围的高端（3735和3729）也会发生相同的情况。希望能澄清
–尼克
20-10-16在16:45
评论
如果SNAP方法“在95％裁剪的直方图之间”的结果与您的“ numpy百分位数2.5-97.5”方法的结果不同，则可能意味着SNAP正在另一种方式进行，可能基于平均值和标准偏差。也许可以从这里github.com/senbox-org/s2tbx找到它。是的，我也这样做，但是我不确定，因为我没有找到相应的代码。
此gis.stackexchange.com/questions/172721/…还讨论了统计信息。
栅格中是否有任何数据？ percentile_025 = np.percentile（array [array！= nodata]，2.5）等可能会给出不同的值？
即使您在QGIS中使用了相似的可视化参数（95％的裁切率），您处理的图像（GDAL和Orfeo）也从0缩放到255，而SNAP图像从5缩放到255。 SNAP和其他图像（使用栅格计算器），以及该差异图像的直方图？这样，我们可以查看实际图像值，该图像值比可视化图像更易于解释