转载 | 使用 Python 处理高光谱数据（上）

近年来，深度学习和大数据处理技术为高光谱分析带来了巨大的变革，使用深度学习处理高光谱数据已经成为研究的主流热点。ENVI提供了成熟高效的高光谱处理平台，但是相对稳定的软件结构降低了应用的灵活性。而且ENVI对于流行的深度学习架构支持不够完善，难以满足现在复杂的研究需要。本文将介绍使用Python处理高光谱数据的思路，并提供相关代码方便读者复现。

Python

Python是一门计算机语言。编程人员通过编写Python代码命令计算机执行指令，完成特定的任务。编写代码时需要遵循Python独特的风格，Python内部才能正确将这些代码转换为计算机可以理解的指令。Python有着非常活跃的社区，一般问题都能通过搜索引擎找到答案，这大大提高了使用Python进行编程的效率。 

由于其简单开放的特性，Python受到科学界的广泛支持，成为科研工作中Z 流行的语言。Python本身不具备科学计算的能力，而且运行效率较低。为此，科学家和工程师开发了完整的基于Python的科学计算生态系统。该生态系统包括矩阵库Numpy、表格库Pandas、科学计算库Scipy、绘图库Matplotlib以及代码编辑运行工具Jupyter等工具包。这些工具包不仅降低了使用Python进行科学计算的门槛，还通过底层优化，大大提高了运算速度。所以在使用Python进行高光谱图像进行处理时要同时学习这些工具，学习的优先级和掌握程度参考右图。

光谱库文件的读取

Python独立安装包内不包含科学计算的工具包，且需要繁琐的步骤进行开发环境的配置，所以推荐使用Anaconda进行Python程序的安装。Anaconda中包含了大部分科学计算工具，可以快速搭建高光谱数据处理平台。

下载地址：https://www.anaconda.com/products/individual

安装完成后重启。

打开 Anaconda Prompt Shell，输入以下命令，将目录切换到工作目录下，运行 Jupyter Lab，并依次点击 File -> New -> Notebook 新建文件：

将ASCII格式的光谱库文件拷贝到工作目录下。

首先，利用下面代码，通过逐行读取的方式将数据读取到内存中。

import csv

f='ascii.txt'

with open(f) as csvfile:

     data = list(csv.reader(csvfile)) 

     csvfile.close()

按 Shift + Enter 运行代码块并创建新的代码块，在新的代码块中输入变量名可直接查看该变量的值，这里输入并运行。观察数据的组织结构可以发现，该数据由两个部分组成分别是文件头和数据主体两个部分，这两个部分需要分开处理。

data

文件头通常由几行字符串组成。Python 本身对于字符串有很强大的处理能力，这里展示了两种文件头的处理方式：

ls1 = [] 

ls2 = [] 

# 使用字符串拆分手段获取列名 

for i in range(2,7): 

    s = data[i][0] 

    class_name = s.split(': ')[-1].split('~~')[0]      

    ls1.append(class_name) 

print(ls1) 

# 使用正则表达式匹配方法获取列名 

import re 

for i in range(2,7): 

    s = data[i][0] 

    class_name = re.search(r'.*: (.*?)~.*', s).group(1) 

    ls2.append(class_name) 

print(ls2)

# 需要注意数据中的空格

输出结果： 

['tree', 'soil', 'water', 'maize', 'road'] 

['tree', 'soil', 'water', 'maize', 'road']

本文件的数据主体部分由固定宽度的几列数字组成，对于这种组织形式可以使用 Pandas.read_fwf() 进行解析。

import pandas as pd

colspecs = [(0,13),(13,23),(23,33),(33,43),(43,53),(53,63)] 

table = pd.read_fwf('ascii.txt',skiprows=7,header=None,colspecs=colspecs) 

#在数据存储中推荐使用以逗号分隔符进行数据存储，这样数据更易于解析。例如Python本身就对逗号分割的数据具有解析功能，而固定宽度的数据需要调用复杂的函数，提高了数据传递的难度。

# 调整数据的显示

table = table.set_index(0)

table.columns = ls1

table.index.name = 'wavelength'

table

输出结果：

至此便完成了ENVI波谱库的解析工作。使用Python的科学计算库可以完成所有常规的高光谱处理工作， 如绘图，插值，统计，回归等。

# 绘图

table.plot()

输出结果：

# 四则运算与统计 

t_s = table['tree'] - table['soil'] 

table['tree'] = table['soil'] + t_s 

table['tree'] = table['tree']*100

table['tree'] = table['tree']/100

table.describe()

输出结果：

# 上下包络线计算 https://stackoverflow.com/questions/34235530/how-to-get-high-and- low-envelope-of-a-signal/34245942#34245942 

from matplotlib import pyplot as plt 

def hl_envelopes_idx(s, dmin=1, dmax=1, split=False): 

    """

    Input : 

    s: 1d-array, data signal from which to extract high and low envelopes dmin, dmax: int, optional, size of chunks, use this if the size of the input signal is too big 

    split: bool, optional, if True, split the signal in half along its mean, might help to generate the envelope in some cases

    Output : 

    lmin,lmax : high/low envelope idx of input signal s 

    """ 

    # locals min 

    lmin = (np.diff(np.sign(np.diff(s))) > 0).nonzero()[0] + 1 

    # locals max 

    lmax = (np.diff(np.sign(np.diff(s))) < 0).nonzero()[0] + 1 

    if split: 

      # s_mid is zero if s centered around x-axis or more generally mean of signal 

      s_mid = np.mean(s) 

      # pre-sorting of locals min based on relative position with respect to s_mid

      lmin = lmin[s[lmin]<s_mid] 

      # pre-sorting of local max based on relative position with respect to s_mid 

      lmax = lmax[s[lmax]>s_mid] 

    # global min of dmin-chunks of locals min  

    lmin = lmin[[i+np.argmin(s[lmin[i:i+dmin]]) for i in range(0,len(lmin),dmin)]] 

    # global max of dmax-chunks of locals max 

    lmax = lmax[[i+np.argmax(s[lmax[i:i+dmax]]) for i in range(0,len(lmax),dmax)]] 

return lmin,lmax 

y='tree' 

l_en,h_en = hl_envelopes_idx(table[y].values) plt.plot(table.iloc[h_en].index,table.iloc[h_en][y],label = 'High Envelope') plt.plot(table.index,table[y],label = 'data',color = 'black') plt.plot(table.iloc[l_en].index,table.iloc[l_en][y],label = 'Low Envelope') plt.legend()

输出结果：

小结

本文介绍了使用Python处理高光谱数据的优势，开发环境的安装，以及光谱库文件的读取与处理。千里之行，始于足下。高光谱图像是由网格化的单根光谱曲线组成的，其处理思路与单个光谱曲线的处理有很多一致性。所以有必要理解如何对单根光谱曲线进行处理，并将这些处理方法应用于高光谱影像。此外，高光谱影像作为地理影像，包含单根光谱曲线没有的空间信息，其数据量也达到了大数据的级别， 因此高光谱图像的处理又有其特殊的地方。