IDL开发(一)——水体营养指数

IDL开发——水体营养指数

---

前言

利用IDL实现水质评价，包括：叶绿素浓度的反演、悬浮物浓度、总磷浓度、总氮浓度、高猛酸盐浓度的反演，将反演后的叶绿素 a 浓度、高锰酸盐浓度、总磷浓度、总氮浓度和透明度 5 个指标进行营养状态指数计算。
参考文献：谭小琴,罗勇,赵铮,甘先霞,张洪吉,陈建华.基于高分遥感的河流水质反演研究——以金马河温江段为例[J].环境生态学,2020,2(07):29-36.

---

一、实现功能

输入叶绿素浓度、悬浮物浓度文件，自动计算总磷浓度、总氮浓度、高猛酸盐浓度的反演，将反演后的叶绿素 a 浓度、高锰酸盐浓度、总磷浓度、总氮浓度和透明度 5 个指标进行营养状态指数计算。最后按照营养状态分级显示水质评价结果

二、IDL代码

1.主窗体

在主窗体中写明各base成分，标注event。

wButton3 = WIDGET_BUTTON(wtemp3, value = '水质评价', event_pro = 'TLI')

2.TLI event事件代码

pro TLI,ev
  COMPILE_OPT IDL2
  ENVI,/RESTORE_BASE_SAVE_FILES
  ENVI_Batch_init

  filename=ENVI_PICKFILE(title='请输入叶绿素a浓度文件',filter='*.tif')
  ENVI_OPEN_FILE,filename,r_fid=fid,/no_realize
  ENVI_FILE_QUERY,fid,ns = numSample,nl = numLine,nb = numBand,$
    data_type=dataType,dims = dims
    
  filename2=ENVI_PICKFILE(title='请输入悬浮物浓度文件',filter='*.tif')
  ENVI_OPEN_FILE,filename2,r_fid=fid2,/no_realize
  ENVI_FILE_QUERY,fid2,ns = numSample,nl = numLine,nb = numBand,$
    data_type=dataType,dims = dims2

  mapInfo = ENVI_GET_MAP_INFO(FID=fid)
  imageData = FLTARR(numSample,numLine,numBand)
  
  imageData[*,*,0] = ENVI_GET_DATA(fid = fid,dims = dims,pos = 0)
  imageData_SS = FLTARR(numSample,numLine,numBand)
  imageData_SS[*,*,0] = ENVI_GET_DATA(fid = fid2,dims = dims,pos = 0)
  ;****************************************计算各个指数*******************
  TP1= FLTARR(numSample,numLine,1)
 ; imageData[*,*,0] = ENVI_GET_DATA(fid = fid,dims = dims,pos = 0)
  TP1 = (-0.00078)*float(imageData[*,*,0]*1.0)+0.0417
  ;*********************************************************************
  TN1= FLTARR(numSample,numLine,1)
  ;imageData[*,*,0] = ENVI_GET_DATA(fid = fid,dims = dims,pos = 0)
  TN1 = 3.166-0.03479*float(imageData[*,*,0]*1.0)
  ;*********************************************************************
  CODmn1= FLTARR(numSample,numLine,1)
  ;imageData[*,*,0] = ENVI_GET_DATA(fid = fid,dims = dims,pos = 0)
  CODmn1 = 0.050*float(imageData[*,*,0]*1.0)+4.543
  ;*********************************************************************
  TLI_Chla= FLTARR(numSample,numLine,1)
  TLI_CODMn= FLTARR(numSample,numLine,1)
  TLI_TP= FLTARR(numSample,numLine,1)
  TLI_TN= FLTARR(numSample,numLine,1)
  TLI_SD= FLTARR(numSample,numLine,1)
  TLI1=FLTARR(numSample,numLine,1)
  TLI_PLOT=BYTARR(numSample,numLine,1)
  TLI_Chla= 10*(2.5+1.086*alog(float(imageData[*,*,0]*1.0)))
  TLI_CODMn= 10*(2.5+1.086*alog(float(CODmn1)))
  TLI_TP= 10*(9.436+1.624*alog(float(TP1)))
  TLI_TN= 10*(5.453+1.694*alog(float(TN1)))
  TLI_SD= 10*(5.118-1.94*alog(imageData_SS[*,*,0]*1.0))
  
  TLI1=0.2662*FLOAT(TLI_Chla)+0.1878*FLOAT(TLI_CODMn)+0.1790*FLOAT(TLI_TP)+0.1834*FLOAT(TLI_TN)+0.1834*FLOAT(TLI_SD)

  max_v=max(TLI1)
  print,max_v
  dst=[[0.0,30.0,50,50,75],$
    [30.0,50,175,175,50],$
    [50.0,60.0,150,200,50],$
    [60.0,70.0,75,200,25],$
    [70.0,max_v,0,175,0]]
    
    
    
   ;**************************************************密度分割****************************************************
   
    
   sz2=size(dst)
   n_ds=sz2[2]

   sz=size(TLI1)
   ns=sz[1]&nl=sz[2]
   result_r=bytarr(ns,nl)
   result_g=bytarr(ns,nl)
   result_b=bytarr(ns,nl)

   ;密度分割
   for i=0,n_ds-1 do begin
     w=where(TLI1 gt dst[0,i] and TLI1 le dst[1,i],count)
     if count gt 0 then begin
       result_r[w]=dst[2,i]
       result_g[w]=dst[3,i]
       result_b[w]=dst[4,i]
     endif
   endfor

   tli_outmap=bytarr(3,ns,nl)
   tli_outmap[0,*,*]=result_r
   tli_outmap[1,*,*]=result_g
   tli_outmap[2,*,*]=result_b

    ;tli_outmap=density_slice(TLI1,dst)
    sz3=size(tli_outmap)
    ns3=sz3[2] & nl3=sz3[3]
 
  W1=window(window_title='水体营养状态分级',dimension=[numSample,numLine*1.25])
  ;i1=image(TLI_PLOT,/current,dimension=[numSample,numLine],$
  i1=image(tli_outmap,dimension=[ns3,nl3],margin=0,window_title='TLI Level',position=[0,0.2,1,1],/order)
   ;tickname=[dst[0,0],transpose(dst[1,*])]
   ;tickname=string(tickname,format='(f6.2)') 
   tickname=['Level 1:0','Level 2:30','Level 3:50','Level 4:60','Level 5:70']
   tab_dst=dst[2:4,*]
  c1=colorbar(position=[0.05,0.1,0.95,0.15],title='TLI',$
   rgb_table=tab_dst,tickname=tickname,taper=0,/border)

  DELVAR,imageData

  o_file=dialog_pickfile(title='图像文件保存为')
  write_image,o_file,'tiff',TLI1
  
end