面临的问题就是我们在进行遥感影像处理的时候，一般影像为16位，理论上范围是0-65535，但是实际上有效值一半位0-1000，1000以后的值几乎没有，直接除以65535进行归一化是不合适的，尤其是在深度学习中。
进行标准化：
（1）错误的方式： x/32767.5

遥感影像数据增强

1 了解遥感影像
- 1.1 高分影像信息统计
- 1.2 了解北京2号影像
- 1.3 统计上面信息的代码
2 数据增强
- 2.1 普通方式（超级耗时，没法用）
- 2.2 改进方式（用np.percentile（）函数）
- 2.3 有黑边占30%的把黑边排除掉
参考：

1 了解遥感影像

1.1 高分影像信息统计

这个高分影像，虽然值范围理论上是0-65535，但是99%的值基本上都小于500。

1.2 了解北京2号影像

（1）TRIPLESAT_3_PSH_L4_20190811025844_0023CBVI_026.pix
这一景是白天拍的，比较亮,它的值有一部分大于1000，比如band1,蓝光波段，40%多的值小于1000。而下面比较暗的图像，80%多的值小于1000，但是他们几乎不超过5000。

（2）TRIPLESAT_3_PSH_L4_20190815030718_0023E3VI_025.pix
这一景可能是黄昏拍的，比较暗

1.3 统计上面信息的代码

import gdal
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from collections import Counter#   一、 读取数据
in_ds = gdal.Open(r"L:\2之前的影像数据\2北京2号数据\2忠县北京2号数据15景\TRIPLESAT_3_PSH_L4_20190815030718_0023E3VI_025.pix")
#  二、 参看数据信息，获取numpy值
print("open tif file succeed")
height = in_ds.RasterYSize  # 获取数据高度
print("高_行：",height)
width = in_ds.RasterXSize  # 获取数据宽度
print("宽_列：",width)
outbandsize = in_ds.RasterCount  # 获取数据波段数
print("波段数：",outbandsize)
im_geotrans = in_ds.GetGeoTransform()  # 获取仿射矩阵信息
print("仿射信息：",im_geotrans)
im_proj = in_ds.GetProjection()  # 获取投影信息
print("投影信息：",im_proj)
datatype = in_ds.GetRasterBand(1).DataType    # 这个数据类型没搞清楚，1，2，3代表什么https://blog.csdn.net/u010579736/article/details/84594742
print("数据类型：",datatype)for i in range(outbandsize):print("band%d"%(i+1))band = in_ds.GetRasterBand(i+1)  # 注意这里从1开始不是从0开始band_data = band.ReadAsArray()  # 获取数据print("数组类型", band_data.dtype)print("数组类型名", band_data.dtype.name)print("数组形状", band_data.shape)print("全部波段最大值：", np.max(band_data))print("全部波段最小值", np.min(band_data))print("全部波段唯一值", np.unique(band_data))num500 = np.sum(band_data < 500)print("小于500的比例：",num500/(band_data.shape[0]*band_data.shape[1]))num1000 = np.sum(band_data < 1000)print("小于1000的比例：", num1000 / (band_data.shape[0] * band_data.shape[1]))num5000 = np.sum(band_data < 5000)print("小于5000的比例：", num5000 / (band_data.shape[0] * band_data.shape[1]))print('==========')

2 数据增强

不同于一般的8位数据，数据库中的遥感图像都是16位的无符号整形数据，而且大部分集中在[0,500]的范围内，对于8位的RGB图像来说根本无法显示正确的图像。因此，必须对灰度级进行压缩，从65536个灰度级压缩到256个灰度级范围内。
但是，仅仅使用线性压缩是无法解决问题的，因为前面也说过，数据主要集中在[0,500]范围内，即它只占用了前面的几百个灰度级，后面的根本没用到。为此，我们必须丢掉出现概率小的灰度级，保留出现概率大的灰度级范围。

2.1 普通方式（超级耗时，没法用）

自己用下面代码试了一下，时间太久，就及时停止了，最少得4个小时


# 将16位遥感影像转换为8位，并进行百分比截断
# https://blog.csdn.net/qq_43814272/article/details/109741119?utm_source=app&app_version=4.17.2&code=app_1562916241&uLinkId=usr1mkqgl919blen
import gdal
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from collections import Counter
import timedef read_img(input_file):"""读取影像Inputs:input_file:16位图像数据的路径Outputs:array_data：保存图像数据的numpy.array数组rows:高度cols:宽度bands:深度"""in_ds = gdal.Open(input_file)rows = in_ds.RasterYSize  # 获取数据高度cols = in_ds.RasterXSize  # 获取数据宽度bands = in_ds.RasterCount  # 获取数据波段数datatype = in_ds.GetRasterBand(1).DataType  # 这个数据类型没搞清楚，1，2，3代表什么https://blog.csdn.net/u010579736/article/details/84594742#print("数据类型：", datatype)array_data = in_ds.ReadAsArray()  # 将数据写成数组，读取全部数据，结果是numpy数组# print("数组类型名", img_array.dtype.name)    # uint16# print("数组形状", array_data.shape)   #  (4, 36786, 37239) ,波段，行，列return array_data, rows, cols, bandsdef cumulativehistogram(array_data, band_min, band_max,low_percent,high_percent):"""累计直方图统计,得到单波段图像的灰度直方图Inputs:array_data:16位图像数据band_min:最小像素值band_max:最大像素值Outputs：cutmax:累积分布函数（CDF）==0.98对应的灰度级cutmin:累积分布函数（CDF）==0.02对应的灰度级"""# 逐波段统计最值print("2统计数组形状", array_data.shape)rows=array_data.shape[0]cols=array_data.shape[1]gray_level = int(band_max - band_min + 1)gray_array = np.zeros(gray_level)for row in range(rows):for col in range(cols):gray_array[int(array_data[row, col] - band_min)] += 1count_percent2 = rows*cols * low_percentcount_percent98 = rows*cols* high_percentcutmax = 0cutmin = 0for i in range(1, gray_level):gray_array[i] += gray_array[i - 1]# 下面这种可能出现gray_array[0]>= count_percent2，这样cutmin = 0，本身酒等于0，所以就不需要另外写elseif (gray_array[i] >= count_percent2 and gray_array[i - 1] <= count_percent2):cutmin = i + band_minif (gray_array[i] >= count_percent98 and gray_array[i - 1] <= count_percent98):cutmax = i + band_minreturn cutmin, cutmaxdef write_img(read_path, img_array):"""read_path:原始文件路径img_array：numpy数组"""read_pre_dataset = gdal.Open(read_path)img_transf = read_pre_dataset.GetGeoTransform()  # 仿射矩阵img_proj = read_pre_dataset.GetProjection()  # 地图投影信息print("1读取数组形状", img_array.shape)if 'uint8' in img_array.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in img_array.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img_array.shape) == 3:img_bands, im_height, im_width = img_array.shapeelse:img_bands, (im_height, im_width) = 1, img_array.shapefilename = read_path[:-4] + '_unit8' + ".tif"driver = gdal.GetDriverByName("GTiff")  # 创建文件驱动# 注意这里先写宽再写高，对应数组的高和宽，这里不对应才对# https://vimsky.com/examples/detail/python-method-gdal.GetDriverByName.htmldataset = driver.Create(filename, im_width, im_height, img_bands, datatype)dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数dataset.SetProjection(img_proj)  # 写入投影# 写入影像数据if img_bands == 1:dataset.GetRasterBand(1).WriteArray(img_array)else:for i in range(img_bands):dataset.GetRasterBand(i + 1).WriteArray(img_array[i])def compress(origin_16,low_per=0.004,high_per=0.996):"""根据任务不同选择的截断点不同，一般不要大于0.02，和小于0.98Input:origin_16:16位图像路径low_per=0.004百分比截断的低点high_per=0.096百分比截断的高点Output:output:8位图像路径"""array_data, rows, cols, bands = read_img(origin_16) # array_data, (4, 36786, 37239) ,波段，行，列print("1读取数组形状", array_data.shape)# 这里控制要输出的是几位compress_data = np.zeros((bands,rows, cols),dtype="uint8")for i in range(bands):band_min = np.min(array_data[i, :, :])band_max = np.max(array_data[i, :, :])cutmin, cutmax = cumulativehistogram(array_data[i,:, :], band_min, band_max,low_percent=low_per,high_percent=high_per)compress_scale = (cutmax - cutmin) / 255for j in range(rows):for k in range(cols):if (array_data[i,j, k] < cutmin):array_data[i,j, k] = cutminif (array_data[i,j, k] > cutmax):array_data[i,j, k] = cutmaxcompress_data[i,j, k] = round((array_data[i,j, k] - cutmin) / compress_scale)write_img(origin_16, compress_data)if __name__ == '__main__':start = time.time()ds = r"I:\5GF2\0rawData\HubeiLuotian5\GF2_PMS1_E115.4_N30.6_20181004_L1A0003497911-PAN1_ORTHO_PSH.img"compress(ds)print("time_cost:",time.time()-start)

2.2 改进方式（用np.percentile（）函数）

亲自测试下面的代码，耗时半个小时，主要是写数据花了有20分钟时间。数据增强只用了15分钟

# 将16位遥感影像转换为8位，并进行百分比截断
# https://blog.csdn.net/qq_43814272/article/details/109741119?utm_source=app&app_version=4.17.2&code=app_1562916241&uLinkId=usr1mkqgl919blen
import gdal
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from collections import Counter
import timedef read_img(input_file):"""读取影像Inputs:input_file:16位图像数据的路径Outputs:array_data：保存图像数据的numpy.array数组rows:高度cols:宽度bands:深度"""in_ds = gdal.Open(input_file)rows = in_ds.RasterYSize  # 获取数据高度cols = in_ds.RasterXSize  # 获取数据宽度bands = in_ds.RasterCount  # 获取数据波段数# 这个数据类型没搞清楚，1，2，3代表什么https://blog.csdn.net/u010579736/article/details/84594742# GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6datatype = in_ds.GetRasterBand(1).DataTypeprint("数据类型：", datatype)array_data = in_ds.ReadAsArray()  # 将数据写成数组，读取全部数据，numpy数组,array_data.shape  (4, 36786, 37239) ,波段，行，列del in_dsreturn array_data, rows, cols, bandsdef write_img(read_path, img_array):"""read_path:原始文件路径img_array：numpy数组"""read_pre_dataset = gdal.Open(read_path)img_transf = read_pre_dataset.GetGeoTransform()  # 仿射矩阵img_proj = read_pre_dataset.GetProjection()  # 地图投影信息print("1读取数组形状", img_array.shape,img_array.dtype.name)# GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,if 'uint8' in img_array.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in img_array.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img_array.shape) == 3:img_bands, im_height, im_width = img_array.shapeelse:img_bands, (im_height, im_width) = 1, img_array.shapefilename = read_path[:-4] + '_unit8' + ".tif"driver = gdal.GetDriverByName("GTiff")  # 创建文件驱动# 注意这里先写宽再写高，对应数组的高和宽，这里不对应才对# https://vimsky.com/examples/detail/python-method-gdal.GetDriverByName.htmldataset = driver.Create(filename, im_width, im_height, img_bands, datatype)dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数dataset.SetProjection(img_proj)  # 写入投影# 写入影像数据if img_bands == 1:dataset.GetRasterBand(1).WriteArray(img_array)else:for i in range(img_bands):dataset.GetRasterBand(i + 1).WriteArray(img_array[i])def compress(origin_16,low_per=0.4,high_per=99.6):"""Input:origin_16:16位图像路径low_per=0.4   0.4%分位数，百分比截断的低点high_per=99.6  99.6%分位数，百分比截断的高点Output:output:8位图像路径"""array_data, rows, cols, bands = read_img(origin_16) # array_data, (4, 36786, 37239) ,波段，行，列print("1读取数组形状", array_data.shape)# 这里控制要输出的是几位compress_data = np.zeros((bands,rows, cols),dtype="uint8")for i in range(bands):# 得到百分比对应的值cutmin = np.percentile(array_data[i, :, :], low_per)cutmax = np.percentile(array_data[i, :, :], high_per)data_band = array_data[i]# 进行截断,用data_band = np.clip(data_band,cutmin,cutmax)也可以data_band[data_band<cutmin] = cutmindata_band[data_band>cutmax] = cutmax# 进行缩放compress_data[i,:, :] = np.around( (data_band[:,:] - cutmin) *255/(cutmax - cutmin) )print("最大最小值：",np.max(compress_data),np.min(compress_data))write_img(origin_16, compress_data)if __name__ == '__main__':start = time.time()ds = r"I:\5GF2\0rawData\HubeiLuotian5\GF2_PMS1_E115.4_N30.6_20181004_L1A0003497911-PAN1_ORTHO_PSH.img"compress(ds)print("time_cost:",time.time()-start)

输出：

目前的效率仍然还有待进一步提高，有想法的小伙伴可以评论区交流

2.3 有黑边占30%的把黑边排除掉

黑边的值一般是0，百分比截断的时候忽略0的影响

# 将16位遥感影像转换为8位，并进行百分比截断
# https://blog.csdn.net/qq_43814272/article/details/109741119?utm_source=app&app_version=4.17.2&code=app_1562916241&uLinkId=usr1mkqgl919blen
import gdal
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from collections import Counter
import time
import osdef read_img(input_file):"""读取影像Inputs:input_file:16位图像数据的路径Outputs:array_data：保存图像数据的numpy.array数组rows:高度cols:宽度bands:深度"""in_ds = gdal.Open(input_file)rows = in_ds.RasterYSize  # 获取数据高度cols = in_ds.RasterXSize  # 获取数据宽度bands = in_ds.RasterCount  # 获取数据波段数# 这个数据类型没搞清楚，1，2，3代表什么https://blog.csdn.net/u010579736/article/details/84594742# GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6datatype = in_ds.GetRasterBand(1).DataTypeprint("数据类型：", datatype)array_data = in_ds.ReadAsArray()  # 将数据写成数组，读取全部数据，numpy数组,array_data.shape  (4, 36786, 37239) ,波段，行，列del in_dsreturn array_data, rows, cols, bandsdef write_img(read_path, img_array):"""read_path:原始文件路径img_array：numpy数组"""read_pre_dataset = gdal.Open(read_path)img_transf = read_pre_dataset.GetGeoTransform()  # 仿射矩阵img_proj = read_pre_dataset.GetProjection()  # 地图投影信息print("1readshape:", img_array.shape,img_array.dtype.name)# GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,if 'uint8' in img_array.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in img_array.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img_array.shape) == 3:img_bands, im_height, im_width = img_array.shapeelse:img_bands, (im_height, im_width) = 1, img_array.shapefilename = read_path[:-4] + '_unit8' + ".tif"driver = gdal.GetDriverByName("GTiff")  # 创建文件驱动# 注意这里先写宽再写高，对应数组的高和宽，这里不对应才对# https://vimsky.com/examples/detail/python-method-gdal.GetDriverByName.htmldataset = driver.Create(filename, im_width, im_height, img_bands, datatype)dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数dataset.SetProjection(img_proj)  # 写入投影# 写入影像数据if img_bands == 1:dataset.GetRasterBand(1).WriteArray(img_array)else:for i in range(img_bands):dataset.GetRasterBand(i + 1).WriteArray(img_array[i])def compress(origin_16,low_per_raw=0.01,high_per_raw=0.96):"""Input:origin_16:16位图像路径low_per=0.4   0.4%分位数，百分比截断的低点high_per=99.6  99.6%分位数，百分比截断的高点Output:output:8位图像路径"""array_data, rows, cols, bands = read_img(origin_16) # array_data, (4, 36786, 37239) ,波段，行，列print("1shape:", array_data.shape)# 这里控制要输出的是几位compress_data = np.zeros((bands,rows, cols),dtype="uint8")for i in range(bands):# 得到百分比对应的值，得到0代表的黑边的占比cnt_array = np.where(array_data[i, :, :], 0, 1)num0 = np.sum(cnt_array)kk = num0 / (rows * cols)    # 得到0的比例# 这里是去掉黑边0值，否则和arcgis的不一样，这样也有偏差，只算剩下的百分比low_per = low_per_raw + kk - low_per_raw * kk  # (A*x-A*KK)/(A-A*KK)=0.01, X = 0.99KK+0.01low_per = low_per * 100high_per = (1 - high_per_raw) * (1 - kk)  # A*x/(A-A*KK) = 0.04, X =  0.04-(0.04*kk)high_per = 100 - high_per * 100print("baifen:", low_per, high_per)cutmin = np.percentile(array_data[i, :, :], low_per)cutmax = np.percentile(array_data[i, :, :], high_per)print("duandian:",cutmin,cutmax)data_band = array_data[i]# 进行截断data_band[data_band<cutmin] = cutmindata_band[data_band>cutmax] = cutmax# 进行缩放compress_data[i,:, :] = np.around( (data_band[:,:] - cutmin) *255/(cutmax - cutmin) )print("maxANDmin：",np.max(compress_data),np.min(compress_data))# 下面可以直接将处理后的遥感图像输出，在神经网络预测的时候可以删掉这一句write_img(origin_16, compress_data)# return compress_data# 1 获取文件的函数
def listdir(path, list_name):"""获取文件夹下的文件，前一个是文件夹，后一个是要保存的列表"""panduan_geshi = [".tif", ".pix", ".img",]for file in os.listdir(path):file_path = os.path.join(path, file)if os.path.isdir(file_path):listdir(file_path, list_name)elif os.path.isfile(file_path):if os.path.splitext(file_path)[1] in panduan_geshi:list_name.append(file_path)else:print("false"*1000)if __name__ == '__main__':lissst = []pathh = r"I:\5GF2\1ALLimg20211130\0rawData"listdir(pathh , lissst)print(len(lissst))#print(lissst)for file_name in lissst:print(file_name)start = time.time()ds = file_namecompress(ds)print("time_cost:",time.time()-start)

参考：

看了很多博客，都和下面这两个一样，数据位数转换看本教程基本就够用了。
上面的案例结合下面的两个函数实现

【1】超级好(11条消息) 遥感影像16位转8位(python)_qq_43814272的博客-CSDN博客
https://blog.csdn.net/qq_43814272/article/details/109741119?utm_source=app&app_version=4.17.2&code=app_1562916241&uLinkId=usr1mkqgl919blen

【2】(11条消息) 遥感影像32位转8位（python）_海里的羊的博客-CSDN博客
https://blog.csdn.net/weixin_42990464/article/details/107152799?utm_source=app&app_version=4.17.2&code=app_1562916241&uLinkId=usr1mkqgl919blen

7. gdal进行遥感影像的16位转8位和百分比截断增强（看这篇就够了）相关推荐

python/gdal处理遥感影像（读取、投影转换、裁剪、建立图像金字塔等）
python/gdal处理遥感影像(读取.投影转换.裁剪.建立图像金字塔等) gdal库简单介绍 python使用gdal 一.安装python环境二.安装gdal库三.使用gdal处理遥感影像 ...
gdal进行遥感影像读写_如何使用遥感影像进行矿物勘探
gdal进行遥感影像读写 Meet Jose Manuel Lattus, a geologist from Chile. In the latest Soar Cast, he discusses ...
8影像计算ndvi landsat_使用GDAL读取遥感影像的信息
读取影像数据集的元数据 GDAL已经提供了足够方便的函数,可以读取影像的一些元数据信息, 从而方便对数据进行处理.GDAL一般是以字典的形式对元数据进行组织的, 但是对于不同的栅格数据类型,元数据的类 ...
原创程序|基于GDAL的遥感影像批量处理工具介绍（三）
本文主要介绍基于C#+GDAL-Python实用工具开发的遥感影像批量处理工具,该工具目前主要包括影像批量切片生成KML文件和影像批量切片生成Tiles文件.该工具.Net框架版本为4.0,GDAL版 ...
Python地学分析 — GDAL对遥感影像重投影
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
基于c语言的遥感图像处理,基于GDAL的遥感影像显示（C#版）
在matlab中实现遥感影像和shp文件的结合显示 clc;close all;clear; road=shaperead('boston_roads.shp'); %读取shape文件 figure ...
位运算全面总结，关于位运算看这篇就够了
文章目录 1 位运算概述 2 位运算的性质 2.1 运算符的优先级 2.2 位运算符的运算律 3 位运算高级操作 4 负数的位运算 5 位运算的一些应用 6 位运算例题 6.1 更新二进制位 6.2 ...
影像科dsa为什么必须买维修保险_了解什么是DSA，看这篇就够了
相对于超声.CT.磁共振等,DSA可能对大部分人而言都是比较陌生的,目前国内DSA市场基本都是进口品牌主导,国内能生产这种医疗设备的屈指可数,也就是说,目前我们在医院看到的DSA设备基本都来自&quo ...
深度学习笔记-遥感影像转为tensor前的检查及线性拉伸
@TOC 前言-为什么要做线性拉伸近期在开展语义分割任务,用到的数据是经过SNAP预处理的Sentinel-1的SAR数据.工作需要将相关的Image和label送入模型,进行目标地物的分割. 这里 ...
遥感影像16位转8位
现在常用卫星影像基本上都是16位影像,如GF,ZY3,Landsat8,WV等,有时我们需要将16位影像降到8位影像,这样不仅减小了数据量,也便于后期处理. 通常的软件在处理降位时会存在一些问题,如曝 ...

7. gdal进行遥感影像的16位转8位和百分比截断增强（看这篇就够了）