IDL 程序设计笔记.docx
《IDL 程序设计笔记.docx》由会员分享,可在线阅读,更多相关《IDL 程序设计笔记.docx(31页珍藏版)》请在冰豆网上搜索。
IDL程序设计笔记
Partone:
文件读写操作
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卫星数据
proread_AWX
file='E:
\zengsk\IDL\Cha06\data\FY2C_TBB_IR1_OTG_20061130_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,JPEG2000,DICOM,TIFF等
相关函数:
文件查询:
query_image(fname,dimensions=dimensions,channels=nb);channels返回波段数、dimensions返回图像的行列数
读取图像:
read_image(fname);结果返回图像数组
写入图像:
write_image,fname(可带路径),‘jpg/bmp/tiff’,data,/order;order关键字用于设置图像的纵坐标从上往下算起,默认为从下往上。
*注:
、以上三个函数对于IDL支持的图像格式都适用。
、显示图像时:
tv用法:
tv,data,true=1/2/3;其中关键字true=1表示数据格式为(3,m,n)、true=2表示数组格式为(m,3,n)、true=3表示格式为(m,n,3)
、tvrd:
屏幕拷贝函数,返回当前直接图形窗口或设备的指定矩形部分的内容。
、对遥感图像而言,其数据存放方式有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)orread_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的反射率数据集
protest_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数据格式的例子
proread_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(networkCommonDataForm,网络通用数据格式),常用于气象科学数据存储。
;下面介绍一个读取ECMWF数据的例子,ECMWF为*.nc格式数据
;
protest_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格式读取:
Parttwo:
图形绘制
(一)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绘制散点图
proplot_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='Observedvalue',ytitle='Estimatedvalue',$
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函数绘制散点图
proplot_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='ObservedValue',ytitle='EstimatedValue',$;再坐标轴
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功能函数绘制柱状图
protest_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)
fori=0,nl-1dobegin
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=ar