IDL 程序设计笔记.docx

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