叶面积指数遥感估算

pro LAICX
  ;打开文件
  cd,'E:\xxsx\idl\idl\sx9'
  fn_ndvi='NDVI'
  fn_landcover='Landcover'
  
  ;读取数据
  envi_open_file,fn_ndvi,r_fid=fid_ndvi
  envi_open_file,fn_landcover,r_fid=fid_landcover
  envi_file_query,fid_ndvi,ns=ns,nl=nl,nb=nb,dims=dims
  map_info=envi_get_map_info(fid=fid_ndvi)
  ndvi=envi_get_data(fid=fid_ndvi,dims=dims,pos=0);读取一条波段，此处第一条
  landcover=envi_get_data(fid=fid_landcover,dims=dims,pos=0)

  ;计算
  lai=laijs(ndvi,landcover)
  ;分别类型平均值
  ;res=mean(lai)
  ; print,'res:'+res
  ;保存结果
  envi_write_envi_file,lai,out_name='lairesult',map_info=map_info
  
  
  w1=where(landcover eq 2)
  w2=where(landcover eq 3)
  print,mean(lai[w1])
  print,mean(lai[w2])
  ;
  ;直方图以w1weibiaoz

  h1=histogram(lai[w1],min=0,max=10,binsize=0.1,locations=locations)
  h2=histogram(lai[w2],min=0,max=10,binsize=0.1)
  pl1=plot(locations,h1,xtitle='lai',ytitle='shumu',color=[0,0,255])
  pl2=plot(locations,h2,color=[255,0,255],/overplot)

 end

function laijs,ndvi,landcover
  ;创建ndvi的数组
  szn=size(ndvi)
  result=make_array(size=szn)

  ;根据landcover计算0
  ;1ndvi++++++0.125+++++0.825,制备1
  w=where(landcover eq 2 and ndvi ge 0.125 and ndvi le 0.825,count)
  if count gt 0 then result[w]=0.1836*exp(4.37*ndvi[w])
  ;ndvi+++++825
  w=where(landcover eq 2  and ndvi ge 0.825,count)
  if count gt 0 then result[w]=6.606
  zb1=mean(w)
  
  ;leixing2
  w=where(landcover eq 3 and ndvi ge 0.125 and ndvi le 0.825,count)
  if count gt 0 then result[w]=0.0884*exp(4.96*ndvi[w])
  ;ndvi+++++825
  w=where(landcover eq 3 and   ndvi ge 0.825,count)
  if count gt 0 then result[w]=6.091
 
  ;p1=plot(result,count,xtilte='lai',ytilte='shumu',color=[0,0,255])
  ;print,'zb2:'+zb2
  return,result
  
 end