forked from drdaphos/LIDAR
-
Notifications
You must be signed in to change notification settings - Fork 2
/
lid_horace_read_test.pro
executable file
·348 lines (263 loc) · 8.73 KB
/
lid_horace_read_test.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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
pro lid_horace_read, core=core, horace=horace, $
press_alt=press_alt, no_gin=no_gin, aimms=use_aimms, verbose=verbose
;
; default behaviour is try reading core and if it does not exist read horace
; if /horace, ignore core; if /core, ignore horace
;
@lid_settings.include
if (n_elements(flno) NE 1) then $
message, 'No flight selected. lid_flight_select must be called first.'
openw, lgf, logfln, /get_lun, /append
printf, lgf, '--> LID_HORACE_READ'
tmp = {generalflightinfo, $
tim:0.0D, alt:0.0D, lat:0.0D, lon:0.0D, ptc:0.0D, rll:0.0D, $
rhg:0.0D, ofn:0.0D, con:0.0D, dis:0.0D, osf:0.0D, rsf:0.0D, $
invalid:0B}
; bypass core/horace data if ground based
if (gnd_based) then begin
n = long(landing - takeoff + 1.0D)
hor = replicate({generalflightinfo}, n)
hor.tim = takeoff + dindgen(n)
hor.alt = gnd_alt
hor.lat = gnd_lat
hor.lon = gnd_lon
hor.con = cos(gnd_offzenith * !dpi / 180.0D)
hor.rhg = replicate(_hundef_, n)
info_str = string($
format='(%"Ground based - Altitude %dm")', gnd_alt)
printf, lgf, info_str
if (keyword_set(verbose)) then print, info_str
goto, skip
endif
; determine which file we are going to use
msg = ''
coretried = 0
if (~keyword_set(horace)) then begin
coretried = 1
corefile_search = string( $
format='(%"core_faam_%04d%02d%02d_v004_r?_%s_1hz.nc")', $
yy, mth, dd, strlowcase(flno_core))
pushd, core_path
corefile = file_search(corefile_search, count=cnt, /fold_case)
popd
if (cnt GT 0) then begin
corefile = corefile[cnt-1]
endif else begin
corefile = corefile_search
endelse
if (cnt LE 0 || $
~file_test(core_path + corefile, /read, /regular)) then begin
msg += corefile + ' '
if (keyword_set(core)) then begin
message, msg + 'not found or not readable'
endif else begin
horace = 1
endelse
endif
endif
if (keyword_set(horace)) then begin
horacefile = string(format='(%"horace_%04d_%02d_%02d.dat")',yy,mth,dd)
if (~file_test(horc_path + horacefile, /read, /regular)) then begin
msg += horacefile + ' '
message, msg + 'not found or not readable'
endif else if (coretried) then begin
info_str = msg + 'not found; using ' + horacefile
printf, lgf, info_str
message, info_str, /continue
endif
endif
if (~keyword_set(horace)) then begin
; option 1 - read 1 Hz core data file downloaded from Faam
info_str = 'Reading aircraft data from ' + corefile
printf, lgf, info_str
if (keyword_set(verbose)) then print, info_str
corefln = corefile
params = [612, 610, 611, 618, 617, 616, $
582, 580, 581, 562, 561, 560, 575, 578]
status = mrfread(core_path + corefln, params, data, flags, $
start=takeoff, stop=landing, time=tim2)
info_str = ''
if (keyword_set(no_gin)) then begin
info_str = 'No-GIN option selected. Will use backup GPS.'
data[0:5,*] = data[6:11,*]
endif else begin
for i=0, 5 do begin
index = where(flags[i,*] NE 0 $
OR ~finite(data[i,*]), cnt)
if (cnt GE 1) then begin
info_str = 'Some navigation data missing. ' $
+ 'Will try to fill gaps with backup GPS.'
data[i,index] = data[i+6, index]
endif
endfor
endelse
if (info_str NE '') then begin
printf, lgf, info_str
message, info_str, /continue
endif
if (keyword_set(press_alt)) then begin
data[0,*] = data[13,*]
info_str = 'Pressure altitude option selected instead of GPS.'
printf, lgf, info_str
message, info_str, /continue
endif
n = n_elements(tim2)
hor = replicate({generalflightinfo}, n)
hor.tim = double(tim2)
hor.alt = reform(double(data[0,*]))
hor.lat = reform(double(data[1,*]))
hor.lon = reform(double(data[2,*]))
hdg = reform(double(data[3,*]))
hor.ptc = reform(double(data[4,*]))
hor.rll = reform(double(data[5,*]))
hor.rhg = reform(double(data[12,*]))
index = where(flags[12,*] NE 0, cnt)
if (cnt GE 1) then hor[index].rhg = _hundef_
endif else begin
; option 2 - read horace data file produced by the on-board viewer lidardisplay
; data will be interpolated to 1 Hz
info_str = 'Reading aircraft data from ' + horacefile
printf, lgf, info_str
if (keyword_set(verbose)) then print, info_str
readhoracedata, dir=horc_path, takeoffdate=[dd,mth,yy], $
tim2, alt2press, alt2, lat2, lon2, ptc2, rll2, rhg2, $
/report_always
if (keyword_set(press_alt)) then begin
alt2 = alt2press
info_str = 'Pressure altitude option selected instead of GPS.'
printf, lgf, info_str
message, info_str, /continue
endif
idx = where(finite(tim2) AND finite(alt2), cnt)
t1 = floor(tim2[idx[0]])
t2 = ceil(tim2[idx[cnt-1]])
n = long(t2 - t1 + 1.0D)
hor = replicate({generalflightinfo}, n)
hor.tim = t1 + dindgen(n)
hor.alt = interpol(alt2[idx], tim2[idx], hor.tim)
hor.lat = interpol(lat2[idx], tim2[idx], hor.tim)
hor.lon = interpol(lon2[idx], tim2[idx], hor.tim)
hor.ptc = interpol(ptc2[idx], tim2[idx], hor.tim)
hor.rll = interpol(rll2[idx], tim2[idx], hor.tim)
hor.rhg = interpol(rhg2[idx], tim2[idx], hor.tim)
endelse
if (keyword_set(use_aimms)) then begin
info_str = 'Using AIMMS GPS for aircraft position.'
printf, lgf, info_str
if (keyword_set(verbose)) then print, info_str
lid_aimms_read, lgf=lgf, aimms=aimms, verbose=verbose
hor.alt = aimms.alt
hor.lat = aimms.lat
hor.lon = aimms.lon
;hor.ptc = aimms.ptc ; not validated
;hor.rll = aimms.rll ; not validated
endif
latint = hor.lat
lonint = hor.lon
; flag and fix invalid data
index = where(~finite(hor.alt), cnt, complement=valid, ncomplement=nvalid)
if (nvalid LE 0) then message, 'No valid altitude data!'
if (cnt GT 1) then begin
hor[index].alt = $
interpol(hor[valid].alt, hor[valid].tim, hor[index].tim)
hor[index].invalid = 1
info_str = 'Some altitude data are undefined. Interpolated.'
printf, lgf, info_str
message, info_str, /continue
endif
index = where(~finite(hor.ptc) OR ~finite(hor.rll), cnt)
if (cnt GT 1) then begin
hor[index].ptc = 0.0D
hor[index].rll = 0.0D
hor[index].invalid = 1
info_str = 'Some pitch/roll data are undefined. Set to zero.'
printf, lgf, info_str
message, info_str, /continue
endif
index = where(~finite(hor.lat) OR ~finite(hor.lon) $
OR (hor.lat EQ 0.0D AND hor.lon EQ 0.0D), cnt, $
complement=valid, ncomplement=nvalid)
if (cnt GT 1) then begin
latint[index] = interpol(hor[valid].lat, hor[valid].tim, hor[index].tim)
lonint[index] = interpol(hor[valid].lon, hor[valid].tim, hor[index].tim)
hor[index].lat = !values.d_nan
hor[index].lon = !values.d_nan
hor[index].invalid = 1
info_str = 'Some latitude/longitude data are undefined.'
printf, lgf, info_str
message, info_str, /continue
endif
index = where(hor.alt LT ymin, cnt)
if (cnt GT 0) then begin
info_str = 'Some altitudes are too small!'
printf, lgf, info_str
message, info_str, /continue
endif
index = where(hor.alt GT ymax, cnt)
if (cnt GT 0) then begin
info_str = 'Some altitudes are too large!'
printf, lgf, info_str
message, info_str, /continue
endif
; off-nadir angle
hor.con = cos((hor.ptc-pitch_offset)*!dpi/180.0D) * cos(hor.rll*!dpi/180.0D)
; geoid correction
avg_lat = mean(hor.lat, /nan)
avg_lon = mean(hor.lon, /nan)
;if (keyword_set(geoid_corr)) then begin
if (!version.os_family EQ 'unix') then begin
cmd = './' + intpt_ux
tmpfl = '/tmp/intpt.inp'
endif else begin
cmd = '.\' + intpt_win
tmpfl = tmpdir + 'intpt.inp'
endelse
geoid = 0.0D
openw, dat, tmpfl, /get_lun
printf, dat, format='(%"%12.8f %12.8f")', avg_lat, avg_lon
free_lun, dat
pushd, f77_path
spawn, cmd, geoid
popd
geoid = double(geoid[0])
geoid = read_ascii(tmpdir + 'intpt.out')
hor.alt -= geoid
info_str = string(format='(%"AVGLAT=%7.3f AVGLON=%7.3f ' $
+ 'GEOID=%6.1fm")', avg_lat, avg_lon, geoid)
;endif else begin
; info_str = string(format='(%"AVGLAT=%7.3f AVGLON=%7.3f ' $
; + '*** geoid correction disabled ***")', avg_lat, avg_lon)
;endelse
printf, lgf, info_str
if (keyword_set(verbose)) then print, info_str
; compute along-track distance
inc_distance = sqrt((latint[1:*] - latint)^2 $
+ cos(!dpi * (latint + latint[1:*]) / 360.0D)^2 $
* (lonint[1:*] - lonint)^2) * 60.0D * 1852.0D
idx = where(~finite(inc_distance) $
OR inc_distance GT (hor[1:*].tim - hor.tim) * 3.0D * default_speed, cnt)
if (cnt GT 0) then $
inc_distance[idx] = (hor[idx+1].tim - hor[idx].tim) * default_speed
hor[0].dis = 0.0D
hor[1:*].dis = total(inc_distance, /cumulative)
skip:
n = n_elements(hor)
; orography
if (orogfln NE '') then begin
info_str = string(format='(%"Orography: %s")', orogfln)
printf, lgf, info_str
if (keyword_set(verbose)) then print, info_str
; the restore produces the arrays orog_lat, orog_lon, orog_surf
restore, orog_path + orogfln
hor.osf = interpolate2d(hor.lon, hor.lat, $
orog_lon, orog_lat, orog_surf, outbound=_hundef_)
endif else begin
hor.osf = _hundef_
endelse
; compute some other flight info
hor.rsf = hor.alt - hor.rhg ; surface (radar)
index = where(hor.rhg LE -2000.0D OR hor.alt GE 2000.0D)
hor[index].rsf = _hundef_
hor.ofn = acos(hor.con) * 180.0D / !dpi
free_lun, lgf
end