IDL个人复习笔记

IDL-TAE

• 输出Hello IDL字样

(1)输出字样

pro HelloIDL

  tmp=dialog_message('HelloIDL',/information)

end

(2)彩色条显示

pro color_map

  data=[0:255]

  help,data

  data=rebin(data,256,100)

  window,1,xsize=280,ysize=110

  loadct,13

  tv,data

end

(3)彩色图显示

Pro setting_color

device,decomposed=0

file=filepath('r_seeberi.jpg',$

  subdirectory=['examples','data'])

read_jpeg,file,img,/grayscale

help,img

window,1,xsize=280,ysize=195;,title='loadct'

tv,img

window,2,xsize=280,ysize=195

loadct,13

tv,img

end 

(4) 读取ACSII

pro readasciifile

  asciiFile = dialog_pickfile(title='输入ascii.txt',filter='*.txt',path=file_dirname(routine_filepath('readasciifile')))

openr,lun,asciifile,/get_lun

if lun eq -1 then begin

  void=dialog_message('文件错误!',/error)

  return

endif

;逐行读取

;tmp="

;while(~eof(lun)) do begin

;  readf,lun,tmp

;  print,tmp

;endwhile

;基于数据类型读取

  tmp=intarr(3)

  str=strarr(2)

  data=fltarr(2,4)

  openr,lun,asciifile,/get_lun

  readf,lun,tem

  readf,lun,str

  readf,lun,data

  print,data

  FREE_LUN,lun

end

二、利用IDL语言，提取数据的采集时间和文件后缀

(1)信息提取

pro timeandfile

  file='C:\Users\Administrator\Desktop\geocoded\SI_20210226_1_disp_frst_est_geo_ql.tif'

  a=strmid(file,43,8)

  print,a

  b=strmid(file,74,4)

  print,b

  viod=dialog_message('print successful please check',/information)

end

(2)时间信息提取

pro fetch_information

file='E:\IDLWorkspace\Default\202012254\data\LC08_L1TP_124036_20130410_20170505_01_T1_MTL.txt'

   openr,lun,file,/get_lun

   tmp=''

  str='DATE';更改此处

  ;str1='    DATE_ACQUIRED = 2020-10-17'

  ;print,strmid(str1,4,4)

  while(~eof(lun))do begin

readf,lun,tmp

tmp1=strtrim(tmp,1)

    ;  help,tmp

    str_tmp=strmid(tmp1,0,4)

    if strmatch(str_tmp,str) eq 1 then begin

      print,tmp

    endif

    ;  print,tmp

  endwhile

  viod=dialog_message('fetch successful please check',/information)

end

(3)Tiff裁剪

pro tiffcut

  tiffile='E:\IDLWorkspace\Default\202012254\data\color_half_world.tif'

  data=read_image(tiffile)

  arr=read_tiff(tiffile)

  help,arr

  arr=read_tiff(tiffile,sub_rect=[99,99,100,100])

  write_image,'E:\IDLWorkspace\Default\202012254\data\Rcut.tif','tiff',arr

  viod=dialog_message('Cut successful please check Data folder',/information)

end

三、以风云二号C卫星（FY-2C）为例读取信息并且输出

(1)图像显示

pro fengyunxianshi

file='E:\IDLWorkspace\Default\202012254\data\FY2C_TBB_IR1_OTG_20061130_AOAD.AWX'

  openr,file_lun,file,/get_lun

  help,file

  point_lun,file_lun,20

  headline=indgen(3)

  readu,file_lun,headline

  data=bytarr(headline[2],headline[0])

  point_lun,file_lun,headline[0]*headline[1]

  readu,file_lun,data

  window,0,xsize=600,ysize=600

  tv,congrid(data,600,600)

end

(2)卫星名读取

pro fenyunname

file='E:\IDLWorkspace\Default\202012254\data\FY2C_TBB_IR1_OTG_20061130_AOAD.AWX'

  openr,file_lun,file,/get_lun

  help,file

  point_lun,file_lun,40

  name=bytarr(8)

  readu,file_lun,name

  name=string(name)

  viod=dialog_message('Read successful',/information)

  print,name

end

四、微波光学遥感影像读取

(1)SAR数据

pro weibo_cc_duqu

  file=' E:\IDLWorkspace\Default\021107_030116_cc'

  openr,lun,file,/get_lun

  data=make_array(1001,800,type=4)

  readu,lun,data

  free_lun,lun

  tvscl,data,/order

end

（2）光学

pro guangxue_img_duqu

  file='E:\Software\ENVIsoftware\ENVI51\classic\data\can_tmr.img'

  openr,lun,file,/get_lun

  data=make_array(640,400,6,type=1)

  readu,lun,data

  free_lun,lun

  blue=data[*,*,0]

  ;help,blue

  tvscl,blue,/order

end

五、利用IDL语言进行多图绘制

(1) Obeverd-theory

pro kexuehuitu1

  theory=sin(2.0*findgen(200)*!pi/25.0)*exp(-0.02*findgen(200))

  obeverd=theory+randomu(seed,201)*0.4-0.2

  p1=plot(obeverd,/overplot,'c2',name='obsvered',color='MEDIUM BLUE')

  p2=plot(theory,/overplot,'r3',name='theory',title='CGJ-polar',color='light blue')

  l=legend(target=[p1,p2],position=[140,0.9],/data)

end

(2) Sinewave- cosine

Pro kexuehuitu2

sinewave=sin(2.0*findgen(200))*exp(-0.02*findgen(200))

cosine=cos(2.0*findgen(200))*exp(-0.02*findgen(200))

p=plot(sinewave,'3r',AXIS_STYLE=1,POSITION=[.075,.075,.90,.90],title='CGJ-Plot')

p=plot(cosine,'b',AXIS_STYLE=1,/CURRENT,POSITION=[.60,.60,.90,.90])

end

(3)

pro kexuehuitu3 

theory = SIN(2.0*FINDGEN(200)*!PI/25.0)*EXP(-0.02*FINDGEN(200))

observed=theory+randomu(seed,201)*0.4-0.2

  p=plot(theory,'-3b',axis_style=1,position=[.075,.075,.90,.90],title='CGJ-Plot')

  p=plot(observed,'-1r',axis_style=1,/current,position=[.60,.60,.90,.90])

end

(4)分离四个窗口绘制

Pro kexuehuitu4

  data=findgen(100)

  !p.multi=[4,2,2,0,0]

  p=plot(data,title='polar1')

  p=plot(sin(data/12),cos(data/12),title='polar2')

  p=plot(data,data,/polar,title='polar3')

  p=plot(sin(data/12),xtitle='x',ytitle='y',title='polar4')

  !p.multi=0

end

(5)综合一个窗口绘制

pro kexuehuitu5

  oribkcolor=!p.BACKGROUND

  !p.BACKGROUND=!p.COLOR

  !p.COLOR=oribkcolor

  data=findgen(100)

  window,0,xsize=800,ysize=600

  !p.multi=[6,2,2,0,0]

    plot,data,title='polar1'

    plot,sin(data/12),cos(data/12),title='polar2'

    plot,data,data,/polar,title='polar3'

    plot,sin(data/12),psym=4,xtitle='x',ytitle='y',title='polar4'

  !p.multi=0

end

(6)DEM_Tiff绘制显示

pro dem_xianshi

  file='E:\IDLWorkspace\Default\202012254\data\PDSTIF.tif'

  dem=read_binary(file,data_dims=[64,64])

  c1=contour(dem,rgb_table=30,/fill,planar=0,title='202012254陈国军')

  cbar=colorbar(target=c1,orientation=1,position=[0.07,0.2,0.1,0.75])

end

六、IDL实现地图显示，绘制河南省行政边界和平顶山市行政边界，显示河流

(1)

pro map_xianshi

  ant_map=MAP('STEREOGRAPHIC',     $

    ;limit=[33,113,34,114],        $

    limit=[31,110,37,117],         $

    center_latitude=33,            $

    center_longitude=113,          $

    title='Henan Pingdingshan Map',$

    fill_color ='light blue')

  mc = mapcontinents(/countries,fill_color='ivory')

  rivers = mapcontinents(/rivers,color='navy')

  file = 'E:\IDLWorkspace\Default\202012254\data\Shp\河南省.shp'

  henan = mapcontinents(file,fill_color='aqua',color='black')

  file_1='E:\IDLWorkspace\Default\202012254\data\Shp\平顶山地区.shp'

  pingdingshan=mapcontinents(file_1,fill_color='misty rose',color='black',thick='2')

  rivers = mapcontinents(/rivers,color='blue',thick='2')

end

(2)

pro map_view

ant_map=map('stereographic',limit=[10,80,50,130],center_latitude=30,center_longitude=105,title='map example',fill_color='light blue')

  mc=mapcontinents(/countries,fill_color='LAVENDER BLUSH')

  rivers=mapcontinents(/rivers,color='c',thick='2')

  file='E:\IDLWorkspace\Default\202012254\data\Shp\河南省.shp';添加路径

  beijing=mapcontinents(file,color='b',thick='2')

end

(3)

pro map_guojie

  ant_map=map('stereographic',$

    limit=[10,80,50,130],     $

    center_latitude=30,       $

    center_longitude=105,     $

    title='Map Example1',     $

    fill_color='light blue')

mc=mapcontinents(/countries,fill_color='beige')

rivers=mapcontinents(/rivers,color='cyan')

file='D:\3program bags\ENVI 5.3\ENVI53\data\natural_earth_vectors\countries.shp'

lushan=mapcontinents(file,color='red')

end

七、度分秒转换。

（1）

pro dufenmiao_convert

  for i=0,19 do begin

    y=['33:30:30.963920','33:48:34.952265','33:48:34.934094','33:48,34.915104','33:48,34.887888',$

'33:48:34.860366','33:48:34.834770','33:48:34.801785','33:47:40.993152','33:47:40.676262',$

'33:47:40.375581','33:47:40.080723','33:47:38.315148','33:47:38.606982','33:47:38.906466',$

'33:47:39.207777','33:47:39.505308','33:47:39.806259','33:47:40.107705','33:47:40.402122']

    x=y[i]

    a=float(strmid(x,0,2))

    b=float(strmid(x,3,2))

    c=float(strmid(x,6,9))

    num=a+b/60+c/3600

    print,num

z=['113:17:19.563873','113:17:15.669447','113:17:11.781438','113:17:07.887462','113:17:03.963768',$

'113:17:00.027015','113:16:56.136387','113:16:52.050585','113:17:05.532684','113:17:05.779527',$

'113:17:06.014463','113:17:06.245475','113:17:07.678428','113:17:07.441710','113:17:07.206432',$

'113:17:06.970227','113:17:06.740709','113:17:06.509004','113:17:06.276372','113:17:06.045063']

    n=z[i]

    a1=double(strmid(n,0,3));float也可，这里使用双精度效果更贴近于原始数据

    b1=double(strmid(n,4,2))

    c1=double(strmid(n,7,9))

    num1=a1+b1/60+c1/3600

    print,num1

    jieguo=[num,num1]

    openw,lun,'E:\IDLWorkspace\Default\calresult.txt',/get_lun

    printf,lun,jieguo

    free_lun,lun

  endfor

end

（2）pro dufenmiao_recover

  for i=0,8 do begin

y=['103,35,11','111,20,45','104,05,11','102,34,43','102,33,21','111,20,45','104,05,11','102,34,43','102,33,21']

    x=y[i]

    a=strmid(x,0,3)

    b=strmid(x,4,2)

    c=strmid(x,7,2)

    d=float(a)

    e=float(b)/60

    f=float(c)/3600

    num=d+e+f

    print,num

  endfor

end

（3）

pro dufenmiao_convert1

  file='E:\IDLWorkspace\Default\JinWeidu.txt'

  openr,lun,file,/get_lun

  range=dblarr(1,20)

  readf,lun,range

  print,range

  free_lun,lun

  i=0

  n=dblarr(1,20)

  while i le 20 do begin

    m=range[*,i]

    a=m[0]

    b=m[1]

    c=m[2]

    d=a

    e=b/60

    f=c/3600

    num=d+e+f

    n[i]=num

    i=i+1

    print,num

    print,'--**--'

  endwhile

  print,n

end

八、编写IDL程序，利用二元二次多项式模型，计算下列图像坐标和地面坐标之间的变换系数。

（1）

pro jisuanxishu

  c='E:\IDLWorkspace\Default\Test.txt'

  openr,lun,c,/get_lun

  range=dblarr(4,15)

  readf,lun,range

  Lx=range[0,*];像点x坐标

  Ly=range[1,*];像点y坐标

  Ax=range[2,*];地面点x坐标

  Ay=range[3,*];地面点y坐标

  A=dblarr(6,15);系数矩阵

  a_=dblarr(1,7);变换系数

  b_=dblarr(1,7);变换系数

  free_lun,lun

  ;给系数矩阵赋值

  A[0,*]=1

  A[1,*]=Ax

  A[2,*]=Ay

  A[3,*]=Ax*Ay

  A[4,*]=Ax*Ax

  A[5,*]=Ay*Ay

  ;数组运算

  LxT=reverse(Lx,1)

  LyT=transpose(Ly)

  AT=transpose(A);A转置

  AT_A=A#AT;计算AT*A

  AT_A_=INVERT(AT_A);计算AT_A的逆

  a_1=AT#AT_A_

  a_=LxT#a_1;计算变换系数a

  b_=LyT#a_1;计算变换系数b

  print,a_,b_

end

九、几何矫正

pro jihejiaozheng

  file='E:\IDLWorkspace\Default\202012254\data\几何校正.txt'

  openr,lun,file,/get_lun

  skip_lun,lun,1,/lines

  arr=fltarr(4,16)

  readf,lun,arr

  ;print,arr

  u=arr[0,*]

  v=arr[1,*]

  x=arr[2,*]

  y=arr[3,*]

  a=transpose([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1])

  u_v_all=u##transpose(v)

  u_v =diag_matrix(u_v_all)

  u_u_all=u##transpose(u)

  u_u=diag_matrix(u_u_all)

  v_v_all=v##transpose(v)

  v_v=diag_matrix(v_v_all)

  u_v_n=[a,u,v,transpose(u_v),transpose(u_u),transpose(v_v)]

  a_n=invert(transpose(u_v_n)##u_v_n)##transpose(u_v_n)##x

  b_n=invert(transpose(u_v_n)##u_v_n)##transpose(u_v_n)##y

  xishu=[a_n,b_n]

  openw,lun,'E:\IDLWorkspace\Default\result.txt',/get_lun

  printf,lun,xishu

  free_lun,lun

end

附加简单的练习程序：

1、数组与数组运算

pro z_shuzu_zufu

;数组的运算

arr1=indgen(5)

print,arr1

arr2=arr1+3

print,arr2

arr3=indgen(4)

print,arr3

arr4=indgen(5)

print,arr4

sum1=arr3+arr4

print,sum1

arr4=[1,3,5,7]

print,arr4

arr5=[2,4,6,8]

print,arr5

print,arr4+arr5

;数组的合并

a=indgen(2,5)

b=indgen(4,5)

c=[a,b]

help,c

d=indgen(2,3)

e=[[a],[d]]

help,e

;矩阵的求逆

arr=indgen(5,5)

invert_arr=invert(arr)

arr1=long(invert_arr)

help,arr1

print,arr1

;字符串反转

a='Xiao Chen is a little sleepy.'

print,a

b=byte(a)

C=reverse(b)

d=string(c)

print,d

;文件名截取

file='E:\IDLWorkspace\Default\202012254\data\FY2C_TBB_IR1_OTG_20061130_AOAD.AWX'

basename=file_basename(file)

print,strmid(basename,17,8)

;度分秒取整转换

lat=45.3438

IDL-TAE

• 输出Hello IDL字样

(1)输出字样

pro HelloIDL

  tmp=dialog_message('HelloIDL',/information)

end

(2)彩色条显示

pro color_map

  data=[0:255]

  help,data

  data=rebin(data,256,100)

  window,1,xsize=280,ysize=110

  loadct,13

  tv,data

end

(3)彩色图显示

Pro setting_color

device,decomposed=0

file=filepath('r_seeberi.jpg',$

  subdirectory=['examples','data'])

read_jpeg,file,img,/grayscale

help,img

window,1,xsize=280,ysize=195;,title='loadct'

tv,img

window,2,xsize=280,ysize=195

loadct,13

tv,img

end 

(4) 读取ACSII

pro readasciifile

  asciiFile = dialog_pickfile(title='输入ascii.txt',filter='*.txt',path=file_dirname(routine_filepath('readasciifile')))

openr,lun,asciifile,/get_lun

if lun eq -1 then begin

  void=dialog_message('文件错误!',/error)

  return

endif

;逐行读取

;tmp="

;while(~eof(lun)) do begin

;  readf,lun,tmp

;  print,tmp

;endwhile

;基于数据类型读取

  tmp=intarr(3)

  str=strarr(2)

  data=fltarr(2,4)

  openr,lun,asciifile,/get_lun

  readf,lun,tem

  readf,lun,str

  readf,lun,data

  print,data

  FREE_LUN,lun

end

二、利用IDL语言，提取数据的采集时间和文件后缀

(1)信息提取

pro timeandfile

  file='C:\Users\Administrator\Desktop\geocoded\SI_20210226_1_disp_frst_est_geo_ql.tif'

  a=strmid(file,43,8)

  print,a

  b=strmid(file,74,4)

  print,b

  viod=dialog_message('print successful please check',/information)

end

(2)时间信息提取

pro fetch_information

file='E:\IDLWorkspace\Default\202012254\data\LC08_L1TP_124036_20130410_20170505_01_T1_MTL.txt'

   openr,lun,file,/get_lun

   tmp=''

  str='DATE';更改此处

  ;str1='    DATE_ACQUIRED = 2020-10-17'

  ;print,strmid(str1,4,4)

  while(~eof(lun))do begin

readf,lun,tmp

tmp1=strtrim(tmp,1)

    ;  help,tmp

    str_tmp=strmid(tmp1,0,4)

    if strmatch(str_tmp,str) eq 1 then begin

      print,tmp

    endif

    ;  print,tmp

  endwhile

  viod=dialog_message('fetch successful please check',/information)

end

(3)Tiff裁剪

pro tiffcut

  tiffile='E:\IDLWorkspace\Default\202012254\data\color_half_world.tif'

  data=read_image(tiffile)

  arr=read_tiff(tiffile)

  help,arr

  arr=read_tiff(tiffile,sub_rect=[99,99,100,100])

  write_image,'E:\IDLWorkspace\Default\202012254\data\Rcut.tif','tiff',arr

  viod=dialog_message('Cut successful please check Data folder',/information)

end

三、以风云二号C卫星（FY-2C）为例读取信息并且输出

(1)图像显示

pro fengyunxianshi

file='E:\IDLWorkspace\Default\202012254\data\FY2C_TBB_IR1_OTG_20061130_AOAD.AWX'

  openr,file_lun,file,/get_lun

  help,file

  point_lun,file_lun,20

  headline=indgen(3)

  readu,file_lun,headline

  data=bytarr(headline[2],headline[0])

  point_lun,file_lun,headline[0]*headline[1]

  readu,file_lun,data

  window,0,xsize=600,ysize=600

  tv,congrid(data,600,600)

end

(2)卫星名读取

pro fenyunname

file='E:\IDLWorkspace\Default\202012254\data\FY2C_TBB_IR1_OTG_20061130_AOAD.AWX'

  openr,file_lun,file,/get_lun

  help,file

  point_lun,file_lun,40

  name=bytarr(8)

  readu,file_lun,name

  name=string(name)

  viod=dialog_message('Read successful',/information)

  print,name

end

四、微波光学遥感影像读取

(1)SAR数据

pro weibo_cc_duqu

  file=' E:\IDLWorkspace\Default\021107_030116_cc'

  openr,lun,file,/get_lun

  data=make_array(1001,800,type=4)

  readu,lun,data

  free_lun,lun

  tvscl,data,/order

end

（2）光学

pro guangxue_img_duqu

  file='E:\Software\ENVIsoftware\ENVI51\classic\data\can_tmr.img'

  openr,lun,file,/get_lun

  data=make_array(640,400,6,type=1)

  readu,lun,data

  free_lun,lun

  blue=data[*,*,0]

  ;help,blue

  tvscl,blue,/order

end

五、利用IDL语言进行多图绘制

(1) Obeverd-theory

pro kexuehuitu1

  theory=sin(2.0*findgen(200)*!pi/25.0)*exp(-0.02*findgen(200))

  obeverd=theory+randomu(seed,201)*0.4-0.2

  p1=plot(obeverd,/overplot,'c2',name='obsvered',color='MEDIUM BLUE')

  p2=plot(theory,/overplot,'r3',name='theory',title='CGJ-polar',color='light blue')

  l=legend(target=[p1,p2],position=[140,0.9],/data)

end

(2) Sinewave- cosine

Pro kexuehuitu2

sinewave=sin(2.0*findgen(200))*exp(-0.02*findgen(200))

cosine=cos(2.0*findgen(200))*exp(-0.02*findgen(200))

p=plot(sinewave,'3r',AXIS_STYLE=1,POSITION=[.075,.075,.90,.90],title='CGJ-Plot')

p=plot(cosine,'b',AXIS_STYLE=1,/CURRENT,POSITION=[.60,.60,.90,.90])

end

(3)

pro kexuehuitu3 

theory = SIN(2.0*FINDGEN(200)*!PI/25.0)*EXP(-0.02*FINDGEN(200))

observed=theory+randomu(seed,201)*0.4-0.2

  p=plot(theory,'-3b',axis_style=1,position=[.075,.075,.90,.90],title='CGJ-Plot')

  p=plot(observed,'-1r',axis_style=1,/current,position=[.60,.60,.90,.90])

end

(4)分离四个窗口绘制

Pro kexuehuitu4

  data=findgen(100)

  !p.multi=[4,2,2,0,0]

  p=plot(data,title='polar1')

  p=plot(sin(data/12),cos(data/12),title='polar2')

  p=plot(data,data,/polar,title='polar3')

  p=plot(sin(data/12),xtitle='x',ytitle='y',title='polar4')

  !p.multi=0

end

(5)综合一个窗口绘制

pro kexuehuitu5

  oribkcolor=!p.BACKGROUND

  !p.BACKGROUND=!p.COLOR

  !p.COLOR=oribkcolor

  data=findgen(100)

  window,0,xsize=800,ysize=600

  !p.multi=[6,2,2,0,0]

    plot,data,title='polar1'

    plot,sin(data/12),cos(data/12),title='polar2'

    plot,data,data,/polar,title='polar3'

    plot,sin(data/12),psym=4,xtitle='x',ytitle='y',title='polar4'

  !p.multi=0

end

(6)DEM_Tiff绘制显示

pro dem_xianshi

  file='E:\IDLWorkspace\Default\202012254\data\PDSTIF.tif'

  dem=read_binary(file,data_dims=[64,64])

  c1=contour(dem,rgb_table=30,/fill,planar=0,title='202012254陈国军')

  cbar=colorbar(target=c1,orientation=1,position=[0.07,0.2,0.1,0.75])

end

六、IDL实现地图显示，绘制河南省行政边界和平顶山市行政边界，显示河流

(1)

pro map_xianshi

  ant_map=MAP('STEREOGRAPHIC',     $

    ;limit=[33,113,34,114],        $

    limit=[31,110,37,117],         $

    center_latitude=33,            $

    center_longitude=113,          $

    title='Henan Pingdingshan Map',$

    fill_color ='light blue')

  mc = mapcontinents(/countries,fill_color='ivory')

  rivers = mapcontinents(/rivers,color='navy')

  file = 'E:\IDLWorkspace\Default\202012254\data\Shp\河南省.shp'

  henan = mapcontinents(file,fill_color='aqua',color='black')

  file_1='E:\IDLWorkspace\Default\202012254\data\Shp\平顶山地区.shp'

  pingdingshan=mapcontinents(file_1,fill_color='misty rose',color='black',thick='2')

  rivers = mapcontinents(/rivers,color='blue',thick='2')

end

(2)

pro map_view

ant_map=map('stereographic',limit=[10,80,50,130],center_latitude=30,center_longitude=105,title='map example',fill_color='light blue')

  mc=mapcontinents(/countries,fill_color='LAVENDER BLUSH')

  rivers=mapcontinents(/rivers,color='c',thick='2')

  file='E:\IDLWorkspace\Default\202012254\data\Shp\河南省.shp';添加路径

  beijing=mapcontinents(file,color='b',thick='2')

end

(3)

pro map_guojie

  ant_map=map('stereographic',$

    limit=[10,80,50,130],     $

    center_latitude=30,       $

    center_longitude=105,     $

    title='Map Example1',     $

    fill_color='light blue')

mc=mapcontinents(/countries,fill_color='beige')

rivers=mapcontinents(/rivers,color='cyan')

file='D:\3program bags\ENVI 5.3\ENVI53\data\natural_earth_vectors\countries.shp'

lushan=mapcontinents(file,color='red')

end

七、度分秒转换。

（1）

pro dufenmiao_convert

  for i=0,19 do begin

    y=['33:30:30.963920','33:48:34.952265','33:48:34.934094','33:48,34.915104','33:48,34.887888',$

'33:48:34.860366','33:48:34.834770','33:48:34.801785','33:47:40.993152','33:47:40.676262',$

'33:47:40.375581','33:47:40.080723','33:47:38.315148','33:47:38.606982','33:47:38.906466',$

'33:47:39.207777','33:47:39.505308','33:47:39.806259','33:47:40.107705','33:47:40.402122']

    x=y[i]

    a=float(strmid(x,0,2))

    b=float(strmid(x,3,2))

    c=float(strmid(x,6,9))

    num=a+b/60+c/3600

    print,num

z=['113:17:19.563873','113:17:15.669447','113:17:11.781438','113:17:07.887462','113:17:03.963768',$

'113:17:00.027015','113:16:56.136387','113:16:52.050585','113:17:05.532684','113:17:05.779527',$

'113:17:06.014463','113:17:06.245475','113:17:07.678428','113:17:07.441710','113:17:07.206432',$

'113:17:06.970227','113:17:06.740709','113:17:06.509004','113:17:06.276372','113:17:06.045063']

    n=z[i]

    a1=double(strmid(n,0,3));float也可，这里使用双精度效果更贴近于原始数据

    b1=double(strmid(n,4,2))

    c1=double(strmid(n,7,9))

    num1=a1+b1/60+c1/3600

    print,num1

    jieguo=[num,num1]

    openw,lun,'E:\IDLWorkspace\Default\calresult.txt',/get_lun

    printf,lun,jieguo

    free_lun,lun

  endfor

end

（2）pro dufenmiao_recover

  for i=0,8 do begin

y=['103,35,11','111,20,45','104,05,11','102,34,43','102,33,21','111,20,45','104,05,11','102,34,43','102,33,21']

    x=y[i]

    a=strmid(x,0,3)

    b=strmid(x,4,2)

    c=strmid(x,7,2)

    d=float(a)

    e=float(b)/60

    f=float(c)/3600

    num=d+e+f

    print,num

  endfor

end

（3）

pro dufenmiao_convert1

  file='E:\IDLWorkspace\Default\JinWeidu.txt'

  openr,lun,file,/get_lun

  range=dblarr(1,20)

  readf,lun,range

  print,range

  free_lun,lun

  i=0

  n=dblarr(1,20)

  while i le 20 do begin

    m=range[*,i]

    a=m[0]

    b=m[1]

    c=m[2]

    d=a

    e=b/60

    f=c/3600

    num=d+e+f

    n[i]=num

    i=i+1

    print,num

    print,'--**--'

  endwhile

  print,n

end

八、编写IDL程序，利用二元二次多项式模型，计算下列图像坐标和地面坐标之间的变换系数。

（1）

pro jisuanxishu

  c='E:\IDLWorkspace\Default\Test.txt'

  openr,lun,c,/get_lun

  range=dblarr(4,15)

  readf,lun,range

  Lx=range[0,*];像点x坐标

  Ly=range[1,*];像点y坐标

  Ax=range[2,*];地面点x坐标

  Ay=range[3,*];地面点y坐标

  A=dblarr(6,15);系数矩阵

  a_=dblarr(1,7);变换系数

  b_=dblarr(1,7);变换系数

  free_lun,lun

  ;给系数矩阵赋值

  A[0,*]=1

  A[1,*]=Ax

  A[2,*]=Ay

  A[3,*]=Ax*Ay

  A[4,*]=Ax*Ax

  A[5,*]=Ay*Ay

  ;数组运算

  LxT=reverse(Lx,1)

  LyT=transpose(Ly)

  AT=transpose(A);A转置

  AT_A=A#AT;计算AT*A

  AT_A_=INVERT(AT_A);计算AT_A的逆

  a_1=AT#AT_A_

  a_=LxT#a_1;计算变换系数a

  b_=LyT#a_1;计算变换系数b

  print,a_,b_

end

九、几何矫正

pro jihejiaozheng

  file='E:\IDLWorkspace\Default\202012254\data\几何校正.txt'

  openr,lun,file,/get_lun

  skip_lun,lun,1,/lines

  arr=fltarr(4,16)

  readf,lun,arr

  ;print,arr

  u=arr[0,*]

  v=arr[1,*]

  x=arr[2,*]

  y=arr[3,*]

  a=transpose([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1])

  u_v_all=u##transpose(v)

  u_v =diag_matrix(u_v_all)

  u_u_all=u##transpose(u)

  u_u=diag_matrix(u_u_all)

  v_v_all=v##transpose(v)

  v_v=diag_matrix(v_v_all)

  u_v_n=[a,u,v,transpose(u_v),transpose(u_u),transpose(v_v)]

  a_n=invert(transpose(u_v_n)##u_v_n)##transpose(u_v_n)##x

  b_n=invert(transpose(u_v_n)##u_v_n)##transpose(u_v_n)##y

  xishu=[a_n,b_n]

  openw,lun,'E:\IDLWorkspace\Default\result.txt',/get_lun

  printf,lun,xishu

  free_lun,lun

end

附加简单的练习程序：

1、数组与数组运算

pro z_shuzu_zufu

;数组的运算

arr1=indgen(5)

print,arr1

arr2=arr1+3

print,arr2

arr3=indgen(4)

print,arr3

arr4=indgen(5)

print,arr4

sum1=arr3+arr4

print,sum1

arr4=[1,3,5,7]

print,arr4

arr5=[2,4,6,8]

print,arr5

print,arr4+arr5

;数组的合并

a=indgen(2,5)

b=indgen(4,5)

c=[a,b]

help,c

d=indgen(2,3)

e=[[a],[d]]

help,e

;矩阵的求逆

arr=indgen(5,5)

invert_arr=invert(arr)

arr1=long(invert_arr)

help,arr1

print,arr1

;字符串反转

a='Xiao Chen is a little sleepy.'

print,a

b=byte(a)

C=reverse(b)

d=string(c)

print,d

;文件名截取

file='E:\IDLWorkspace\Default\202012254\data\FY2C_TBB_IR1_OTG_20061130_AOAD.AWX'

basename=file_basename(file)

print,strmid(basename,17,8)

;度分秒取整转换

lat=45.3438

a=strmid(lat,6,2)

m=(lat-a) *60

b=strmid(m,6,2)

n=round(m-b)*60

c=strmid(n,10,2)

print,string (a+'°' )+string(b+'′' )+string(c+'″')

end

a=strmid(lat,6,2)

m=(lat-a) *60

b=strmid(m,6,2)

n=round(m-b)*60

c=strmid(n,10,2)

print,string (a+'°' )+string(b+'′' )+string(c+'″')

end