python+gdal对影像降位深

python+gdal对影像降位深

示例对影像从16位降为8位

对于几十个G的影像可参考作者分块处理影像文章中读取数据的方式

def compress(origin_file, output_file):
    """
    16位压缩成8位
    :param origin_file:
    :param output_file:
    :return:
    """
    img_array, projection, transform, rows, cols = read_tiff(origin_file)
    new_img_array = reduce_dimension(img_array)
    write_tiff(img_array=new_img_array, projection=projection, transform=transform, target_file=output_file)

def read_tiff(input_file):
    """
        读取影像
        :param input_file:输入影像
        :return: 影像数据，投影信息， 仿射变换参数
    """
    img = gdal.Open(input_file)  # 用完记得释放，是暂时存储在内存里的
    cols = width = img.RasterXSize
    rows = height = img.RasterYSize

    img_band = img.RasterCount
    img_array = img.ReadAsArray(0, 0, width, height)
    projection = img.GetProjection()
    transform = img.GetGeoTransform()
    img = None  # todo 释放内存，只有强制为None才可以释放干净
    del img
    gc.collect()
    return img_array, projection, transform, rows, cols

def reduce_dimension(img_array):
    """
    把RGB值转到0-255之间
    :param img_array:  tiff文件的数组
    :return: 一个RGB值在0-255之间的数组
    """
    new_img_array = np.empty(img_array.shape, img_array.dtype)
    # 此处为了减少循环中计算量，提前将图像高宽获取写死，不用在循环中再多次计算，提升性能
    band = len(img_array[:, 0, 0])
    height = len(img_array[0, :, 0])
    width = len(img_array[0, 0, :])
    for i in range(band):
        print('共%d个波段，正在处理第%d个波段，请稍候...' % (band, i + 1))
        for m in range(height):
            for n in range(width):
                new_img_array[i, m, n] = round(img_array[i, m, n] / 256)
    return new_img_array

def write_tiff(img_array, projection, transform, target_file):
    """
    将文件保存成8bit
    :param img_array: tiff文件的数组
    :param projection: 影像的投影
    :param transform: 影像的仿射变换
    :param target_file: 存储到的位置
    :return:
    """

    datatype = gdal.GDT_Byte

    if len(img_array.shape) == 3:
        img_bands, img_height, img_width = img_array.shape
    else:
        img_bands, (img_height, img_width) = 1, img_array.shape

    # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(target_file, img_width, img_height, img_bands, datatype)
    if dataset:
        dataset.SetGeoTransform(transform)  # 写入仿射变换参数
        dataset.SetProjection(projection)  # 写入投影
    for i in range(img_bands):
        dataset.GetRasterBand(i + 1).WriteArray(img_array[i])
    del dataset