IDL实现火灾监测与火点提取

---

基于白天的MODIS数据提取森林火点信息。MODIS数据为6个波段的ENVI格式文件，分别为第1、2、7波段表观反射率以及第21、31、32波段表观辐亮度数据，需要从MOD021KM或者MYD021KM数据中做波段提取与合成。

---

火点提取流程图.jpg

代码及过程：

主过程文件：

Pro Zhushan_fire_detection
 

  fn=dialog_pickfile(title='选择MODIS数据',get_path=work_dir)
  cd,work_dir
  envi_open_file,fn,r_fid=fid
  envi_file_query,fid,ns=ns,nl=nl,nb=nb,dims=dims,$
    data_type=data_type,interleave=interleave,offset=offset
  map_info=envi_get_map_info(fid=fid)

  ;read 1 2 7 bands
  ;rad21=envi_get_data(fid=fid,dims=dims,pos=0)
  ;rad21=envi_get_data(fid=fid,dims=dims,pos=0)
  ;rad21=envi_get_data(fid=fid,dims=dims,pos=0)
  ;rad21=envi_get_data(fid=fid,dims=dims,pos=0)

  ;rad31=envi_get_data(fid=fid,dims=dims,pos=1)
  ;rad32=envi_get_data(fid=fid,dims=dims,pos=2)
  ref1=envi_get_data(fid=fid,dims=dims,pos=0)
  ref2=envi_get_data(fid=fid,dims=dims,pos=1)
  ref7=envi_get_data(fid=fid,dims=dims,pos=2)
  ;read 21 31 32 bands
  rad21=envi_get_data(fid=fid,dims=dims,pos=3)
  rad31=envi_get_data(fid=fid,dims=dims,pos=4)
  rad32=envi_get_data(fid=fid,dims=dims,pos=5)

  ;brightness temporature
  wv=[3.99157,11.0121,12.0159]
  tcs=[0.9998646,0.995608,0.9997256]
  tci=[0.09262664,0.1302699,0.07181833]

  tb21=cal_tb(temporary(rad21),wv[0],tcs[0],tci[0])
  tb31=cal_tb(temporary(rad31),wv[1],tcs[1],tci[1])
  tb32=cal_tb(temporary(rad32),wv[2],tcs[2],tci[2])


  ;Cloud detect
  cloud=ref1+ref2 gt 0.9 or tb32 lt 265 $
    or (ref1+ref2 gt 0.7 and tb32 lt 285)
  ;Water detect
  NDVI=(ref2-ref1)/(ref2+ref1)
  water=ref2 lt 0.15 and ref7 lt 0.05 and ndvi lt 0
  ;Water
  mask=cloud eq 0 and water eq 0
  fire=bytarr(ns,nl)

  ;extract fire
  w=where((Tb21 gt 310 or (Tb21 gt 300 and Tb21-tb31 gt 20))$
    and mask eq 1,count)
  if count gt 0 then fire[w]=1

  pfire_locations=where(Tb21 gt 310 and Tb21-Tb31 gt 10 and $
    ref2 lt 0.3 and mask eq 1 and fire eq 0,count)

  if count gt 0 then begin
    window_size=7
    for i=0,count-1 do begin
      fire_flg=fire_dect(Tb21,Tb31,mask,pfire_locations[I],$
        window_size)
      if fire_flg eq 1 then fire[pfire_locations[i]]=1
    endfor
  endif

  ;save fire point to *.evf'
  fire_location=where(fire eq 1,count)
  if count gt 0 then begin

    print,'The Number of Fire Points:',count

    xf=fire_location mod ns
    yf=fire_location/ns

    envi_convert_file_coordinates,fid,xf,yf,xmap,ymap,/to_map


    ;creat '.evf'
    o_fn=dialog_pickfile(title='Save the Result of Fire Points Monitoring')
    evf_ptr=envi_evf_define_init(o_fn+'.evf',data_type=4,$
      projection=map_info.proj,layer_name='Detected fires')

    ;add fire points to evf file
    pts_fires=[transpose(xmap),transpose(ymap)]
    for i=0,count-1 do begin
      envi_evf_define_add_record,evf_ptr,pts_fires[*,I]
    endfor

    ;save and close
    evf_id=envi_evf_define_close(evf_ptr,/return_id)
    envi_evf_close,evf_id

    ;edit atributes of *.evf file
    attributes=replicate({Name:'Fire',ID:0L},count)

    for i=0,count-1 do attributes[i].id=I+1
    envi_write_dbf_file,o_fn+'.dbf',attributes

  endif else begin

    print,'Did not Detect the Fire Points'

  endelse

end

火点检测函数文件：

function fire_dect,Tb21,Tb31,mask,pfire_location,window_size
  ;fire points or not
  sz=size(mask)
  ns=sz[1] & nl=sz[2]
  
  ;coordinate transport
  colu=pfire_location mod ns;column
  row =pfire_location/ns    ;row
  ;the beginning and ending of background windows
  col_s=colu-window_size/2 > 0
  col_e=colu+window_size/2 < ns-1
  row_s=row-window_size/2 > 0
  row_e=row+window_size/2 < ns-1

  Tb21_pfire=Tb21[colu,row]
  Tb31_pfire=Tb31[colu,row]
  dif_Tb_pfire=Tb21_pfire-Tb31_pfire

  tTb21=Tb21[col_s:col_e,row_s:row_e]
  tTb31=Tb31[col_s:col_e,row_s:row_e]
  dif_Tb=tTb21-tTb31
  tmask=mask[col_s:col_e,row_s:row_e]

  mask_bfire=Tb21 gt 325 and dif_Tb gt 20
  w_bfire=where(mask_bfire eq 1,nbfire)
  w_background=where(tmask eq 1 and mask_bfire eq 0,nbackground)

  if nbackground gt 0 then begin

    mean_Tb21=mean(tTb21[w_background])
    mean_Tb31=mean(tTb31[w_background])
    mean_dif_Tb=mean(dif_Tb[w_background])

    MAD_Tb21=mean(abs(tTb21[w_background]-mean_Tb21))
    MAD_Tb31=mean(abs(tTb31[w_background]-mean_Tb31))
    MAD_dif_Tb=mean(abs(dif_Tb[w_background]-mean_dif_Tb))
    if nbfire gt 0 then begin
      MAD_Tb21_bfire=mean(abs(dif_Tb[w_bfire]-mean(tTb21[w_bfire])))
    endif else begin
      MAD_Tb21_bfire=0

    endelse

    if dif_Tb_pfire gt (mean_dif_Tb+3.5*MAD_dif_Tb)and $
      dif_Tb_pfire gt (mean_dif_Tb+6) and $
      Tb21_pfire gt (mean_Tb21+3*MAD_Tb21)and $
      (Tb31_pfire gt (mean_Tb31+MAD_Tb31-4) or MAD_Tb21_bfire gt 5)$
      eq 1 then begin
      return,1
    endif else begin
      return,0
    endelse

  endif else begin
    return,0
  endelse
end

亮温计算函数文件：

;coculate the lightness temperature of thermal 
function cal_tb,radiance,wv,tcs,tci

  C1=1.1910659e8
  C2=1.438833e4
  x1= wv*alog(1+c1/(wv^5*radiance))
  result=(C2/x1 - tci ) / tcs
  return,result
end

CSDN分享

简单书写，

希望你十分美好！