基于IDL的遥感影像波段运算

基于IDL的遥感影像波段运算1.IDL语言特征在IDL中，矩阵是按照列的方式处理的，即矩阵是以n列、m行的方式表示的，这一点与Fortran语言一样，但与C语言（按m行、n列的维度来标记）表示矩阵的方式不同。只有熟悉IDL的数据存储、处理的特点，才能方便地进行分析。以300列、200行的遥感数据band1为例，其数据的记录特征如下：IDL图形坐标band1DN值坐标图1.IDL中像元坐标与亮度矩阵的对应关系（1）数据按列存储；最先存储的是第一列的数据（图象第一列的像元亮度值），接着是第二列，直至最后一列。（show:print,data[0,*]）（2）采用顺序显示数据的方式。将第一行的亮度值按顺序从左到右以行的方式显示或存储，不管band1矩阵有多少行、多少列，总是先显示、存储完第一行，然后是下一行。由于IDL的默认显示宽度是80列，所以在输出窗口中要用几行才能显示完矩阵的第一行亮度值。然后另起一行显示亮度矩阵的下一行值，直到最后一行。这一点对于图象显示和分析相当重要。复习：IDL矩阵的显示方式，data（show:print,data[*,0]）（3）默认情况下，图象数组的第一行和第一个元素显示在屏幕的左下角。后面的行将从下往上显示。在图形窗口显示图象时，图象的每一个元素在屏幕中显示为一个像素[1]。在IDL中图象是从下往上画的。即把band1亮度矩阵的第一行值显示在图象最下端的一行像元中，把band1亮度矩阵的最后一行值显示在图象最上端的一行像元中。这是图象显示最重要的特征。区别两幅图象的显示效果：（tvscl,data）与（tvscl,data,300，0，order=1）（4）矩阵转置与图象的关系像元坐标亮度矩阵1,12,1……300,19294……901,22,2……300,29692……91…………………………1,2002,200……300,20098101……90将矩阵band1转置，(print,(transpose(data))[*,0])，打印转置矩阵的第一行，它在band1中是哪行或哪列？(print,transpose(data)tvscl,data,200,0,order=1结论：按行显示的数据DN，按列显示；按列显示的DN数据，按行显示2.遥感数据的处理与分析31通用的栅格数据存储格式BSQ：波段顺序存储格式。每行数据后面紧接着同一波谱波段的下一行数据。这种格式最适于对单个波谱波段中任何部分的空间（X，Y）存取。每个tile是单个波段的一个空间子集。BIP：波段按像元交叉格式。图像按顺序存储第1个像元所有的波段，接着是第2个像元的所有波段，然后是第3个像元的所有波段，等等，交叉存取直到像元总数为止。这种格式为图像数据波谱（Z）的存取提供最佳性能。每个tile是一幅图象中所有波段行的所有像素。BIL：波段按行交叉格式。按BIL格式存储的图像先存储第一个波段的第一行，接着是第二个波段的第一行，然后是第三个波段的第一行，交叉存取直到波段总数为止。每个波段随后的行按照类似的方式交叉存取。这种格式提供了空间和波谱处理之间一种折衷方式，它是大多数ENVI处理任务中所推荐的文件格式。每个tile是一幅图像所有波段的一行。32栅格数据输入以BSQ格式的图象为例。假设有一幅BSQ格式的7波段TM影象“2003_sw.img”，需要计算各波段亮度值的基本统计信息（包括：最大值、最小值、均值、均方差、倾斜度和平坦度）。遥感数据其实就是像元的亮度值（也叫做DN值），这与它有多少个波段、是多光谱数据还是高光谱数据都没有关系。这里分析遥感数据的目的，是为了帮助理解、纠正易于混淆的栅格数据存储方法和数据分析方法。影象“subset_2003.img”的大小为300×200（即：lines:200，Samples：300）。以band1波段为例，在IDL的图形显示系统中，左上角坐标为（1，1）。利用程序把band1的亮度值读入IDL环境中，构造成一个亮度矩阵（见图1）。在IDL的图形系统中，图象左上角的坐标为（1，1），右下角（X,Y）坐标值最大。在图1中，图象band1的像元坐标左上角为（1，1），右下角为（300，200），坐标对中的X为列号，Y为行号，像元数量为300列200行。与像元一一对应，亮度矩阵是band1的亮度值，每个像元一个值，左上角为92，右下角为90。亮度矩阵中的第一行值，与图象中的第一行像元从左到有一一对应，亮度矩阵中的第二行值，与图象中的第二行像元从左到有一一对应，依次类推，直到最后一行。在IDL中打开band1的亮度矩阵并且显示在IDL绘图环境中的时候，要正确地显示图象（见图2），必须要确保亮度矩阵与图象坐标系统的正确关系，如亮度矩阵中的左下三角（或右上三角）一定要严格对应图象像元矩阵的左下三角（或右上三角），否则就会发生图象翻转、倒置（见图3），其实质就是将波段亮度矩阵进行了转置或行列倒置。练习：程序调用。调用band4的ASCII文件的代码如下：file=filepath('band1_matrix.txt',subdir='examples\data')openr,lun,file,/get_lundata=lonarr(300,200);事先知道图象大小readf,lun,dataprint,size(data,/n_elements)；data就是band4图象的亮度矩阵free_lun,lun3.波段运算（BandMath）BandMathTM功能允许你处理导致单个波段输出的复杂表达式。这些数学表达式也可以应用于一个多波段文件中的所有波段，providing“FileMath”。关于使用波段运算的更多信息，请参阅ENVIProgrammer’sGuide第29页的“BandMathBasics”。1．可利用的波段运算功能（AvailableBandMathFunctions）BandMath功能为用户提供一个灵活的图像处理工具，其中许多功能是无法在任何其它的图像处理系统中获得的。该功能的能力与IDL语言的能力直接相关。可用的函数包括但不仅限于表4-2中列出的数学表达式。表4-2:一些可用的波段运算函数。SeriesandScalar数学三角函数其它波段运算选项加（+）正弦（sin（x））关系运算符（EQ、NE、LE、LT、GE、GT）减（-）余弦（cos（x））逻辑运算符（AND、OR、XOR、NOT）乘（*）正切（tan（x））类型转换函数（byte，fix，long，float，double，complex）除（/）反正弦（asin（x））IDL返回数组结果的函数最小运算符（）反余弦（acos（x））IDL返回数组结果的程序最大运算符（）反正切（atan（x））UserIDL函数和程序绝对值（abs（x））双曲正弦（sinh（x））平方根（sqrt（x））双曲余弦（cosh（x））指数（^）双曲正切（tanh（x））自然指数（exp（x））自然对数（alog（x））以10为底的对数（alog10（x））注意一些有效的IDL表达式要求整个输入数组存在于内存中，它可以不必与ENVItiling操作相兼容。2．BandMath对话框（1）.选择BasicToolsBandMath.将出现BandMath对话框。假如运算结果是一个二维数组，它将接受任何有效的IDL数学表达式、函数或程序。（2）.在标签为“Enteranexpression:”的文本框内，输入变量名（将被赋值到整个图像波段或可能应用到一个多波段文件中的每个波段）和所需要的数学运算符。变量名必须以字符“b”或“B”开头，后面跟着5个以内的数字字符。实例：若你想计算三个波段的平均值，则在文本框“Enteranexpression:”内输入数学方程式：（（float（b4）-float（b3））/（（float（b4）+float（b3））这时，变量b4、b3自动跳入”Previousbandmathexpression”对话框中，可以输入到文本框中。该表达式中使用的三个变量，“b3”是第一个变量，“b4”是第二个变量，注意，在本例中，IDL的浮点型函数用来防止计算时出现字节溢出错误。（3）.输入一个有效的表达式，点击“OK”处理。将出现Variable/BandNamePairings对话框。请参见以下部分。4．要重新使用、保存或取消任何以前应用的数学表达式：（1）.点击显示在“PreviousExpression:”列表中的任何表达式，把它导入到“Enteranexpression:”文本区中。（2）.一旦被导入，点击“OK”，把该表达式应用到一组新的波段。将出现Variable/BandNamePairings对话框。请见下列的详细向导。5．要把表达式保存到一个输出文件，点击“Save”，然后当出现EnterOutputFilename对话框时，键入输出文件名。为了保持一致，输出文件名应该指定扩展名为.exp。（Tips:可以在记事本中先输入复杂的公式，再导入公式文件）。6．要恢复原先保存的表达式，点击“Restore”，然后选择适当的文件名。该表达式将显示在“PreviousExpression:”列表中。·要清除所有原先的表达式，点击“Clear”。图4-13:BandMath对话框。Variable/BandNamePairings对话框Variable/BandNamePairings对话框允许你从一个输入波段列表中，把波段赋值给输入在“Enteranexpression:”文本框中的变量。·要把一个值赋给原先实例中的变量“b1”：1.在标签为“Variablesusedinexpression:”的文本框内，点击表达式“B1”。2.在标签为“可利用波段列表:”的列表中，点击所需要的波段。注意，一旦第一个波段被选择，只有那些相同空间大小的波段被显示在波段列表中。3.按照同种方法，为“B2”、“B3”等赋予一个值。·要把一个多波段图像赋值给一个或所有变量：1.点击“MapvariabletoInputfile”。2.使用标准的ENVI文件选择步骤，选择一个文件（这可视为“FileMath”）。所选择的文件可以是波谱子集，但是若一个以上的文件被使用，它们必须有相同的波段数。通过数学表达式修改的文件数学（filemath），一个多波段输出图像产生。3.一旦所有变量被定义，标准的ENVI输出对话框显示在Variable/BandsPairings对话框的底部。可以把一个波段运算结果作为新的变量输入下一个运算过程。·要选择一个空间子集：1.点击“SpatialSubset”。2.将出现标准的FileSpatialSubset对话框（系统默认值被设置为处理整个空间场景。·要把结果输出到一个文件或内存，选择“File”或“Memory”切换按钮。·若选择输出到一个文件，键入一个输出文件名，或使用“Choose”按钮选择一个文件名，然后点击“OK”。结果图像被显示在可利用波段列表中。将用公式计算的NDVI和用菜单命令直接计算的NDVI进行比较。4.运用IDL程序进行波段运算由于ENVI为你提供对IDL性能的访问，你可以使用内置的IDL功能部件的能力、IDL用户函数，或书写你自己的程序执行自定义的操作。这些函数的唯一要求是它们接受一个或多个图像阵列作为输入，并且输出一个单波段二维数组的计算结果。这些函数必须保存在IDL路径列表内的一个目录下，以便它们将自动编译。通过使用ENVI主菜单System下拉菜单下的CompileModule选项（见ENVIProgrammer’sGuide第23页的“IncorporatingNewRoutines”），也可以对它们进行编译。以下是用户波段运算功能的一些简单的实例。波段运算函数1下面的实例是一个非常简单的自定义波段运算函数，它把两个波段相加。下面的程序文本可以在一个文本编辑器中输入，并用文件名user_bm1.pro来保存：实例：FUNCTIONuser_bm1，b1，b2