Python 标准差拉伸（更新）

16bit转为8bit

2022.9.24 更新: 排除 0 值,代码如下.

# -*- coding: utf-8 -*-
from osgeo import gdal
import numpy as np
import os


def bit8(data, d_min, d_max, mean, std, Kn):
    uc_max = mean + Kn * std
    uc_min = mean - Kn * std
    k = (d_max - d_min) / (uc_max - uc_min)
    b = (uc_max * d_min - uc_min * d_max) / (uc_max - uc_min)
    if uc_min <= 0:
        uc_min = 0
    # 跳过 0 值 2022.9.24
    data = np.select(
        [data == 0, data == d_min, data <= uc_min, data >= uc_max, k * data + b < d_min, k * data + b > d_max,
         (k * data + b > d_min) & (k * data + b < d_max)],
        [0, d_min, d_min, d_max, d_min, d_max, k * data + b], data)
    return data


def convert2bit8(in_path):
    if not os.path.exists(in_path):
        print('输入路径不存在！')
        return ""
    for root, dirs, files in os.walk(in_path):
        for file in files:
            '''筛选tif文件'''
            if not file.split('.')[-1] == 'tif':
                continue
            '''输入路径的文件名'''
            file_name = os.path.join(root, file)
            '''创建输出路径文件名'''
            out_file = os.path.join(root, os.path.splitext(file)[0] + '_bit8.tif')
            '''文件存在则删除文件重新生成'''
            if os.path.exists(out_file):
                os.remove(out_file)
            in_tif = gdal.Open(file_name)
            width = in_tif.RasterXSize
            height = in_tif.RasterYSize

            '''跳过8bit'''
            if in_tif.ReadAsArray().dtype.name == 'uint8':
                continue

            geo_transform = in_tif.GetGeoTransform()
            print('左上角坐标 %f %f' % (geo_transform[0], geo_transform[3]))
            print('像素分辨率 %f %f' % (geo_transform[1], geo_transform[5]))
            '''True Color 1,2,3波段'''
            out_tif = gdal.GetDriverByName('GTiff').Create(out_file, width, height, in_tif.RasterCount, gdal.GDT_Byte)
            out_tif.SetProjection(in_tif.GetProjection())
            out_tif.SetGeoTransform(geo_transform)

            for i in range(1, int(in_tif.RasterCount) + 1):
                band = in_tif.GetRasterBand(i)
                data = band.ReadAsArray()
                # 排除 0 值 2022.9.24
                not0 = np.where(data > 0, data, np.nan)
                mean = np.nanmean(not0)
                std = np.nanstd(not0, ddof=1)
                # std = np.std(data)
                out_band = bit8(data, 0, 255, mean, std, 5)
                out_tif.GetRasterBand(i).WriteArray(out_band)

            out_tif.FlushCache()
            for i in range(1, int(in_tif.RasterCount) + 1):
                out_tif.GetRasterBand(i).ComputeStatistics(False)
            # out_tif.BuildOverviews('average', [2, 4, 8, 16, 32])
            del out_tif
    return out_file


if __name__ == '__main__':
    '''调用选择文件夹对话框获取输入路径'''
    path = r"F:\16bit\1"

    convert2bit8(path)
    print("Done!")