-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathm3_photom_i_over_f.pro
179 lines (179 loc) · 5.53 KB
/
m3_photom_i_over_f.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
pro m3_photom_i_over_f, event, infile, obsfile, outfile, r_fid=r_fid, flag=flag, sd=sd
;
; Calculates I/F as:
;
; I/F = (radiance * pi)/solar_spectrum
;
; where solar_spectrum is divided by moon-sun distance (in AU) squared.
;
;---------------------------------------------------------
compile_opt strictarr
;
; simple (standard) error catching mechanism
;
CATCH, error
IF (error NE 0) THEN BEGIN
catch, /cancel
help, ns, nl, nb
ok = m3_error_message()
return
ENDIF
;
; set a flag for failure or success of this routine
;
flag=0
;
; Get input file
;
if n_elements(infile) eq 1 then begin
envi_open_file, infile, r_fid=rad_fid
endif else begin
envi_select, fid=rad_fid, title='Select Radiance File', /no_spec, /no_dims, /file_only
endelse
if rad_fid eq -1 then return
envi_file_query, rad_fid, ns=ns, nl=nl, nb=nb, data_type=dt, wl=wl, xstart=xstart, $
ystart=ystart, fname=infile
rad_pos = lindgen(nb)
inherit=envi_set_inheritance(rad_fid, [-1L,0,ns-1,0,nl-1], /spatial)
;
; get output file if not passed in
;
if n_elements(outfile) eq 0 then begin
tlb = widget_auto_base(title = 'M3 Level 2')
ofw = widget_outf(tlb, uvalue='outf', func='m3_file_test', /auto)
result = auto_wid_mng(tlb)
if result.accept eq 0 then return
outfile = result.outf
endif
;
; is the input file global or target?
;
case m3_get_imaging_mode(infile) of
'T':mode='target'
'G':mode='global'
else: message, 'Could not determine if input data is global or target data.'
endcase
;
; open the solar spectrum
;
solar_spec_file = filepath(root_dir=m3_programrootdir(), subdir = ['resources', 'level2'], $
'solar_spectrum_' + mode + '.sli')
envi_open_file, solar_spec_file, r_fid=ss_fid, /no_realize
envi_file_query, ss_fid, ns=ns_ss, nl=nl_ss, nb=nb_ss, wl=wl_ss, spec_names=ss_name
nb=n_elements(wl_ss)
;
; read the solar spectrum, reform to the size of a BIL tile
;
sli_dims = [-1,0,ns_ss - 1,0,nl_ss - 1]
solar_spectrum = double(rebin(transpose(envi_get_data(fid=ss_fid, dims=sli_dims, pos=0)), ns, nb))
envi_file_mng, id=ss_fid, /remove
;
; determine solar distance if not passed in
;
if n_elements(sd) eq 0 then begin
;
; Try to find the OBS file
;
if n_elements(obsfile) eq 1 then begin
envi_open_file, obsfile, r_fid=obs_fid
endif else begin
indir = file_dirname(infile, /mark_directory)
complete_basename = file_basename(infile, '.IMG')
basename = strmid(complete_basename, 0, strpos(complete_basename, "_", /reverse_search) + 1)
obsfile_test = indir + basename + 'OBS.IMG'
if file_test(obsfile_test) then begin
obsfile = obsfile_test
envi_open_file, obsfile, r_fid=obs_fid
endif else begin
title='Select OBS File (OR CANCEL TO USE 1AU FOR SOLAR DISTANCE)'
envi_select, fid=obs_fid, title=title, /no_spec, /no_dims, /file_only
endelse
endelse
if obs_fid eq -1 then solar_distance_squared=1.0d
;
; get the solar distance from the band name of the solar distance band in the OBS file,
; convert it to double, then square it
;
if obs_fid gt 0 then begin
envi_file_query, obs_fid, bnames=obs_bnames
pos=5 ;where(stregex(obs_bnames, 'To-Sun Path Length', /boolean),count)
;if count eq 0 then message, 'Did not find the To-Sun Path Length band'
idx=stregex(obs_bnames[pos], '\(au-(.+)\)', length=len, /subexpr)
solar_distance_squared=double((STRMID(obs_bnames[pos], idx, len))[1])^2 ;square sd here
if solar_distance_squared eq 0 then begin
msg='Solar distance not found. Use 1.0 AU?'
yesno=dialog_message(msg, /question)
if yesno eq 'Yes' then solar_distance_squared = 1.0d else return
endif
endif
endif else solar_distance_squared=double(sd)
;
; divide out the solar distance from the solar spectrum tile
;
sd_norm_ss = solar_spectrum/solar_distance_squared
;
; set up progress bar
;
rstr = ['Processing input file: ' + file_basename(infile), 'This may take several minutes...']
envi_report_init, rstr, title="M3 I/F", base=repbase, /interupt
;
; open output file and setup tiling
;
rad_tile=envi_init_tile(rad_fid, rad_pos, num_tiles=num_tiles, interleave=1) ; BIL tiles
envi_report_inc, repbase, num_tiles
openw, lun, outfile, /get_lun
;
; tile loop
;
for i = 0L, num_tiles-1 do begin
envi_report_stat, repbase, i, num_tiles, cancel = cancel
if cancel then begin
envi_report_init, base=repbase, /finish
return
endif
;
; get radiance tile
;
radiance = double(envi_get_tile(rad_tile, i))
;
; calculate I/F - now with PI!
;
i_over_f = (radiance*!dpi)/sd_norm_ss
;
; convert to float and write to file
;
writeu, lun, float(i_over_f)
;
endfor
;
; tiling clean up
;
envi_tile_done, rad_tile
envi_report_init, base=repbase, /finish
;
; close file, setup envi header
;
free_lun,lun
descrip='M3 I/F Cube'
bnames='I/F Band ' + strtrim(indgen(nb) + 1,2)
envi_setup_head, fname=outfile, ns=ns, nl=nl, nb=nb, offset=0, interleave=1, $
data_type=4, inherit=inherit, descrip=descrip, wl=wl_ss, xstart=xstart, $
ystart=ystart, bnames=bnames, /write
envi_open_file, outfile, r_fid=r_fid
;
; add reflectance version and solar spectrum version to envi header
;
envi_assign_header_value, fid=r_fid, keyword='m3_IoverF', value='A'
envi_assign_header_value, fid=r_fid, keyword='m3_solar_spectrum_used', value=ss_name
;
; close all the files
;
envi_file_mng, id=rad_fid, /remove
envi_file_mng, id=obs_fid, /remove
envi_file_mng, id=r_fid, /remove
;
; set flag to success
;
flag=1
;
end