2023年IDL程序设计笔记.docx

目录
Part one：文件读写操作 1
1、格式化输入与输出（read\print） 1
2、ASCII文件读写 2
3、二进制文件读写 5
4、图像格式文件读写 7
5、科学数据格式读写*** 8
Part two：图形绘制 11
（一）plot过程绘制 11
（二）plotg()函数绘制 12
（三）散点图 12
（四）柱状图、条形图 15
（五）绘制颜色条colorbar 18
Part three: ENVI二次开发 20
1、自定义波段运算函数 20
2、IDL中常见波段运算 21
3、ENVI classic 创建自定义菜单 22
4、ENVI classic+IDL 二次开发 23

Part one：文献读写操作
1、格式化输入与输出（read\print）
1) Read:
该函数用于从键盘键入数据，默认数据类型为浮点型，若要输入其他类型数据，需要先定义数据类型。
Eg:

2) Print：
该函数用于将数据输出打印到控制台，这里重要讲格式化控制输出format关键字。
表01 常用格式化控制符
格式代码
输出效果
aN
字符或者字符串按照N个字符宽度输出
fn.m
浮点型数组按照N个字符宽度输出，小数点保存M位
dn.m
double型数组按照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,/append

4） 读写操作：
IDL中运用readf函数读取文本文献、运用printf函数写入ASCII文献
① Readf, lun, var01, var02….
② printf, lun, var01, var02….
5） 文献关闭：
在对文献操作完毕之后，需要关闭文献的逻辑设备号，eg: freee_lun，lun

EG:

*如何获取某个文献你的列数:
先读取一行数据，然后运用空格进行字符串拆分；最后运用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:\zengsk\IDL\Cha06\data\FY2C_TBB_IR1_OTG_20231130_AOAD.AWX'
openr, lun, file, /get_lun
;跳过前20个字节
point_lun,lun,20
Headline = indgen(3);定义三个整型数据，分别读取记录长度、头文献记录数、数据记录数
readu, lun, Headline
print,headline[0], headline[1],headline[2]

;定位到第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(Headline[2],Headline[0])
;定位到数据部分，其头文献结束的字节位置为“头文献记录数*记录长度”
point_lun, lun, Headline[0] * Headline[1]
print,Headline[0] * Headline[1]
;读取数据部分
readu, lun, data
help,data
;显示数据
window, 0, xsize = 500, ysize = 500
tv, congrid(data, 500, 500)
end
4、图像格式文献读写
IDL自带了丰富的多种图像读写函数，如BMP，JPG，PNG, JPEG2023， 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、对遥感图像而言，其数据存放方式有BSQ\BIP\BIL，以一个三维m列n行的遥感数据为例：BSQ的表达方式为[m, n, 3]、BIP表达方式为[3, m ,n]、BIL的表达方式为[m, 3, n]，常用的图像文献一般是BIP方式存储，遥感文献一般是BSQ。

1） JPEG格式：读read_jpeg\写write_jpeg
Eg01:

Eg02:
2） BMP格式：(不采用任何压缩)，数据存储量一般很大
读：read_imge(fname) or read_bmp(fname, /rgb)
写：write_bmp, fname, data, /rgb
3） 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:\ENVI\IDL下的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:\ENVI\IDL\HJ-1.tif', Bandata_BIP
end

2） 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,value
end

4）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、 坐标轴设立：
x\ytitle设立标题；x\ycharsize设立字体大小；
x\yrange用于设立坐标轴范围，注意加上x\ystyle=1来强行设定
坐标轴范围
x\yticks和x\yminor分别设立主刻度和最小刻度间隔数目
x\ytickname设立刻度的名称
4、 颜色设立：
color: 设立图像颜色；一般用十六进制来表达颜色；’FFFFFF’xl表达白色；’000000’xl为黑色
background：设立背景颜色
5、 添加标注： xyouts过程用于在图像窗口中添加标注信息。
示例：xyouts, 0.18, 0.80, ‘RMSE=’, color=’FFFFFF’xl, charsize=1.2, /normal
6、 图形保存为文献：一般方法为运用TVRD函数拷贝图形窗口的内容，在写入图像文献。
img = tvrd(x, y, channel=value, true=1\2\3, /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 = data[0,*];观测数据
y = data[1,*];估算数据

plot, x, y, psym=2,xrange=[5,15],yrange=[5,15],$
xtitle='Observed value', ytitle='Estimated value', $
background='FFFFFF'xl, color='oooooo'xl, $
xstyle=1, ystyle=1, $
charsize=1.3, /nodata
oplot, x, y, psym=5, color='000000'xl
;绘制基准线
x2=[5, 15] & y2=[5, 15]
oplot, x2, y2, color='FF0000'xl

;***********计算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='000000'xl, charsize=1.2, /data
xyouts, 6, 12.3, RMSE_label, color='000000'xl, charsize=1.2, /data

;保存散点图为图像文献，格式为.png
img=tvrd(true=1)

;法一：
o_fn=dialog_pickfile(title='文献保存为：');键入的文献名要包含后缀名
write_png, o_fn, img
;法二
;write_png, 'E:\ENVI\scatter.png', img
end

Eg02: 读取文本文献中的数据并运用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=data[0, *]
Es=data[1, *]
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_fn
end
（四）柱状图、条形图
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
area[i]=strmid(var[i], 0, 8)
number[i]=fix(strmid(var[i], 11, 3))
endfor
help, area
print, area
help, number
print, number[3]

;##################画柱状图####
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='EEEEEE'xl, fill_color='AOAOFF'xl, $
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
area[i]=strmid(var[i], 0, 8)
number[i]=fix(strmid(var[i], 11, 3))
endfor
help, area
print, area
help, number
print, number[3]
;##################画柱状图################
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='EEEEEE'xl, fill_color='AOAOFF'xl, $
margin=0.2, /horizontal)
;############文献保存为jpg#####################
o_fn=dialog_pickfile(title='图片保存为：')+'.jpg';注意后面加一个后缀名字符串更加方便
His.save, o_fn
End
运营结果：

（五）绘制颜色条colorbar
Eg:
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=[5001000],TICKVALUES=[550,650,750,850,950], TICKNAME=['A','B','C','D','E']

6) ;创建一个colorbar,设立垂直色条——orientation=1，范围0-100——range，锥形端——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 = rgb
IDL> rgb = rgb[25 + 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], RGB_TABLE=rgb, $
> TICKNAME=tickname, BORDER=1)

10）;带渐缩端的离散颜色条。;与填充轮廓图中使用的类似
IDL> cb = COLORBAR(POSITION=[0.05,0.15,0.7,0.2], RGB_TABLE=rgb, $
> TICKNAME=tickname, BORDER=1, TAPER=1)

11）;创建一个离散的颜色条，并且在左端（最小值处）渐缩——taper=2。
IDL> cb = COLORBAR(POSITION=[0.05,0.05,0.45,0.1], RGB_TABLE=rgb,
> TICKNAME=tickname, BORDER=1, TAPER=2)

;创建一个离散的颜色条，并且在右端端（最大值处）渐缩——taper=3。

Part three: ENVI二次开发
1、自定义波段运算函数
编写波段运算函数跟在envi中直接输入数学表达式相比，在操作上更灵活，函数中可以添加更多的控制，如运算中的错误检测等。波段运算函数的编写非常简朴，接受输入波段，解决数据，并返回结果。以下面的方式定义：

Function bm_func, b1, [b2,…,bn, parameters and keywords]
processing steps!
return, result
end
编写好function之后，点击编译，在控制台输入：
IDL> save, filename='E:\ENVI\cal_ndvi.sav',/routine
然后把.sav文献拷贝到envi的安装目录之下。

波谱运算：
波谱运算的操作环节跟波段运算基本一致，不同之处有：
1) 变量命名以“S”或“s”开头；
2) 执行时 ENVI 主菜单下，选择[Basic Tools]→[Spectral Math]

2、IDL中常见波段运算

*注：在IDL中NaN值如何得到：

3、ENVI classic 创建自定义菜单
;自动创建ENVI菜单过程
; ENVI_DEFINE_MENU_BUTTON, ButtonInfo [, /DISPLAY], $ EVENT_PRO=string | /MENU | UVALUE=string [, POSITION=long integer or