Python读取16位的多光谱遥感图像

# 参考《python地理空间分析指南》
from osgeo import gdal_array
import operator
from functools import reduce
import matplotlib.pyplot as plt
import cv2

# 打开文件查看
src = 'phase2/20170218_UH_CASI_S4_NAD83.pix'  # 50波段遥感数据
# 来自DFC2018 Houston
# 下载https://hyperspectral.ee.uh.edu/?page_id=1075
pix_arr = gdal_array.LoadFile(src)
print(pix_arr.shape)  # 查看当前pix的形状信息
# 常量
NUMS = 65536  # 16位的图像
# 该数据的50个波段中RGB的索引
R = 21
G = 12
B = 7

(50, 1202, 4172)

# 显示其中三个波段合成的图像（RGB波段）
img_arr = cv2.merge((pix_arr[R]/float(NUMS), pix_arr[G]/float(NUMS), pix_arr[B]/float(NUMS)))  # 合并波段，缩小到0-1以显示
print(img_arr.shape)
plt.figure(figsize=(18.03, 62.58))
plt.imshow(img_arr)
plt.show()

(1202, 4172, 3)

# 计算直方图
def histogram(ima, bins=list(range(0, NUMS))):
    flat = ima.flat
    n = gdal_array.numpy.searchsorted(gdal_array.numpy.sort(flat), bins)
    n = gdal_array.numpy.concatenate([n, [len(flat)]])
    hist = n[1:]-n[:-1]
    return hist

# 直方图均衡化
def stretch(ima):
    hist = histogram(ima)
    lut = []
    for bt in range(0, len(hist), NUMS):
        # 步长尺寸
        step = reduce(operator.add, hist[bt:bt+NUMS]) / (NUMS-1)
        # 创建均衡的查找表
        n = 0
        for i in range(NUMS):
            lut.append(n / step)
            n += hist[i+bt]
	gdal_array.numpy.take(lut, ima, out=ima)
	return ima

# 原图直方图
x = [i for i in range(NUMS)]
hist_r = histogram(pix_arr[R])
hist_g = histogram(pix_arr[G])
hist_b = histogram(pix_arr[B])
plt.figure(figsize=(5, 5))
plt.plot(x, hist_b, color='b')
plt.plot(x, hist_g, color='g')
plt.plot(x, hist_r, color='r')
plt.show()

# 均衡化
stretched_r = stretch(pix_arr[R])
stretched_g = stretch(pix_arr[G])
stretched_b = stretch(pix_arr[B])
stretched_img = cv2.merge((stretched_r/float(NUMS), stretched_g/float(NUMS), stretched_b/float(NUMS)))
plt.figure(figsize=(18.03, 62.58))
plt.imshow(stretched_img)
plt.show()

# 均衡化后直方图
x = [i for i in range(NUMS)]
s_hist_r = histogram(stretched_r)
s_hist_g = histogram(stretched_g)
s_hist_b = histogram(stretched_b)
plt.figure(figsize=(5, 5))
plt.plot(x, s_hist_b, color='b')
plt.plot(x, s_hist_g, color='g')
plt.plot(x, s_hist_r, color='r')
plt.show()