IDL-程序设计笔记.docx

目录Part one：文件读写操作11、格式化输入与输出（readprint）12、ASCII文件读写23、二进制文件读写54、图像格式文件读写75、科学数据格式读写*8Part two：图形绘制11（一）plot过程绘制11（二）plotg()函数绘制12（三）散点图12（四）柱状图、条形图15（五）绘制颜色条colorbar18Part three: ENVI二次开发201、自定义波段运算函数202、IDL中常见波段运算213、ENVI classic 创建自定义菜单224、ENVI classic+IDL 二次开发23 Part one：文件读写操作1、格式化输入与输出（readprint）1) Read:该函数用于从键盘键入数据，默认数据类型为浮点型，若要输入其他类型数据，需要先定义数据类型。Eg:2) Print：该函数用于将数据输出打印到控制台，这里主要讲格式化控制输出format关键字。表01 常用格式化控制符格式代码输出效果aN字符或者字符串按照N个字符宽度输出fn.m浮点型数组按照N个字符宽度输出，小数点保留M位dn.mdouble型数组按照N个字符宽度输出，小数点保留M位en.m按照N个字符宽度的科学计数法输出，小数点后精确到M位Nx输出N个空格字符串/H直接引用字符串输出或者用H直接输出字符串c( )用于输出日期数据Eg：2、ASCII文件读写 常见ASCII文件：txt、matlab程序文件(*.m)、c程序文件（*.c）、遥感影像头文件（如ENVI格式的头文件*.hdr）1） 选择文件：fn = dialog_pickfile(title=选择ASCII文件：,get_path = work_dir)cd,work_dir ;用于跳转到当前工作路径2） 查询文件：file_search()函数Eg:查找当前工作路径下所有的txt文件，并返回文件数目*: File_lines(fname)函数用于查询文本文件的行数。3） 打开文件：在IDL中读写ASCII码或者二进制文件，首先需要将文件与一个逻辑设备号关联起来。语法：openr/openw/openu, lun, fname, /get_lun, width=10,/append4） 读写操作：IDL中利用readf函数读取文本文件、利用printf函数写入ASCII文件 Readf, lun, var01, var02. printf, lun, var01, var02.5） 文件关闭：在对文件操作完成之后，需要关闭文件的逻辑设备号，eg: freee_lun，lunEG:*如何获取某个文件你的列数:先读取一行数据，然后利用空格进行字符串拆分；最后利用n_elements()计算元素个数即为ns对于某些遥感数据的头文件前面几行是说明文件，后一部分才是数据的，读取有两种方式：1） 定义一个临时变量temp，将说明文字用该变量存储2） 直接跳行：skip_lun, lun, 3, /lines ;跳过文件的前三行3、二进制文件读写多数遥感数据的数据文件都是二进制文件。1) 函数readu, lun, var01, var02.用于读取二进制文件*：对于某些数据文件开头是说明信息，读取时可以用point_lun, lun, position直接按字节跳过。Eg:读取IDL文件下的envi.img数据2) 函数writeu, lun, var01, var02用于写入二进制文件。Eg01:将上面读取的envi.img计算其NDVI并保存为二进制文件Eg02:Eg03:;读取cha06下的风云02卫星数据pro read_AWX file = E:zengskIDLCha06dataFY2C_TBB_IR1_OTG_20061130_AOAD.AWX openr, lun, file, /get_lun ;跳过前20个字节 point_lun,lun,20 Headline = indgen(3);定义三个整型数据，分别读取记录长度、头文件记录数、数据记录数 readu, lun, Headline print,headline0, headline1,headline2 ;定位到第58个字节，读取数据日期 point_lun, lun, 58 Begindate = indgen(5);记录开始日期 readu, lun, Begindate Enddate = indgen(5);记录结束日期 readu, lun, Enddate lat_long = indgen(4);记录网格左上和右上角的经纬度 readu, lun, lat_long print,Begindate print,Enddate print,lat_long ;基于Headline信息数组，定义字节数组，其行列数分别为记录长度和数据记录数 data = bytarr(Headline2,Headline0) ;定位到数据部分，其头文件结束的字节位置为“头文件记录数*记录长度” point_lun, lun, Headline0 * Headline1 print,Headline0 * Headline1 ;读取数据部分 readu, lun, data help,data ;显示数据 window, 0, xsize = 500, ysize = 500 tv, congrid(data, 500, 500)end4、图像格式文件读写IDL自带了丰富的多种图像读写函数，如BMP，JPG，PNG, JPEG2000， DICOM，TIFF等相关函数：文件查询：query_image(fname, dimensions = dimensions, channels = nb ) ;channels返回波段数、dimensions返回图像的行列数 读取图像：read_image ( fname );结果返回图像数组 写入图像：write_image, fname(可带路径), jpg/bmp/tiff ，data, /order ; order关键字用于设置图像的纵坐标从上往下算起，默认为从下往上。*注： I、以上三个函数对于IDL支持的图像格式都适用。 II、显示图像时：tv用法：tv, data, true = 1/2/3 ;其中关键字true = 1表示数据格式为（3, m, n）、true = 2表示数组格式为（m, 3, n）、true = 3表示格式为（m, n, 3） III、tvrd:屏幕拷贝函数，返回当前直接图形窗口或设备的指定矩形部分的内容。IV、对遥感图像而言，其数据存放方式有BSQBIPBIL，以一个三维m列n行的遥感数据为例：BSQ的表示方式为m, n, 3、BIP表达方式为3, m ,n、BIL的表达方式为m, 3, n，常用的图像文件一般是BIP方式存储，遥感文件一般是BSQ。1） JPEG格式：读read_jpeg写write_jpegEg01:Eg02:2） BMP格式：(不采用任何压缩)，数据存储量一般很大读：read_imge(fname) or read_bmp(fname, /rgb)写：write_bmp, fname, data, /rgb3） TIFF格式： TIFF可以存储多波段图像，还可以包含投影信息，如landsat-7就是直接用.tiff存储的。Eg01: 5、科学数据格式读写*IDL中支持的科学数据格式包括CDF、HDF、HDF5、HDF-EOS、NetCDF等等。1） HDF4文件： 目前国内外的多种卫星传感器都是将HDF作为标准数据格式，包括EOS/MODIS、EOS/OMI、HJ-01/HSI、FY-3/MERSI。EG01:;读取E:ENVIIDL下的MODIS数据的经纬度以及其1KM的反射率数据集pro test_hdf fname=dialog_pickfile(title =选择数据文件：,get_path=cur_dir) cd, cur_dir hdf_id=hdf_sd_start(fname);打开hdf文件，返回一个文件id hdf_sd_fileinfo,hdf_id,sd_nums,attribute;用于获取HDF文件的数据集数目和属性数目 print,sd_nums,attribute ;读取经度 lat_index=hdf_sd_nametoindex(hdf_id,Latitude);根据数据集名称来获取数据集的索引号 lat_id=hdf_sd_select(hdf_id,lat_index);利用索引号选择数据集，返回一个数据集的id hdf_sd_getinfo,lat_id,name=name,unit=unit;查询已打开的数据集的基本信息 print,name, , unit hdf_sd_getdata,lat_id,lat;读取数据，参数lat用于返回读取的数据 hdf_sd_endaccess,lat_id;关闭已经打开的数据集 help,lat ;读取纬度 lon_index=hdf_sd_nametoindex(hdf_id,Longitude) lon_id=hdf_sd_select(hdf_id,lon_index) hdf_sd_getinfo,lon_id,name=name,unit=unit print,name, , unit hdf_sd_getdata,lon_id,lon hdf_sd_endaccess,lon_id help,lon ;读取1KM反射率数据 ref_index=hdf_sd_nametoindex(hdf_id,EV_1KM_RefSB) ref_id=hdf_sd_select(hdf_id,ref_index) hdf_sd_getinfo,ref_id,name=name,unit=unit print,name, , unit hdf_sd_getdata,ref_id,ref_value hdf_sd_endaccess,ref_id help,ref_value hdf_sd_end, hdf_id;关闭打开的文件 end运行效果2）*.h5、*.hdf5、*.he5格式的读取：;下面介绍一个读取环境卫星*.h5数据格式的例子pro read_hdf5 fname=dialog_pickfile(title =选择数据文件：,get_path=cur_dir) cd, cur_dir ;打开一个图形用户界面查看h5文件 result=h5_browser(fname) ;h5_message=h5_parse(fname);查询文件的基本信息，返回一个结构体变量 ;print,h5_message h5_id=h5f_open(fname);打开一个h5文件，返回一个文件的id sd_id=h5d_open(h5_id,ImageData/BandData);打开相应的数据集,数据集名称可带路径 Bandata=h5d_read(sd_id);读取数据 ;关闭文件 h5d_close,sd_id h5f_close,h5_id help,Bandata ;将Bandata转成BIP格式存储 Bandata_BIP = transpose(Bandata, 2,0,1) help,Bandata_BIP write_tiff,E:ENVIIDLHJ-1.tif, Bandata_BIPend2） netCDF文件读取：netCDF(network Common Data Form, 网络通用数据格式)，常用于气象科学数据存储。;下面介绍一个读取ECMWF数据的例子，ECMWF为*.nc格式数据;pro test_nc fname=dialog_pickfile(title =选择数据文件：,get_path=cur_dir) cd, cur_dir nc_id=ncdf_open(fname);打开netCDF文件 nc_message=ncdf_inquire(nc_id);对打开的文件进行查询，返回文件基本信息 help,nc_message ;变量查询 ;后一个参数为变量索引号或者名称;这里我读取的是它的第四个变量V10 var_message=ncdf_varinq(nc_id,4) help,var_message var_id=ncdf_varid(nc_id, v10);根据变量名称获取对应变量的索引号 print,var_id ;ncdf_varget获取数据 ncdf_varget,nc_id,var_id,value ;文件关闭 ncdf_close, nc_id help,valueend4）Grid格式读取：Part two：图形绘制（一）plot过程绘制1、plot过程：1）plot, xdata, ydata, /nodata 2）oplot, xdata, ydata ;用于在现有窗口中添加新曲线2、线形符号设置：1）关键字psym修改数据点的符号，直方图线形psym=10， psym= - a表示除显示各数据点的符号外，还将各数据点连接起来。2）linestyle关键字设置曲线线形，0为实线、2为虚线、3为点划线。3、 坐标轴设置：xytitle设置标题；xycharsize设置字体大小；xyrange用于设置坐标轴范围，注意加上xystyle=1来强行设定坐标轴范围xyticks和xyminor分别设置主刻度和最小刻度间隔数目xytickname设置刻度的名称4、 颜色设置：color: 设置图像颜色；一般用十六进制来表示颜色；FFFFFFxl表示白色；000000xl为黑色background：设置背景颜色5、 添加标注： xyouts过程用于在图像窗口中添加标注信息。示例：xyouts, 0.18, 0.80, RMSE=, color=FFFFFFxl, charsize=1.2, /normal6、 图形保存为文件：一般方法为利用TVRD函数拷贝图形窗口的内容，在写入图像文件。img = tvrd(x, y, channel=value, true=123, /order)（二）plotg()函数绘制1、plot函数：curve=plot(x, y, /buffer, /current, dimensions=width, height, margin=num, title= , name= , /overplot, window_title= , /nodata)*notice：/buffer指将图形保存在缓存中/current设置在当前窗口中绘制图形dimensions设置窗口的大小 title设置图形标题，name设置图像对象的名称*函数plot的返回结果为一个图像对象，具有多种对象和方法2、符号、线形设置：symbol、linestyle、thick用于改变曲线线宽3、坐标轴设置：见plot过程4、添加标注：函数text用于在图形窗口中添加标注信息，结果返回一个图形对象Eg: label=text(x, y, 标注内容, target=curve, color=red, font_size=1.2)5、添加图例：函数legend用于在图形窗口中添加标注信息，该函数只适用于图形对象（三）散点图Eg01：读取文本文件中的数据并利用plot过程绘制相应的散点图;利用过程plot绘制散点图pro plot_scatter ;*读取数据* fname=dialog_pickfile(title=选择数据：, get_path=work_dir) cd, work_dir nl=file_lines(fname) print, nl temp_var= data=fltarr(2, nl-1) openr,lun,fname,/get_lun result=fstat(lun) help,result readf, lun, temp_var readf, lun, data free_lun, lun print,temp_var print,data ;*绘制散点图* ; x = data0,*;观测数据 y = data1,*;估算数据 plot, x, y, psym=2,xrange=5,15,yrange=5,15,$ xtitle=Observed value, ytitle=Estimated value, $ background=FFFFFFxl, color=ooooooxl, $ xstyle=1, ystyle=1, $ charsize=1.3, /nodata oplot, x, y, psym=5, color=000000xl ;绘制基准线 x2=5, 15 & y2=5, 15 oplot, x2, y2, color=FF0000xl ;*计算MAE与RMSE并在散点图中添加标注* ; MAE=mean(abs(x-y);MAE平均绝对误差 RMSE=sqrt(mean(x-y)2);均方根误差 MAE_label=MAE=+string(MAE, format=(f5.2) RMSE_label=RMSE=+string(RMSE, format=(f5.2) ;添加标注 xyouts, 6, 13.0, MAE_label, color=000000xl, charsize=1.2, /data xyouts, 6, 12.3, RMSE_label, color=000000xl, charsize=1.2, /data ;保存散点图为图像文件，格式为.png img=tvrd(true=1) ;法一： o_fn=dialog_pickfile(title=文件保存为：);键入的文件名要包含后缀名 write_png, o_fn, img ;法二 ;write_png, E:ENVIscatter.png, imgendEg02: 读取文本文件中的数据并利用plot函数绘制相应的散点图;plot函数绘制散点图pro plot_scatter02 ;读取数据 file=dialog_pickfile(title=选择数据：,get_path=work_dir) cd, work_dir var= nl=file_lines(file) data=fltarr(2, nl-1) openr, lun, file, /get_lun readf, lun, var readf, lun, data free_lun, lun print, data ;绘图 ob=data0, * Es=data1, * graphic01=plot(ob, Es, xrang=7,15,yrang=7,15, xminor=5, yminor=5, $;先数据 title=scatter, xtitle=Observed Value, ytitle=Estimated Value, $;再坐标轴 symbol=X, sym_size=1.2, linestyle=none, color=black, $;然后线形 margin=0.1) x2=7,15 & y2=7,15 graphic02=plot(x2, y2, linestyle=1, /overplot) ;计算平均绝对误差和均方根误差并添加标注 MAE=mean(abs(Es-ob); RMSE=sqrt(mean(Es-ob)2) MAE_label=MAE=+string(MAE, format=(f5.2) RMSE_label=RMSE=+string(RMSE, format=(f5.2) t01=text(0.15, 0.8, MAE_label, target=graphic01, font_size=12);fonts_size用于设置标注的字体大小 t01=text(0.15, 0.75, RMSE_label, target=graphic01, font_size=12) ;保存文件 o_fn=dialog_pickfile(title=图形保存为：) graphic01.save, o_fn graphic02.save, o_fnend（四）柱状图、条形图Eg01: *利用bar_plot()函数绘制柱状图（不建议用直接图形法来绘制） *利用error_plot添加误差线;barplot功能函数绘制柱状图pro test_histogram02 ;#读取文本文件数据# file=dialog_pickfile(title=选择文本文件：, get_path=work_dir) cd, work_dir nl=file_lines(file) var=strarr(1,nl) openr, lun, file, /get_lun readf, lun, var free_lun, lun print, var ;拆分 area=strarr(nl) number=intarr(nl) for i=0, nl-1 do begin areai=strmid(vari, 0, 8) numberi=fix(strmid(vari, 11, 3) endfor help, area print, area help, number print, number3 ;#画柱状图# baselines=replicate(65,nl) His=barplot(number, $ bottom_values=baselines, xticklen=0, xrange=-0.5, 5.5, xtickname=area, $ yrange=65, 100, title=Histogram, ytitle=Impact Facter, $ width=0.5, dimensions=900, 500, $ bottom_color=EEEEEExl, fill_color=AOAOFFxl, $ margin=0.2) ;#绘制误差线# error=1.7, 2.7, 2.5, 2.4, 1.2, 1.0 error_lines=errorplot(number, error, linestyle=6, errorbar_capsize=0.2, $ errorbar_color=blue, /overplot) ;#文件保存为jpg# o_fn=dialog_pickfile(title=图片保存为：)+.jpg;注意后面加一个后缀名字符串更加方便 His.save, o_fn error_lines.save, o_fn end运行结果：Eg02: *在利用bar_plot()函数绘制图像时，添加 /horizontal关键字用于将图像横过来，即绘制条形图;bar_plot函数绘制条形图chart;pro test_chart ;#读取文本文件数据# file=dialog_pickfile(title=选择文本文件：, get_path=work_dir) cd, work_dir nl=file_lines(file) var=strarr(1,nl) openr, lun, file, /get_lun readf, lun, var free_lun, lun print, var ;拆分 area=strarr(nl) number=intarr(nl) for i=0, nl-1 do begin areai=strmid(vari, 0, 8) numberi=fix(strmid(vari, 11, 3) endfor help, area print, area help, number print, number3 ;#画柱状图# baselines=replicate(65,nl) His=barplot(number, bottom_values=baselines, $ yticklen=0, yrange=-0.5, 5.5, ytickname=area, $ xrange=65, 100, title=Histogram, xtitle=Impact Facter, $ width=0.5, dimensions=700, 500, $ bottom_color=EEEEEExl, fill_color=AOAOFFxl, $ margin=0.2, /horizontal) ;#文件保存为jpg# o_fn=dialog_pickfile(title=图片保存为：)+.jpg;注意后面加一个后缀名字符串更加方便 His.save, o_fnEnd运行结果：（五）绘制颜色条colorbarEg: 1）创建一个快速可视化窗口IDL w = window(dimensions=600, 400)2）; 创建一个Colorbar，有自定义颜色表，不设置标注位置。IDL cb = colorbar(POSITION = 0.05,0.85,0.7,0.9, RGB_TABLE= 72)3) ; 创建一个Colorbar具有自定义颜色表，设置标标注在顶部texpos=1。IDL cb = colorbar(POSITION = 0.05,0.85,0.7,0.9,RGB_TABLE=72, TEXTPOS=1)4) ; 创建Colorbar具有自定义颜色表, 设置自定义范围自定义范围range=*,*。 IDL cb = colorbar(POSITION = 0.05,0.75,0.7,0.8, RGB_TABLE = 72,RANGE = 0,1)5) ; Colorbar具有自定义范围和自定义标签tickname= A array。IDL cb = COLORBAR(POSITION=0.05,0.55,0.7,0.6, RGB_TABLE=72, RANGE=500,1000,TICKVALUES=550,650,750,850,950, TICKNAME=A,B,C,D,E6) ;创建一个colorbar,设置垂直色条orientation=1，范围0-100range，锥形端tapar=1（左图）。IDL cb = COLORBAR(POSITION=0.8,0.2,0.85,0.9, TAPER=1, RGB_TABLE=72, RANGE=0,100, ORIENTATION=1) *;离散色彩条垂直方向;标注在右边textpos=1，添加轮廓线BORDER=1（见右图）。IDL cb = COLORBAR(POSITION=0.9,0.15,0.95,0.95, RGB_TABLE=rgb,IDL TICKNAME=tickname, BORDER=1, TAPER=1, $ ORIENTATION=1, TEXTPOS=1) 7) ; 加载颜色表的一个子集，注意：这不是一个“离散”的色条因为每个颜色没有刻度标签。IDL LOADCT, 72, RGB_TABLE = rgbIDL rgb = rgb25 + 25*INDGEN(10), *IDL cb = COLORBAR(POSITION=0.05,0.45,0.7,0.5, RGB_TABLE=rgb)8）;创建一个离散颜色条，颜色条的颜色数量与标注的数量匹配，标注集中在颜色之下。并且添加轮廓线BORDER=1。IDL tickname = STRING(10*INDGEN(10), FORMAT=(I0)IDL cb = COLORBAR(POSITION=0.05,0.35,0.7,0.4, RGB_TABLE=rgb, $ TICKNAME=tickname, BORDER=1)9）;离散色条与一个额外的标签相比;颜色数量与这些标签在线下方。;与填充轮廓图中使用的类似IDL tickname = STRING(10*INDGEN(11), FORMAT=(I0)IDL cb = COLORBAR(POSITION=0.05,0.25,0.7,0.3, R