6.3.2 光谱曲线可视化

import warnings
warnings.simplefilter("ignore")

导入功能模块

import sys

load_dataset_module_path = '../../'
sys.path.append(load_dataset_module_path)

from load_hyperspectral_dataset import (load_hyperspectral_data, y_labels,
                                        extract_features,
                                        get_label,
                                        plot_selected_categories,
                                        plot_decision_function)

#%matplotlib inline
#%matplotlib notebook

#%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
from sklearn.utils import shuffle
import seaborn as sns
from ipywidgets import interact_manual
import numpy as np

mpl.style.use('ggplot')

加载数据集

img_file_path = '../../Hyperspectral_Project/dc.tif'
label_file_path = '../../Hyperspectral_Project/dctest.project'

raw_X, raw_y, pixel_position = load_hyperspectral_data(img_file_path,
                                                       label_file_path)

hyperspectral_df =  extract_features(raw_X, raw_y)

y_labels

['Roofs', 'Street', 'Path', 'Grass', 'Trees', 'Water', 'Shadow']

sum(raw_y==0)

3834

画出指定类别，训练模型分类

label_list, region_list = get_label(label_file_path)

Roofs_region_list = [region for label,region in zip(label_list,region_list)
                          if label==0]

len(Roofs_region_list)

38

#ax = plt.gca()
#colors = sns.color_palette("Paired",n_colors=38)
#colors = sns.color_palette("husl", n_colors=38)

marker_list = ['o','s','d','^','*','h']
fig = plt.figure()
for i,r in enumerate(Roofs_region_list):

    is_in_r = np.where(
        (pixel_position[:, 0] >= r[0]) &
        (pixel_position[:, 0] <= r[2]) &
        (pixel_position[:, 1] >= r[1]) &
        (pixel_position[:, 1] <= r[3]))[0]

    #print(i)
    plt.scatter( hyperspectral_df[['MNDWI']].iloc[is_in_r],
                hyperspectral_df[['NDBI']].iloc[is_in_r],
                #cmap=sns.color_palette("husl", n_colors=38,as_cmap=True),
                #c = [colors[i]],
                marker = marker_list[i//7],
                s=20,
                label=i+1,
               alpha=0.6,
                edgecolors='white',
                linewidths=0.1
                )
    ax = plt.gca()
    ax.set(xlim=(0.4, 0.95), ylim=(-0.8, -0.56))

    #plt.show()
#xlim = ax.get_xlim()
#ylim = ax.get_ylim()

#print(xlim,ylim)

plt.legend(loc='lower left',
           #fontsize=3,
           ncol=3,
           frameon=False,
           bbox_to_anchor=(1, -0.01))

plt.savefig('Roofs_MNDWI_NDBI.pdf', bbox_inches='tight')

marker_list = ['o','s','d','^','*','h']

for animation_index in range(38):

    fig = plt.figure()

    for i,r in enumerate(Roofs_region_list):

        is_in_r = np.where(
            (pixel_position[:, 0] >= r[0]) &
            (pixel_position[:, 0] <= r[2]) &
            (pixel_position[:, 1] >= r[1]) &
            (pixel_position[:, 1] <= r[3]))[0]

        if i<= animation_index:
            x = hyperspectral_df[['MNDWI']].iloc[is_in_r]
            y = hyperspectral_df[['NDBI']].iloc[is_in_r]
        else:
            x = []
            y = []

        plt.scatter( x,
                         y,
                        marker = marker_list[i//7],
                        s=20,
                        label=i+1,
                       alpha=0.8,
                        edgecolors='white',
                        linewidths=0.1
                        )
        ax = plt.gca()
        ax.set(xlim=(0.4, 0.95), ylim=(-0.8, -0.56))

    plt.legend(loc='lower left',
               #fontsize=3,
               ncol=3,
               frameon=False,
               bbox_to_anchor=(1, -0.01))

    plt.savefig(f'Roofs_MNDWI_NDBI/Roofs_MNDWI_NDBI-{animation_index+1}.jpg',
                bbox_inches='tight',
                quality=95,
               dpi=300)

import imageio

file_names = [
    f'Roofs_MNDWI_NDBI/Roofs_MNDWI_NDBI-{i+1}.jpg'
    for i in range(38)
]

frames = []

for image_name in file_names:
    frames.append(imageio.imread(image_name))

imageio.mimsave(
    f'Roofs_MNDWI_NDBI_animation.gif',
    frames,
    'GIF',
    duration=0.33)

def plot_Roofs_MNDWI_NDBI(start,stop):

    for i,r in enumerate(Roofs_region_list[start-1:stop]):

        is_in_r = np.where(
            (pixel_position[:, 0] >= r[0]) &
            (pixel_position[:, 0] <= r[2]) &
            (pixel_position[:, 1] >= r[1]) &
            (pixel_position[:, 1] <= r[3]))[0]

        #print(i)
        plt.scatter( hyperspectral_df[['MNDWI']].iloc[is_in_r],
                    hyperspectral_df[['NDBI']].iloc[is_in_r],
                    #s=10,
                    label=i+start,
                    edgecolors='white',
                    )
        ax = plt.gca()
        ax.set(xlim=(0.4, 0.95), ylim=(-0.8, -0.56))


    plt.legend(loc='lower left',
               #fontsize=3,
               ncol=3,
               frameon=False,
               bbox_to_anchor=(1, -0.01))

    plt.savefig(f'Roofs_MNDWI_NDBI_{start}-{stop}.pdf', bbox_inches='tight')

plot_Roofs_MNDWI_NDBI(start=1,stop=7)

plot_Roofs_MNDWI_NDBI(start=8,stop=14)

plot_Roofs_MNDWI_NDBI(start=9,stop=21)

plot_Roofs_MNDWI_NDBI(start=22,stop=28)

plot_Roofs_MNDWI_NDBI(start=29,stop=35)

plot_Roofs_MNDWI_NDBI(start=35,stop=38)

from skimage.io import imread, imsave
import numpy as np
from skimage.exposure import rescale_intensity

img = imread(img_file_path)


#选取第60、27和17波段，组合成RGB图像
img_RGB = np.dstack([img[59, :, :], img[26, :, :], img[16, :, :]])
#对图像进行对比度拉伸，放缩到0-255的数值范围
p_10, p_90 = np.percentile(img_RGB, (10, 90))
img_rescale = rescale_intensity(img_RGB,
                                in_range=(p_10, p_90),
                                out_range=(0, 255))
img_rescale = img_rescale.astype('uint8')

from skimage import color
from skimage.io import imshow
from matplotlib import patches

gray_img = color.rgb2gray(img_rescale)
gray_img = gray_img.reshape(img.shape[1], img.shape[2], 1)
gray_img = gray_img.repeat(3, axis=2)


label_list, region_list = get_label(label_file_path)

#colors = sns.color_palette("husl", n_colors=38)
colors = sns.color_palette('bright')

for i,region in enumerate(Roofs_region_list):
    gray_img[region[0]:region[2], region[1]:region[3]] = colors[3]
    plt.text(region[1],region[0]+4, f'{i+1}',fontsize=1.2,color='white')

imshow(gray_img)
#labels = ['Roofs', 'Street', 'Path', 'Grass', 'Trees', 'Water', 'Shadow']
legend_patches = []
legend_patches.append(patches.Patch(color=colors[3],label='Roofs'))

plt.legend(handles=legend_patches,
           loc='lower left',
           fontsize=3,
           frameon=False,
           bbox_to_anchor=(1, -0.01))


plt.axis('off')
plt.savefig('Roofs_number_image.pdf', bbox_inches='tight', dpi=600)

from load_hyperspectral_dataset import get_wavelength_table

wavelength_table = get_wavelength_table()

y_region = np.array(region_list)

pixel_mean_value_list = []

linestyle_list = ['--','-.']
marker_list = [None,'s','D','^']

for i,region in enumerate(Roofs_region_list):
    region_img = img[:,region[0]:region[2],region[1]:region[3]]
    region_img = region_img.transpose((1,2,0)).reshape(-1,img.shape[0])


    region_mean_value = region_img.mean(axis=0)

    if (i//7)<4:
        marker = marker_list[i//7]
        linestyle= '-'
    else:
        marker=None

        linestyle = linestyle_list[(i//7)-4]



    plt.plot(wavelength_table,
             region_mean_value,
             linewidth=0.5,
             markersize=0.5,
             label=f'Roofs {str(i+1)}',
             linestyle=linestyle,
             marker=marker)



plt.legend(ncol=3,
           loc='lower left',
           bbox_to_anchor=(1, -0.02))
plt.xlabel('Wave length (nm)')
plt.savefig('Roofs_mean_spectral_curve.pdf', bbox_inches='tight')

from load_hyperspectral_dataset import get_wavelength_table

wavelength_table = get_wavelength_table()

y_region = np.array(region_list)


def plot_Roofs_mean_spectral_curve(start,stop):

    pixel_mean_value_list = []


    for i,region in enumerate(Roofs_region_list[start-1:stop]):
        region_img = img[:,region[0]:region[2],region[1]:region[3]]
        region_img = region_img.transpose((1,2,0)).reshape(-1,img.shape[0])


        region_mean_value = region_img.mean(axis=0)


        plt.plot(wavelength_table,
                 region_mean_value,
                 label=f'Roofs {str(start+i)}')



    plt.legend(
               loc='lower left',
               bbox_to_anchor=(1, -0.02))
    plt.xlabel('Wave length (nm)')
    plt.savefig(f'Roofs_mean_spectral_curve_{start}-{stop}.pdf', bbox_inches='tight')

plot_Roofs_mean_spectral_curve(start=1,stop=7)

plot_Roofs_mean_spectral_curve(start=8,stop=14)

plot_Roofs_mean_spectral_curve(start=15,stop=21)

plot_Roofs_mean_spectral_curve(start=22,stop=28)

plot_Roofs_mean_spectral_curve(start=29,stop=35)

plot_Roofs_mean_spectral_curve(start=36,stop=38)