利用python同时显示不同尺度遥感影像的直方图

利用python同时显示不同尺度遥感影像的直方图，以观察不同分辨率图像分布的变换情况。

我想验证的结论是，不同尺度的遥感影像的分布不同。

from osgeo import gdal
from osgeo import gdal_array
import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv

def plothist(image):
    max = np.max(image)
    min = np.min(image)
    length = max - min
    print(max, "", min)
    hist = []

    for value in range(0, 255):
        times = np.sum(image == value)
        hist.append(times)
    return hist

def normalize(image):
    max = np.max(image)
    min = np.min(image)
    diff = max - min
    image = ((image - min)/diff) * 255
    image = image.astype(np.int)
    return image

filepath = "LC81200342017074LGN00\LC81200342017074LGN00_B2.TIF"
ds = gdal.Open(filepath)
datatype = ds.GetRasterBand(1).DataType
image = np.zeros((ds.RasterYSize, ds.RasterXSize, 1), 
    dtype = gdal_array.GDALTypeCodeToNumericTypeCode(datatype))

image_data = ds.GetRasterBand(1)
image[:, :, 0] = image_data.ReadAsArray()
experiment_image = image [2000:3000, 2000:3000, 0].reshape(1000, 1000)
max = np.max(experiment_image)
min = np.min(experiment_image)

norm_image = normalize(experiment_image)
norm_image = norm_image.astype(np.uint8)
dst_2 = cv.resize(norm_image, (500, 500), interpolation = cv.INTER_LINEAR)
dst_4 = cv.resize(norm_image, (250, 250), interpolation = cv.INTER_LINEAR)
dst_8 = cv.resize(norm_image, (125, 125), interpolation = cv.INTER_LINEAR)
hist_1 = plothist(norm_image)
hist_2 = plothist(dst_2)
hist_4 = plothist(dst_4)
hist_8 = plothist(dst_8)

for i in [[hist_1, 'red'], [hist_2, 'green'], [hist_4, 'blue'], [hist_8, 'black']]:
    plt.plot(range(255), i[0], color = i[1])

好，果然分布不同。

上面的结论已经认为是错误的。

下面做补充。上面的结果和结论是有问题的，因为我没考虑像元个数导致下面的直方图出现了压缩的情况。所以得出的不同尺度的遥感影像的分布不同的结论。今天对结果进行修改。

首先对直方图的个数要进行归一化，使之不造成压缩直方图形状的问题。

def pixelnumbernorm(hist, ratio):
    for i in range(256):
        hist[i] = int((hist[i])/(ratio*ratio))
    return hist

下面做显示：

hist_1 = pixelnumbernorm(hist_1, 8)
hist_2 = pixelnumbernorm(hist_2, 4)
hist_4 = pixelnumbernorm(hist_4, 2)
hist_8 = pixelnumbernorm(hist_8, 1)

所以可以看到，不同尺度遥感影像的直方图是类似的，所以结果就是尺度对于遥感影像的分布没什么影响。这个结论与猜测完全相悖，现非常沮丧，需要进一步的实验和研究。

今天有了一个新的想法，继续进行了研究。我认为广凭借图像观看是不合适的，必须拿出数字来。我考虑了一下，决定把直方图当作向量，计算一下原分辨率影像归一化了像元数量之后的直方图与缩小两倍、缩小四倍之间的欧式距离。

好，增加代码。

# 判断一下hist_1 与 hist_2 之间的欧式距离与hist_1与hist_4之间的几何距离之差
def calc_l2_distance(hist_1, hist_2):
    distance_sum = 0
    for i in range(256):
        distance = np.sqrt(np.abs(hist_1[i] - hist_2[i])**2)
        distance_sum += distance
    return distance_sum

l2_distance_1 = calc_l2_distance(hist_1, hist_2)
l2_distance_2 = calc_l2_distance(hist_1, hist_4)
print("The Euclid between hist_1 and hist_2 is: ", l2_distance_1)
print("The Euciid between hist_1 and hist_4 is: ", l2_distance_2)

结果：

The Euclid between hist_1 and hist_2 is:  399.0
The Euciid between hist_1 and hist_4 is:  644.0

好，这下事实证明，在分布上，相似分辨率的影像更接近。也就是说，不同尺度的遥感影像分布还是具有差异性的。这种差异性，是高分辨率遥感影像更丰富的纹理细节等促成的。

那么我们根据这个结论预言，hist_1和hist_8之间的欧式距离应该更大才对。下面我们来实验一下是不是。

l2_distance_3 = calc_l2_distance(hist_1, hist_8)
print("The Euciid between hist_1 and hist_8 is: ", l2_distance_3)

结果：

The Euciid between hist_1 and hist_8 is:  1308.0

这证明了我们的结论是正确的。