使用 ENVI + IDL 开发自定义处理算法 - 如此简单！
19695 对此文章评分：
4.2
使用 ENVI + IDL 开发自定义处理算法 - 如此简单！
作者：匿名 2015年3月19日，星期四
下周我将在我们位于科罗拉多州博尔德的办公室教授一门名为“使用 IDL 扩展 ENVI”的课程。许多人可能不知道，我们的地理空间图像分析软件 ENVI 是建立在一个名为 IDL 的强大编程语言之上的。IDL 被多个学科用于从复杂的数值数据中提取有意义的信息并创建可视化。IDL 的前身版本是为分析来自 NASA 任务（如水手号和国际紫外探测器）的数据而编写的。这些 IDL 的前身版本于 1970 年代在科罗拉多大学博尔德分校的大气与空间物理实验室 (LASP) 开发。IDL 的第一个官方版本于 1977 年由 Research Systems Inc. (RSI) 发布。自 1970 年代创建以来，IDL 已经走过了漫长的道路。我们的工程师不断为 IDL 添加新功能，以使其保持现代化、灵活性和用户友好性。
由于 IDL 作为一种编程语言，其根源面向处理科学数据多维数组的科学家和工程师，因此它非常自然地适用于处理来自遥感系统的栅格数据。1994 年，RSI 发布了 ENVI 的第一个版本。从那时起，ENVI 和 IDL 就像形影不离的好朋友。IDL 能做什么，ENVI 也能做什么。ENVI 有一个用 IDL 编写的 API，可用于创建自定义算法、自定义工作流以及对遥感数据进行批处理。
在我看来，ENVI 能够使用 IDL 进行扩展是其最强大的一个方面。ENVI API 学习曲线平缓，一旦掌握，你处理数据的能力将不受限制。这包括地理数据的强大可视化，例如我使用 CONTOUR 函数创建的下面这个。
ENVI API 更强大的一个方面是能够实现自定义算法。这允许用户应用他们在科学文献中读到的算法，或者试验他们自己的图像分析算法。让我们直接来看看如何使用 IDL 在 ENVI 中实现自定义算法。在本例中，我们将对 Landsat 8 使用缨帽变换。TCT 最初是为 1972 年发射的陆地卫星多光谱扫描仪 (MSS) 设计的。随后，针对陆地卫星专题制图仪 (TM)、陆地卫星增强型专题制图仪 (ETM) 和陆地卫星增强型专题制图仪+ (ETM+) 传感器的 TCT 改编版本已在科学文献中发表。ENVI 开箱即用地支持对 Landsat MSS、TM、ETM 和 ETM+ 图像执行 TCT 的能力。然而，用于陆地卫星操作陆地成像仪 (OLI) 传感器的 TCT 算法直到最近才发表，尚未进入 ENVI 工具箱。但是，由于该算法现在已在科学文献中可用，我们可以使用 ENVI API 将其应用于 Landsat 8 数据。该算法获取 Landsat 8 图像并计算一个包含 6 个波段的新图像文件，每个波段包含关于场景的独特信息，例如反照率（亮度）、植被健康状况（绿度）和水分含量（湿度）。
下面这个示例代码可以在 ENVI + IDL 中实现对 Landsat 8 数据执行 TCT，该算法描述于题为 Derivations of a tasseled cap transformation based on Landsat 8 at-satellite reflectance 的文章中：
; Copyright (c) 1988-2015, Exelis Visual Information Solutions, Inc. 版权所有。
pro tasseledcapidl
; 设置编译选项 - 这是一个标准兼容性行，建议在所有 IDL 代码中使用
compile_opt IDL2
; 使用 ENVI 的当前会话
e=envi(/current)
; 预校准的 Landsat 8 图像文件的路径和文件名
; 输入文件应经过辐射定标，可选择进行大气校正
inputRaster = e.OpenRaster('C:\Data\Naples\LC80160422015028LGN00_RadCal_MS_Subset.dat')
; 对输入栅格进行子集提取 - TCT 计算只需要这 6 个波段
subsetRaster=ENVIsubsetRaster(inputRaster, BANDS=[1,2,3,4,5,6])
; 设置输出栅格文件
outputRaster = ENVIRaster(URI=outputURI, INHERITS_FROM=subsetRaster)
; 创建瓦片 - 创建瓦片允许 ENVI + IDL 迭代处理瓦片，这是一个好主意
; 当你想最小化图像处理对计算机资源的影响时
tiles = subsetRaster.CreateTileIterator(MODE='spectral')
; 遍历瓦片
FOREACH tile, tiles DO BEGIN
; 波段 1
b1 = tile[*,0]
; 波段 2
b2 = tile[*,1]
; 波段 3
b3 = tile[*,2]
; 波段 4
b4 = tile[*,3]
; 波段 5
b5 = tile[*,4]
; 波段 6
b6 = tile[*,5]
dims= size(tile,/DIMENSIONS)
data = make_array(dims,/FLOAT)
; 使用 IDL 进行 TCT 计算
; 计算亮度
data[*,0] = (b1 * 0.3029) + (b2 * 0.2786)+ (b3*0.4733) + (b4*0.5599) + (b5*0.508) + (b6*0.1872)
; 计算绿度
data[*,1] = (b1 * (-0.2941)) + (b2 * (-0.243))+ (b3 * (-0.5424)) + (b4 * 0.7276) + (b5 * 0.0713)
; 计算湿度
data[*,2] = (b1 * 0.1511) + (b2 * 0.1973)+ (b3 * 0.3283) + (b4 * 0.3407) + (b5 * (-0.7117)) + (b6*(-0.4559))
; 计算 TCT4（霾）
data[*,3] = (b1 * (-0.8239)) + (b2 * 0.0849)+ (b3 * 0.4396) + (b4 * (-0.058)) + (b5 * 0.2013) + (b6*(-0.2773))
; 计算 TCT5
data[*,4] = (b1 * (-0.3294)) + (b2 * 0.0557)+ (b3 * 0.1056) + (b4 * 0.1855) + (b5 * (-0.4349)) + (b6 *0.8085)
; 计算 TCT6
data[*,5] = (b1 * 0.1079) + (b2 * (-0.9023))+ (b3 * 0.4119) + (b4 * 0.0575) + (b5 * (-0.0259)) + (b6 *0.0252)
; 写入输出文件
outputRaster.SetTile, data, tiles
ENDFOREACH
; 向 HDR 文件添加适当的波段名称
metadata = outputRaster.METADATA
metadata.UpdateItem, 'BAND NAMES', ['Brightness','Greenness','Wetness','TCT4','TCT5','TCT6']
; 保存更改到输出文件
outputRaster.save
; 在 ENVI 显示中显示输出文件
view=e.getview()
layer=view.createlayer(outputRaster)
end
此示例的输入文件是 Landsat 8 于 2015 年 1 月 28 日在佛罗里达州那不勒斯上空拍摄的。下面来自 ENVI 显示的屏幕截图展示了场景的彩色红外表示，以及从缨帽变换计算出的亮度、绿度和湿度波段。这只是 IDL 的强大功能如何使实现自定义算法变得如此简单的一个例子。
参考文献:
Muhammad Hasan Ali Baig, Lifu Zhang, Tong Shuai & Qingxi Tong (2014) Derivations of a tasseled cap transformation based on Landsat 8 at-satellite reflectance, Remote Sensing Letters, 5:5, 423-431, DOI: 10.1080/2150704X.2014.915434