Processing/BrightnessTemperature: get_atm_temp.pro

File get_atm_temp.pro, 4.1 KB (added by mggr, 8 years ago)

Top level routine

Line 
1pro get_atm_temp,level1b_file,tempfile,emissivity=emissivity,brightness=brightness
2
3; Converts ATM radiances to temperature and outputs it in a copy of the ATM
4; level 1B file with band 11 replaced by scaled temperature.
5; Also replaces band 1 (useless anyway!) with a sinusoidal weighting function
6; from zero at the two ends to 65535 in the middle of the scan.
7; If emissivity is not set, a value of 0.98 typical of water is used.
8; If brightness is set, the emissivity is set to 1 so brightness temperature is
9; output.
10
11scale=100.
12lambda=8000.+findgen(141)*50.
13response=[.355,.601,1.516,3.861,9.228,18.909,32.217,49.088,66.147,74.362, $
14  73.921,74.65,79.395,85.655,89.796,89.961,86.917,83.729,81.554,81.689,82.382, $
15  83.753,85.398,87.673,90.818,93.771,97.239,99.093,99.788,99.303,97.256, $
16  94.367,91.656,89.055,89.012,89.984,91.436,93.076,96.083,97.547,97.839, $
17  96.276,95.595,93.101,90.414,89.376,88.846,88.937,89.45,89.09,88.989,89.273, $
18  89.493,91.365,91.361,90.564,93.651,96.481,94.359,89.362,83.905,79.661, $
19  77.323,73.96,71.288,68.293,63.318,59.337,56.775,54.588,52.397,50.389,48.655, $
20  46.431,44.638,42.024,41.344,39.604,39.065,38.557,38.788,38.903,39.965, $
21  40.435,41.115,41.589,40.722,40.448,39.724,39.7,40.509,39.513,39.043,38.721, $
22  38.377,37.527,36.191,35.059,35.889,32.06,30.039,27.925,25.09,23.672,20.754, $
23  18.254,16.791,15.736,12.703,10.61,8.807,8.775,7.364,6.668,6.233,4.96,3.874, $
24  3.444,3.832,4.096,3.517,3.584,1.962,2.043,2.2,2.027,2.063,1.373,2.123,1.357, $
25  1.182,1.296,.652,1.234,.645,.676,.695,.687,.577,.686,.757]
26if n_elements(emissivity) gt 0 then eps=emissivity[0] else $
27  if keyword_set(brightness) then eps=1. else $
28  eps=.98 ; water emissivity, could make it angle dependent
29if ~hdf_ishdf(level1b_file) then begin
30  ; let's try an ENVI BIP - look for the header
31  if ~file_test(level1b_file+'.hdr') then $
32    message,'Unidentified level 1B file - not HDF or BIP'
33  openr,unit,level1b_file+'.hdr',/get
34  str=''
35  repeat readf,unit,str until strcmp(strmid(str,0,7),'samples',/fold)
36  ns=long(strmid(str,strpos(str,'=')+1))
37  repeat readf,unit,str until strcmp(strmid(str,0,5),'lines',/fold)
38  nl=long(strmid(str,strpos(str,'=')+1))
39  repeat readf,unit,str until strcmp(strmid(str,0,5),'bands',/fold)
40  nb=long(strmid(str,strpos(str,'=')+1))
41  if ~max(nb eq [1,11]) then message,'Need 1 or 11 bands in ENVI file'
42  repeat readf,unit,str until strcmp(strmid(str,0,10),'interleave',/fold)
43  format=strlowcase(strtrim(strmid(str,strpos(str,'=')+1),2))
44  if format ne 'bip' then $
45    message,'Level 1B file format '+format+' not currently supported'
46  bip=1
47  close,unit
48  openr,unit,level1b_file
49  if nb eq 11 then indata=uintarr(11,ns) else dn=uintarr(ns)
50  if ~file_test(tempfile) then begin
51    temp_sdid=hdf_sd_start(tempfile,/create)
52    temp_sdsid=hdf_sd_create(temp_sdid,'ATdata',[ns,nl,11],/dfnt_uint16)
53    hdf_sd_endaccess,temp_sdsid
54    hdf_sd_end,temp_sdid
55  endif
56endif else begin
57  bip=0
58  if ~file_test(tempfile) then file_copy,level1b_file,tempfile
59  level1b_fid=hdf_open(level1b_file)
60  level1b_sdid=hdf_sd_start(level1b_file)
61  index=hdf_sd_nametoindex(level1b_sdid,'ATdata')
62  level1b_sdsid=hdf_sd_select(level1b_sdid,index)
63  hdf_sd_getinfo,level1b_sdsid,dims=dims
64  ns=dims[0]
65  nl=dims[1]
66endelse
67file_chmod,tempfile,/u_write
68weights=uint(sin(dindgen(ns)*!dpi/(ns-1))*65535d0+.5d0)
69temp_sdid=hdf_sd_start(tempfile,/rdwr)
70index=hdf_sd_nametoindex(temp_sdid,'ATdata')
71temp_sdsid=hdf_sd_select(temp_sdid,index)
72for line=0,nl-1 do begin
73  if bip then begin
74    if nb eq 11 then begin
75      readu,unit,indata
76      dn=reform(indata[10,*])
77    endif else readu,unit,dn
78  endif else hdf_sd_getdata,level1b_sdsid,dn,start=[0,line,10],count=[ns,1,1]
79  radiance=dn*.001
80  temp=t_bb(radiance,lambda,response,emissivity=eps,/trunc,t_bb_data=t_bb_data)
81  dn[*]=temp*scale>0<65535
82  hdf_sd_adddata,temp_sdsid,dn,start=[0,line,10],count=[ns,1,1]
83  hdf_sd_adddata,temp_sdsid,weights,start=[0,line,0],count=[ns,1,1]
84  check=check_math()
85  if (check and not 32) ne 0 then stop
86endfor
87hdf_sd_endaccess,temp_sdsid
88hdf_sd_end,temp_sdid
89if bip then free_lun,unit else begin
90  hdf_sd_endaccess,level1b_sdsid
91  hdf_sd_end,level1b_sdid
92  hdf_close,level1b_fid
93endelse
94
95end
96