hdf_sd_attr.pro

文件是基于IDL开发环境的读取EOS的HDF文件格式的代码

function hdf_sd_varread,hdfid,varname,attname,_extra = extra_keywords

;- Check arguments
if (n_params() ne 3) then $
  message, 'Usage: HDF_SD_VARREAD, HDFID, VARNAME, DATA'
if (n_elements(hdfid) eq 0) then $
  message, 'Argument HDFID is undefined'
if (n_elements(varname) eq 0) then $
  message, 'Argument VARNAME is undefined'
;if (arg_present(data) eq 0) then $
;  message, 'Argument DATA cannot be modified'

;- Get index of the requested variable

index = hdf_sd_nametoindex(hdfid, varname)

if (index lt 0) then $
  message, 'SDS was not found: ' + varname

;- Select and read the SDS

varid = hdf_sd_select(hdfid, index)

hdf_sd_getdata, varid, data, _extra=extra_keywords

sizeofdata = size(data)

;window,/free,xsize = sizeofdata(1),ysize = sizeofdata(2),title = varname
;
;tvscl,data

attindex = hdf_sd_attrfind(varid,attname);寻找属性

hdf_sd_attrinfo,varid,attindex,count=c,name = n,data =d,type = t

;help,c,n,d,t


hdf_sd_endaccess,varid


return,data

end

;返回文件的属性有哪些
function hdf_sd_vardir,hdfid

;- Check arguments
if (n_params() ne 1) then $
  message, 'Usage: RESULT = HDF_SD_VARDIR(HDFID)'
if (n_elements(hdfid) eq 0) then $
  message, 'HDFID is undefined'

;- Set default return value
varnames = ''

;- Get file information
hdf_sd_fileinfo, hdfid, nvars, ngatts

;- If variables were found, get variable names
if (nvars gt 0) then begin
  varnames = strarr(nvars)
  for index = 0L, nvars - 1L do begin
    varid = hdf_sd_select(hdfid, index)
    hdf_sd_getinfo, varid, name=name
    hdf_sd_endaccess, varid
    varnames[index] = name
  endfor
endif

;- Return the result
return, varnames

END


function equalvalues,array

result = 0

meanvalue = mean(array)
std_div = STDDEV(array)

if ( (std_div eq 0) and (meanvalue eq array(2,2)) ) then begin
      result = 1
 endif else begin
     result = 0
 endelse

return,result

end


;文件匹配，其中covertype为分类文件数据，emissdata为3d数组(6,240,240),
;anglevalue为day and night角度匹配后所得到的文件
;函数将范返回一个三维数组比辐射率,为(6,17,131),其中131为角度的变化范围
function datamatch,covertype,emissdata,anglevalue

result = fltarr(6,17,131);result为比辐射率返回值，三维数组

value = fltarr(6,17,131);value保存各波段各种地物类型各个角度的比辐射率的总和
count = fltarr(6,17,131);count保存各波段各种地物类型各个角度的比辐射率的数目

result(*,*,*) = 0
value(*,*,*) = 0
count(*,*,*) = 0


for kcur = 0,5 do begin

         for nx = 0,239 do begin
            for ny = 0,239 do begin

              if(anglevalue(nx,ny) le 130) then begin;第一重if选择，用于控制是否具有有效角度值

                      if((covertype(nx,ny)lt 20 )and (emissdata(kcur,nx,ny) ne 0 ))then  begin

                         value(kcur,covertype(nx,ny),anglevalue(nx,ny)) = value(kcur,covertype(nx,ny),anglevalue(nx,ny))+emissdata(kcur,nx,ny)
                         count(kcur,covertype(nx,ny),anglevalue(nx,ny)) ++

                      endif

              endif

      		endfor
		 endfor



	for i =0,16 do begin;第一重for控制类型
   		for j = 0,130 do begin

  			if(count(kcur,i,j)ne 0) then begin
 				result(kcur,i,j) = (0.002*value(kcur,i,j)/count(kcur,i,j))+0.49
			endif else begin
			  result(kcur,i,j) = 0
		  endelse
   		endfor
	endfor

endfor;控制大循环




return,result


end




pro hdf_sd_attr_rhz

;*****************************************************************************
;*****************************************************************************
;               处理covertype  from 1200*1200 to 240*240
;*****************************************************************************
;*****************************************************************************

cd,'J:\'


hdfid_covertpye = hdf_sd_start('MOD12Q1.A2001001.h26v05.004.2004149205632.hdf')

coverfile_att = hdf_sd_vardir(hdfid_covertpye)


data_covertype = hdf_sd_varread(hdfid_covertpye,coverfile_att(0),'valid_range')



covertype = bytarr(240,240)

covertypesize = size(data_covertype)


for ni = 2,covertypesize(1)-3,5 do begin
   for nj =  2,covertypesize(2)-3,5 do begin

      mask = [[data_covertype(ni-2:ni+2,nj-2)],[data_covertype(ni-2:ni+2,nj-1)],[data_covertype(ni-2:ni+2,nj)],$
              [data_covertype(ni-2:ni+2,nj+1)],[data_covertype(ni-2:ni+2,nj+2)]]

     result = equalvalues(mask)

     if(result eq 1) then begin
         covertype(ni/5,nj/5) = data_covertype(ni,nj)
      endif  else begin

          covertype(ni/5,nj/5) = 20

    endelse


 endfor
endfor



;window,0,xsize = 240,ysize = 240,title = 'covertype_proceed'
;tvscl,covertype


;***************************************************************************
;***************************************************************************
;                     Obtain the data of  HDF
;***************************************************************************
;***************************************************************************

;list1 = dialog_pickfile(filter='*.hdf')

cd,'J:\h26v052000'

list = findfile('*.hdf')

list_size = size(list)


name = ['Emis_20','Emis_22','Emis_23','Emis_29','Emis_31','Emis_32']

openw,lun,'newh26v06.txt',/get_lun

for filenum = 0,list_size(1)-1 do begin;list_size(1)-100 do begin


      		hdfid = hdf_sd_start(list(filenum))



			varname = "
			attname = "
			data = bytarr(6,240,240)



			;申请角度变量
			dayview_angle = intarr(240,240)
			nightview_angle = intarr(240,240)
			anglevalue = intarr(240,240)

			attribute = hdf_sd_vardir(hdfid)

			attname = 'long_name';辅助变量，无实际意义

			  printf,lun,list(filenum);写入文件名

			 ;获得白天和夜晚的观测角度
			  dayview_angle = hdf_sd_varread(hdfid,attribute(3),attname)
			  nightview_angle =hdf_sd_varread(hdfid,attribute(7),attname)

				;拟合白天与夜晚的观测角度

					for x = 0,239 do begin
				 		 for y =0,239 do begin

				      		day = dayview_angle(x,y);
				      		night = nightview_angle(x,y)

				;角度属性中，可利用值的范围为0-130，而不可用值都用255表示

				   			 if((day le 130) and (night gt 130) ) then anglevalue(x,y) = day

				   			 if((night le 130) and (day gt 130) ) then anglevalue(x,y) = night

				;以下两句话处理当二者都不为255的情况
    						if((day le 130)and (night le 130) and day le night)then anglevalue(x,y) = day
   						    if((day le 130)and (night le 130) and day gt night)then anglevalue(x,y) = night

				;处理当二者都比有效值130大时赋予255
    						if((day gt 130) and (night gt 130) )then anglevalue(x,y) = 255
     					endfor
  					endfor

     		for kcur = 0,5 do begin;循环一个文件的6个比辐射率

         		 	varname = attribute(kcur+8)
        		  	data[kcur,*,*] = hdf_sd_varread(hdfid,varname,attname)

	        endfor

;*************************************************************************
;                     文件匹配
;*************************************************************************

 			emiss = fltarr(6,17,131)
			emiss =  datamatch(covertype,data,anglevalue);进入文件匹配过程

 				for ki = 0,5 do begin
 				       printf,lun,name(ki)
                    for kj = 0,16 do begin
                    ;打印的格式为列为地物17种，行为角度的比辐射率131种
						printf,lun,emiss(ki,kj,*),format = '(131f10)';;5i5, 第一个5表示每行有5个元素，第2个5表示每两个数之间的距离;5f5.2 表示为double
			     	endfor
				endfor

			hdf_sd_end, hdfid

endfor

free_lun,lun

end