python GDAL遥感影像16位转8位——gdal_translate实现

可以使用-scale重新拉伸像素值。例如，此选项可用于将16位（0至65535）数据缩放为8位（0到255）。选项采用4个参数：src_min src_max dst_min dst_max。前两个描述输入图像的最小值和最大值，后一个描述输出栅格的相同值。如果省略了dst_min和dst_max，则值将重新缩放为8位等效值（0到255）。如果src_min和src_max也被省略，该函数将自动计算输入栅格的值范围。

文档中提到，通过使用-scale_bn（其中bn是波段号），可以多次使用该选项来应用于特定波段。

命令行实现：

gdal_translate -of GTiff  \
-ot Byte  \
-b 3 -b 2 -b 1 \
-scale_1 0 1068 0 255 \
-scale_2 0 1459 0 255 \
-scale_3 0 1228 0 255 \
-co BIGTIFF=YES \
uint16.tif uint8.tif

Python代码实现：

def get_p2_p98(ch, total):
    """
    计算2%，98%
    :param ch: 累计直方图
    :param total: 像素个数
    :return:
    """
    p2 = total * 0.02
    p98 = total * 0.98
    min2 = 0
    max98 = 65535
    for index, v in enumerate(ch):
        if v > p2:
            min2 = index - 1
            break
    for index, v in enumerate(ch):
        if v > p98:
            max98 = index
            break

    return min2, max98


dst_tif_int8 = "int8.tif"
# 新建一个tif
ds: gdal.Dataset = gdal.Open(dst_tif_int16)
# 计算直方图
p2_98_list = []
for i in range(1, 4):
    band: gdal.Band = ds.GetRasterBand(i)
    # 1. 计算直方图
    h = band.GetHistogram(-0.5, 65535.5, 65535, 0, 1)
    # 像素个数
    total = sum(h)
    # 2. 计算累计直方图
    ch = np.cumsum(h)
    # 3. 计算2%，98%
    p2, p98 = get_p2_p98(ch, total)
    p2_98_list.append({"p2": p2, "p98": p98})
print("p2_98_list:", p2_98_list)
# 使用-scale_1 -scale_2 -scale_3 拉伸每个波段, -b指定波段顺序
options = f"-of GTiff " \
          f"-ot Byte " \
          f"-b 3 -b 2 -b 1 " \
          f"-scale_1 {p2_98_list[2]['p2']} {p2_98_list[2]['p98']} 0 255 " \
          f"-scale_2 {p2_98_list[1]['p2']} {p2_98_list[1]['p98']} 0 255 " \
          f"-scale_3 {p2_98_list[0]['p2']} {p2_98_list[0]['p98']} 0 255 " \
          f"-co BIGTIFF=YES "
translate_options = gdal.TranslateOptions(options=options)
gdal.Translate(dst_tif_3857_int8, dst_tif_3857, options=translate_options)