gdal读取tif文件

gdal读取tif文件

import gdal
import numpy as np


def stretch_16to8(bands, lower_percent=2, higher_percent=98):
    out = np.zeros_like(bands, dtype=np.uint8)
    n = bands.shape[0]
    for i in range(n):
        a = 0  # np.min(band)
        b = 255  # np.max(band)
        c = np.percentile(bands[i, :, :], lower_percent)
        d = np.percentile(bands[i, :, :], higher_percent)
        t = a + (bands[i, :, :] - c) * (b - a) / (d - c)
        t[t < a] = a
        t[t > b] = b
        out[i, :, :] = t
    return out.astype(np.uint8)


def main():

    img_path = "xxxx.tif"
    save_path = 'save.tif'

    dataset = gdal.Open(img_path)
    if dataset == None:
        print(img_path + "文件无法打开")
        return
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_bands = dataset.RasterCount  # 波段数
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 获取数据
    im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
    im_proj = dataset.GetProjection()  # 获取投影信息

    # 16位转8位
    image = stretch_16to8(im_data)

    # =============================
    # 测试
    import cv2
    image_ = np.transpose(image, (1, 2, 0))[:, :, 0:3]
    cv2.imwrite('image.png', image_)
    # =============================

    # 数据保存
    driver = gdal.GetDriverByName("GTiff")  # 创建文件
    dataset = driver.Create(save_path, im_width, im_height, im_bands, gdal.GDT_Byte)
    if (dataset != None):
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i + 1).WriteArray(image[i])
    del dataset


if __name__ == "__main__":
    main()