cal_val_plots.py 30.3 KB
Newer Older
1
#!/usr/bin/env python3.8
2

3
4
from cal_val import  AERI_lib
from cal_val import bbcal
William Roberts's avatar
William Roberts committed
5
import matplotlib
6

7
matplotlib.use('Agg')
8
import faulthandler
9
import matplotlib.pylab as plt
10
import numpy as np
William Roberts's avatar
William Roberts committed
11
import argparse
12
13
import os, logging, math, numbers
from textwrap import wrap
Marco Kurzynski's avatar
Marco Kurzynski committed
14
from datetime import datetime
Marco Kurzynski's avatar
Marco Kurzynski committed
15
from utils import file_finder
16
faulthandler.enable()
William Roberts's avatar
William Roberts committed
17

18
19
20
21
22
23
24
matplotlib_logger = logging.getLogger('matplotlib')
# set WARNING for Matplotlib.
matplotlib_logger.setLevel(logging.WARNING)

logger = logging.getLogger(__name__)


William Roberts's avatar
William Roberts committed
25
def to_int(arg):
26
    try:
William Roberts's avatar
William Roberts committed
27
        return int(arg)
28
    except ValueError:
William Roberts's avatar
William Roberts committed
29
        return arg
30
31


32
def get_invalid_records(c1_file, c2_file, excluded_records, records):
33
34
35
36
    missing_data1 = AERI_lib.getMetaValues(c1_file, 'missingDataFlag')[0]
    missing_data2 = AERI_lib.getMetaValues(c2_file, 'missingDataFlag')[0]
    record = 1
    invalid_records = []
William Roberts's avatar
William Roberts committed
37
    for data1, data2 in zip(missing_data1, missing_data2):
38
        if data1 != 0 or data2 != 0:
39
            if record not in excluded_records and record in records:
40
41
42
                invalid_records.append(record)
        record += 1
    return invalid_records
William Roberts's avatar
William Roberts committed
43
44


45
46
def format_time(time):
    try:
William Roberts's avatar
William Roberts committed
47
        return '{:02}:{:02.0f}'.format(math.floor(time), time * 60 % 60)
48
    except ValueError:
William Roberts's avatar
William Roberts committed
49
50
        return str(time)

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
def partition_ticks(ax, font_size, num_xticks=None, num_yticks=None, discrete=False):
    xmin = np.nanmin(ax.get_xlim())
    xmax = np.nanmax(ax.get_xlim())
    ymin = np.nanmin(ax.get_ylim())
    ymax = np.nanmax(ax.get_ylim())

    if num_xticks is not None:
        if discrete is True:
            xrange = xmax - xmin
            ax.set_xticks(np.arange(xmin, xmax + 1, np.ceil(xrange / num_xticks)).astype(int))
        else:
            ax.set_xticks(np.linspace(xmin, xmax, num_xticks))
        ax.tick_params(axis='x', labelsize=font_size)

    if num_yticks is not None:
        if discrete is True:
            yrange = ymax - ymin
            ax.set_yticks(np.arange(ymin, ymax + 1, np.ceil(yrange / num_yticks)).astype(int))
        else:
            ax.set_yticks(np.linspace(ymin, ymax, num_yticks))
        ax.tick_params(axis='y', labelsize=font_size)


def wrap_text(text, num_chars=11):
    return '\n'.join(wrap(text, num_chars))
77
78


William Roberts's avatar
William Roberts committed
79
80
def get_plot_strings(title_body, name_body, start_record, end_record, c1_file):
    # Note: This relies on the RNC file name format.
81
    dateYYMMDD = c1_file[-12:-6]
William Roberts's avatar
William Roberts committed
82
83
    timeshhmm = [format_time(time) for time in AERI_lib.getMetaValues(c1_file, 'Time')[0]]
    instrumentUnitNumber = to_int(AERI_lib.getMetaValues(c1_file, 'instrumentUnitNumber')[0][start_record - 1])
84
85

    # Title strings:
86
    title_instrument_info = 'AERI-{0} {1}'.format(instrumentUnitNumber, title_body)
87
    title_date = 'Date: 20{0}-{1}-{2}'.format(dateYYMMDD[:2], dateYYMMDD[2:4], dateYYMMDD[4:])
88
    title_records = 'Record: ' + str(start_record)
William Roberts's avatar
William Roberts committed
89
    title_times = 'Time: ' + timeshhmm[start_record - 1] + ' UTC'
90
91

    # file name strings:
92
    name_record_string = 'record_' + str(start_record)
93
94

    # Change strings if not seperating records.
95
96
    if start_record != end_record:
        title_records = 'Records: {0} - {1}'.format(start_record, end_record)
97
        title_times = 'Time: {0} - {1} UTC'.format(timeshhmm[start_record - 1], timeshhmm[end_record - 1])
98
99
100
        name_record_string = ''
    # Seperated by 8 spaces.
    title = '        '.join([title_instrument_info, title_date, title_records, title_times])
101
102
103
    name = "aeri_{0}.png".format(
        '_'.join([str(string) for string in [instrumentUnitNumber, name_body, name_record_string, dateYYMMDD]
                  if string]))
104
105
106
    return title, name


107
def scatter_plot(x, y, color, excluded_records, **kwargs):
108
109
110
111
112
113
114
115
    x = np.array(x)
    y = np.array(y)
    for record in excluded_records:
        if record < max(x):
            index = np.searchsorted(x, record)
            x = np.insert(x, index, record)
            y = np.insert(y, index, np.nan)
    # Make widths larger or smaller based on size of graph.
116
117
    plt.plot(x, y, color + 'o', markersize=2 + 4 / np.sqrt(len(x)), antialiased=True, **kwargs)
    return plt.plot(x, y, color, linewidth=.75 + 2 / np.sqrt(len(x)), antialiased=True, **kwargs)[0]
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


def set_axis(ax, xstart=None, xstop=None, ystart=None, ystop=None):
    x = [line.get_xdata() for line in ax.lines]
    y = [line.get_ydata() for line in ax.lines]

    axis_range = list(plt.axis())
    if xstart == 'fit':
        axis_range[0] = np.nanmin(x)
    elif xstart is not None:
        axis_range[0] = xstart
    if xstop == 'fit':
        axis_range[1] = np.nanmax(x)
    elif xstop is not None:
        axis_range[1] = xstop
    if ystart == 'fit':
        axis_range[2] = np.nanmin(y)
    elif ystart is not None:
        axis_range[2] = ystart
    if ystop == 'fit':
        axis_range[3] = np.nanmax(y)
    elif ystop is not None:
        axis_range[3] = ystop
    plt.axis(axis_range)


def make_therm_temps_plot(c1_file, NBB, IBB, records, output_path, excluded_records):
    logger.debug('Making thermtemps plot')

    if isinstance(records, numbers.Number):
        records = [records]

    plot_title, filename = get_plot_strings('BB Temps & Diffmeans', 'therm_temps', records[0], records[-1], c1_file)

    lookup = ['Record', 'HBBapexTemp', 'HBBbottomTemp', 'HBBtopTemp', 'ABBapexTemp', 'ABBbottomTemp', 'ABBtopTemp',
              'NBBapexTemp', 'NBBbottomTemp', 'NBBtopTemp', 'IBBapexTemp', 'IBBbottomTemp', 'IBBtopTemp']
    metadata = np.array(AERI_lib.getMetaValues(c1_file, lookup))[:, np.array(records) - 1]

    formatter = plt.FormatStrFormatter('%3.2f')
    font_size = 12
    tick_font_size = 10
William Roberts's avatar
William Roberts committed
159

160
161
    fig = plt.figure(figsize=(13, 7), dpi=240)
    legend = {}
162
163
164
165
166
167
168
169
170

    rows = 3
    if len(NBB) > 0:
        rows = rows + 1
    if len(IBB) > 0:
        rows = rows + 1

    ax1 = plt.subplot(rows, 1, 1)
    ax1.yaxis.set_major_formatter(formatter)
171
172
173
    legend['Apex'] = scatter_plot(records, metadata[1], 'r', excluded_records)
    legend['Bottom'] = scatter_plot(records, metadata[2], 'k', excluded_records)
    legend['Top'] = scatter_plot(records, metadata[3], 'b', excluded_records)
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
    plt.ylabel(wrap_text('HBB Temp (K)'), fontsize=font_size)
    set_axis(ax1, xstart='fit', xstop='fit', ystart='fit', ystop='fit')
    partition_ticks(ax1, tick_font_size, num_xticks=20, discrete=True)
    partition_ticks(ax1, tick_font_size, num_yticks=5)
    ax1.set_xlim(records[0], records[-1])
    plt.grid()

    plt.title(plot_title, fontsize=font_size)

    ax2 = plt.subplot(rows, 1, 2, sharex=ax1)
    ax2.yaxis.set_major_formatter(formatter)
    scatter_plot(records, metadata[4], 'r', excluded_records)
    scatter_plot(records, metadata[5], 'k', excluded_records)
    scatter_plot(records, metadata[6], 'b', excluded_records)
    plt.ylabel(wrap_text('ABB Temp (K)'), fontsize=font_size)
    set_axis(ax2, ystart='fit', ystop='fit')
    partition_ticks(ax2, tick_font_size, num_yticks=5)
    plt.grid()

    ax3 = plt.subplot(rows, 1, 3, sharex=ax1)
    ax3.yaxis.set_major_formatter(formatter)
    scatter_plot(records, metadata[7], 'r', excluded_records)
    scatter_plot(records, metadata[8], 'k', excluded_records)
    scatter_plot(records, metadata[9], 'b', excluded_records)
    plt.ylabel(wrap_text('NBB Temp (K)'), fontsize=font_size)
    set_axis(ax3, ystart='fit', ystop='fit')
    partition_ticks(ax3, tick_font_size, num_yticks=5)
    plt.grid()

William Roberts's avatar
William Roberts committed
203
    # Diffmean plots.
204
205
206
207
208
    if len(NBB) > 0:
        NBB = np.array(NBB)
        ax4 = plt.subplot(rows, 1, 4, sharex=ax1)
        ax4.yaxis.set_major_formatter(formatter)
        # Channel 1 diffmean.
209
        legend['Channel 1'] = scatter_plot(NBB[:, 0], NBB[:, -2], 'm', excluded_records)
210
211
212
        plt.hlines([.079, -.079], [records[0], records[0]], [records[-1], records[-1]],
                   'm', linewidth=1, linestyles='dashed')
        # Channel 2 diffmean.
213
        legend['Channel 2'] = scatter_plot(NBB[:, 0], NBB[:, -1], 'g', excluded_records)
214
215
216
217
218
219
220
221
222
223
224
225
        plt.hlines([.083, -.083], [records[0], records[0]], [records[-1], records[-1]],
                   'g', linewidth=1, linestyles='dashed')
        plt.ylabel(wrap_text('NBB Diffmean (K)'), fontsize=font_size)
        set_axis(ax4, ystart=-.2, ystop=.2)
        partition_ticks(ax4, tick_font_size, num_yticks=5)
        plt.grid()

    if len(IBB) > 0:
        IBB = np.array(IBB)
        ax5 = plt.subplot(rows, 1, 5, sharex=ax1)
        ax5.yaxis.set_major_formatter(formatter)
        # Channel 1 diffmean.
226
        legend['Channel 1'] = scatter_plot(IBB[:, 0], IBB[:, -2], 'm', excluded_records)
227
228
229
        plt.hlines([.237, -.237], [records[0], records[0]], [records[-1], records[-1]],
                   'm', linewidth=1, linestyles='dashed')
        # Channel 2 diffmean.
230
        legend['Channel 2'] = scatter_plot(IBB[:, 0], IBB[:, -1], 'g', excluded_records)
231
232
233
234
235
236
237
238
239
        plt.hlines([.359, -.359], [records[0], records[0]], [records[-1], records[-1]],
                   'g', linewidth=1, linestyles='dashed')
        plt.ylabel(wrap_text('IBB Diffmean (K)'), fontsize=font_size)
        set_axis(ax5, ystart=-1, ystop=1)
        partition_ticks(ax5, tick_font_size, num_yticks=5)
        plt.grid()

    plt.xlabel('Record', fontsize=font_size)
    plt.subplots_adjust(hspace=.3)
240
    leg = fig.legend(legend.values(), legend.keys(), fontsize=font_size,
William Roberts's avatar
William Roberts committed
241
                     loc='upper center', ncol=len(legend), mode="expand")
242
243
    for line in leg.legendHandles:
        line.set_linewidth(2)
244
245
246
247

    plt.savefig(os.path.join(output_path, filename), bbox_inches='tight', dpi=240)
    plt.close('all')

248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
    if len(IBB) > 0:
        plot_title, filename = get_plot_strings('IBB Diffmean', 'icebody_diffmean', records[0], records[-1], c1_file)
        fig = plt.figure(figsize=(13, 7), dpi=240)
        legend = {}

        IBB = np.array(IBB)
        ax = plt.subplot(1, 1, 1)
        ax.yaxis.set_major_formatter(formatter)
        # Channel 1 diffmean.
        legend['Channel 1'] = scatter_plot(IBB[:, 0], IBB[:, -2], 'm', excluded_records)
        plt.hlines([.237, -.237], [records[0], records[0]], [records[-1], records[-1]],
                   'm', linewidth=1, linestyles='dashed')
        # Channel 2 diffmean.
        legend['Channel 2'] = scatter_plot(IBB[:, 0], IBB[:, -1], 'g', excluded_records)
        plt.hlines([.359, -.359], [records[0], records[0]], [records[-1], records[-1]],
                   'g', linewidth=1, linestyles='dashed')
William Roberts's avatar
William Roberts committed
264
        plt.ylabel('IBB Diffmean (K)', fontsize=font_size)
265
266
267
268
269
270
271
272
        set_axis(ax, xstart='fit', xstop='fit', ystart=-1, ystop=1)
        partition_ticks(ax, tick_font_size, num_xticks=20, discrete=True)
        partition_ticks(ax, tick_font_size, num_yticks=11)
        plt.grid()

        plt.title(plot_title, fontsize=font_size)
        plt.xlabel('Record', fontsize=font_size)
        leg = fig.legend(legend.values(), legend.keys(), fontsize=font_size,
William Roberts's avatar
William Roberts committed
273
                         loc='upper center', ncol=len(legend), mode="expand")
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
        for line in leg.legendHandles:
            line.set_linewidth(2)

        plt.savefig(os.path.join(output_path, filename), bbox_inches='tight', dpi=240)
        plt.close('all')

    if len(NBB) > 0:
        plot_title, filename = get_plot_strings('NBB Diffmean', '3rdbody_diffmean', records[0], records[-1], c1_file)
        fig = plt.figure(figsize=(13, 7), dpi=240)
        legend = {}

        NBB = np.array(NBB)
        ax = plt.subplot(1, 1, 1)
        ax.yaxis.set_major_formatter(formatter)
        # Channel 1 diffmean.
        legend['Channel 1'] = scatter_plot(NBB[:, 0], NBB[:, -2], 'm', excluded_records)
        plt.hlines([.079, -.079], [records[0], records[0]], [records[-1], records[-1]],
                   'm', linewidth=1, linestyles='dashed')
        # Channel 2 diffmean.
        legend['Channel 2'] = scatter_plot(NBB[:, 0], NBB[:, -1], 'g', excluded_records)
        plt.hlines([.083, -.083], [records[0], records[0]], [records[-1], records[-1]],
                   'g', linewidth=1, linestyles='dashed')
William Roberts's avatar
William Roberts committed
296
        plt.ylabel('NBB Diffmean (K)', fontsize=font_size)
297
298
299
300
301
302
303
304
        set_axis(ax, xstart='fit', xstop='fit', ystart=-.2, ystop=.2)
        partition_ticks(ax, tick_font_size, num_xticks=20, discrete=True)
        partition_ticks(ax, tick_font_size, num_yticks=11)
        plt.grid()

        plt.title(plot_title, fontsize=font_size)
        plt.xlabel('Record', fontsize=font_size)
        leg = fig.legend(legend.values(), legend.keys(), fontsize=font_size,
William Roberts's avatar
William Roberts committed
305
                         loc='upper center', ncol=len(legend), mode="expand")
306
307
308
309
310
311
        for line in leg.legendHandles:
            line.set_linewidth(2)

        plt.savefig(os.path.join(output_path, filename), bbox_inches='tight', dpi=240)
        plt.close('all')

312
313
314
315
316
317
318
319
320

def make_NBB_plot(c1_file, lw, sw, records, output_path):
    logger.debug('Making NBB plot')

    if isinstance(records, numbers.Number):
        records = [records]

    formatter = plt.FormatStrFormatter('%3.2f')
    font_size = 12
321
322
    tick_font_size = 10

323
    lwv = lw[0]
324
325
    obs_lw = lw[3]
    pred_lw = lw[4]
326
327

    swv = sw[0]
328
329
330
331
332
333
    obs_sw = sw[3]
    pred_sw = sw[4]

    for key, val in {'longwave numbers': lwv, 'observed longwave': obs_lw, 'predicted longwave': pred_lw}.items():
        if np.isnan(val).all():
            logger.warning(
334
                'All nans in NBB {0}, records {1}-{2}. Skipping plot generation'.format(key, records[0], records[-1]))
335
336
            return
        elif np.isnan(val).any():
William Roberts's avatar
William Roberts committed
337
            logger.warning('nan values found in NBB {0}, records {1}-{2}, in the following wave numbers: {3}'.format(
338
                key, records[0], records[-1], ', '.join([str(lwv[i]) for i in np.argwhere(np.isnan(val)).flatten()])))
339
340
341
342

    for key, val in {'shortwave numbers': swv, 'observed shortwave': obs_sw, 'predicted shortwave': pred_sw}.items():
        if np.isnan(val).all():
            logger.warning(
343
                'All nans in NBB {0}, records {1}-{2}. Skipping plot generation'.format(key, records[0], records[-1]))
344
345
            return
        elif np.isnan(val).any():
William Roberts's avatar
William Roberts committed
346
            logger.warning('nan values found in NBB {0}, records {1}-{2}, in the following wave numbers: {3}'.format(
347
                key, records[0], records[-1], ', '.join([str(swv[i]) for i in np.argwhere(np.isnan(val)).flatten()])))
William Roberts's avatar
William Roberts committed
348

349
350
351
352
353
354
355
    y_min = np.nanmin(np.concatenate((pred_lw, pred_sw)))
    y_max = np.nanmax(np.concatenate((pred_lw, pred_sw)))
    y_diff = math.ceil(y_max) - math.floor(y_min)
    # Give at least a 20% margin on the top and bottom.
    y_min = math.floor(y_min - y_diff * .2)
    y_max = math.ceil(y_max + y_diff * .2)

356
    # Create 3rd Body plot.
357
358
359
    fig = plt.figure(figsize=(13, 7), dpi=240)
    legend = {}

360
361
    ax1 = plt.subplot(2, 1, 1)
    ax1.yaxis.set_major_formatter(formatter)
362
    # Placeholders for legend.
363
    phr = plt.plot([0, 1], [0, 0], 'r', linewidth=1.5)
364
    ph = plt.plot([0, 1], [0, 0], linewidth=0)
365
    phk = plt.plot([0, 1], [0, 0], 'k', linewidth=1.5)
366
    ph2 = plt.plot([0, 1], [0, 0], linewidth=0)
367

368
    p1 = plt.plot(lwv, obs_lw, 'r', linewidth=.1, antialiased=True)
369
    legend['Mean Observed NBB Brightness Temperature'] = p1[0]
370
    p3 = plt.plot(swv, obs_sw, 'r', linewidth=.1, antialiased=True)
371
372
373
    p2 = plt.plot(lwv, pred_lw, 'k', linewidth=.5, antialiased=True)
    legend['Mean Predicted NBB Brightness Temperature'] = p2[0]
    p4 = plt.plot(swv, pred_sw, 'k', linewidth=.5, antialiased=True)
374

375
    plt.axis([500, 3000, y_min, y_max])
376
    partition_ticks(ax1, tick_font_size, num_xticks=11, num_yticks=11)
William Roberts's avatar
William Roberts committed
377
    plt.ylabel('Temp (K)', fontsize=font_size)
378

William Roberts's avatar
William Roberts committed
379
380
    plt.grid()

381
    plot_title, filename = get_plot_strings('3rd Body Verification', '3rdbody', records[0], records[-1], c1_file)
382
383
    plt.title(plot_title, fontsize=font_size)

384
385
386
    # NBB difference plot.
    ax2 = plt.subplot(2, 1, 2, sharex=ax1)
    ax2.yaxis.set_major_formatter(formatter)
387
    phg = plt.plot([0, 1], [0, 0], 'g', linewidth=1.5)
388
    phx = plt.plot([0, 1], [0, 0], linewidth=0)
389
390
391
392
    bounds = plt.hlines([.079, -.079, .083, -.083], [900, 900, 2100, 2100], [1100, 1100, 2200, 2200], 'g', linewidth=2)
    legend['Predicted 3 Sigma'] = bounds
    p1 = plt.plot(lwv, obs_lw - pred_lw, 'k', linewidth=.15, antialiased=True)
    p2 = plt.plot(swv, obs_sw - pred_sw, 'k', linewidth=.15, antialiased=True)
393
394

    plt.axis([500, 3000, -.2, .2])
395
    partition_ticks(ax2, tick_font_size, num_xticks=11, num_yticks=9)
William Roberts's avatar
William Roberts committed
396
    plt.ylabel('Delta Temp (K)', fontsize=font_size)
397
398
399

    plt.grid()

400
    plt.xlabel('Wavenumber', fontsize=font_size)
401
    leg = fig.legend(legend.values(), legend.keys(), fontsize=font_size,
William Roberts's avatar
William Roberts committed
402
                     loc='upper center', ncol=len(legend), mode="expand")
403
404
    for line in leg.legendHandles:
        line.set_linewidth(2)
405

406
407
    plt.savefig(os.path.join(output_path, filename), bbox_inches='tight', dpi=240)
    plt.close('all')
William Roberts's avatar
William Roberts committed
408

William Roberts's avatar
William Roberts committed
409

410
411
412
413
414
415
416
def make_IBB_plot(c1_file, lw, sw, records, output_path):
    logger.debug('Making IBB plot')

    if isinstance(records, numbers.Number):
        records = [records]

    formatter = plt.FormatStrFormatter('%3.2f')
417
418
419
    font_size = 14
    tick_font_size = 10

420
421
422
423
424
425
426
427
    lwv = lw[0]
    obs_lw = lw[1]
    pred_lw = lw[2]

    swv = sw[0]
    obs_sw = sw[1]
    pred_sw = sw[2]

428
429
430
    for key, val in {'longwave numbers': lwv, 'observed longwave': obs_lw, 'predicted longwave': pred_lw}.items():
        if np.isnan(val).all():
            logger.warning(
431
                'All nans in IBB {0}, records {1}-{2}. Skipping plot generation'.format(key, records[0], records[-1]))
432
433
            return
        elif np.isnan(val).any():
William Roberts's avatar
William Roberts committed
434
            logger.warning('nan values found in IBB {0}, records {1}-{2}, in the following wave numbers: {3}'.format(
435
                key, records[0], records[-1], ', '.join([str(lwv[i]) for i in np.argwhere(np.isnan(val)).flatten()])))
436
437
438
439

    for key, val in {'shortwave numbers': swv, 'observed shortwave': obs_sw, 'predicted shortwave': pred_sw}.items():
        if np.isnan(val).all():
            logger.warning(
440
                'All nans in IBB {0}, records {1}-{2}. Skipping plot generation'.format(key, records[0], records[-1]))
441
442
            return
        elif np.isnan(val).any():
William Roberts's avatar
William Roberts committed
443
            logger.warning('nan values found in IBB {0}, records {1}-{2}, in the following wave numbers: {3}'.format(
444
                key, records[0], records[-1], ', '.join([str(swv[i]) for i in np.argwhere(np.isnan(val)).flatten()])))
William Roberts's avatar
William Roberts committed
445

446
447
448
449
450
451
452
    y_min = np.nanmin(np.concatenate((pred_lw, pred_sw)))
    y_max = np.nanmax(np.concatenate((pred_lw, pred_sw)))
    y_diff = math.ceil(y_max) - math.floor(y_min)
    # Give at least a 20% margin on the top and bottom.
    y_min = math.floor(y_min - y_diff * .2)
    y_max = math.ceil(y_max + y_diff * .2)

453
    # Create Ice Body plot.
454
455
456
457
    fig = plt.figure(figsize=(13, 7), dpi=240)
    legend = {}

    ax1 = plt.subplot(2, 1, 1)
458
    ax1.yaxis.set_major_formatter(formatter)
459
    # Placeholders for legend.
460
    phr = plt.plot([0, 1], [0, 0], 'r', linewidth=1.5)
461
    ph = plt.plot([0, 1], [0, 0], linewidth=0)
462
    phk = plt.plot([0, 1], [0, 0], 'k', linewidth=1.5)
463
    ph2 = plt.plot([0, 1], [0, 0], linewidth=0)
464

465
    p1 = plt.plot(lwv, obs_lw, 'r', linewidth=.1, antialiased=True)
466
    legend['Mean Observed IBB Brightness Temperature'] = p1[0]
467
    p3 = plt.plot(swv, obs_sw, 'r', linewidth=.1, antialiased=True)
468
469
470
    p2 = plt.plot(lwv, pred_lw, 'k', linewidth=.5, antialiased=True)
    legend['Mean Predicted IBB Brightness Temperature'] = p2[0]
    p4 = plt.plot(swv, pred_sw, 'k', linewidth=.5, antialiased=True)
471

472
    plt.axis([500, 3000, y_min, y_max])
473
    partition_ticks(ax1, tick_font_size, num_xticks=11, num_yticks=11)
William Roberts's avatar
William Roberts committed
474
    plt.ylabel('Temp (K)', fontsize=font_size)
475
476

    plt.grid()
William Roberts's avatar
William Roberts committed
477

478
    plot_title, filename = get_plot_strings('Ice Body Verification', 'icebody', records[0], records[-1], c1_file)
479
480
481
    plt.title(plot_title, fontsize=font_size)

    # IBB difference plot
482
    ax2 = plt.subplot(2, 1, 2, sharex=ax1)
483
    ax2.yaxis.set_major_formatter(formatter)
484
    phg = plt.plot([0, 1], [0, 0], 'g', linewidth=1.5)
485
    phx = plt.plot([0, 1], [0, 0], linewidth=0)
486
487
488
489
    bounds = plt.hlines([.237, -.237, .359, -.359], [900, 900, 2100, 2100], [1100, 1100, 2200, 2200], 'g', linewidth=2)
    legend['Predicted 3 Sigma'] = bounds
    p1 = plt.plot(lwv, obs_lw - pred_lw, 'k', linewidth=.15, antialiased=True)
    p2 = plt.plot(swv, obs_sw - pred_sw, 'k', linewidth=.15, antialiased=True)
490

491
    plt.axis([500, 3000, -1, 1])
492
    partition_ticks(ax2, tick_font_size, num_xticks=11, num_yticks=9)
William Roberts's avatar
William Roberts committed
493
    plt.ylabel('Delta Temp (K)', fontsize=font_size)
494
495
496

    plt.grid()

497
    plt.xlabel('Wavenumber', fontsize=font_size)
498
    leg = fig.legend(legend.values(), legend.keys(), fontsize=font_size,
William Roberts's avatar
William Roberts committed
499
                     loc='upper center', ncol=len(legend), mode="expand")
500
501
    for line in leg.legendHandles:
        line.set_linewidth(2)
502
503
504

    plt.savefig(os.path.join(output_path, filename), bbox_inches='tight', dpi=240)
    plt.close('all')
William Roberts's avatar
William Roberts committed
505

506

Marco Kurzynski's avatar
Marco Kurzynski committed
507
508
509
def main(input_dir, output_dir=None, record_range=None, excluded_records=None, separate=None, no_spectral=None):
    c1_file = os.path.join(input_dir, ( os.path.basename(input_dir)[2:8] + "C1.RNC" ))
    c2_file = os.path.join(input_dir, ( os.path.basename(input_dir)[2:8] + "C2.RNC" ))
Marco Kurzynski's avatar
Marco Kurzynski committed
510
    if not os.path.isfile(c1_file) or not os.path.isfile(c2_file):
Marco Kurzynski's avatar
Marco Kurzynski committed
511
        logger.error(f"Invalid input_dir: {input_dir} -- YYMMDDC1.RNC or YYMMDDC2.RNC not in input_dir")
Marco Kurzynski's avatar
Marco Kurzynski committed
512
513
        exit(1)

Marco Kurzynski's avatar
Marco Kurzynski committed
514
515
516
    if output_dir is not None:
        output_path = os.path.realpath(output_dir)
    else:
Marco Kurzynski's avatar
Marco Kurzynski committed
517
        output_path = input_dir
Marco Kurzynski's avatar
Marco Kurzynski committed
518

519
    if not os.path.isdir(output_path):
Marco Kurzynski's avatar
Marco Kurzynski committed
520
        logger.info('Creating output directory: {0}'.format(output_path))
521
        os.mkdir(output_path)
522

Marco Kurzynski's avatar
Marco Kurzynski committed
523
    times = AERI_lib.getMetaValues(str(c1_file), 'Time')[0]
524
    # 83=NBB, 84=IBB.
Marco Kurzynski's avatar
Marco Kurzynski committed
525
    sceneMirrorPositions = [to_int(pos) for pos in AERI_lib.getMetaValues(str(c1_file), 'sceneMirrorPosition')[0]]
526

Marco Kurzynski's avatar
Marco Kurzynski committed
527
    if record_range is None:
528
        records = range(1, len(times) + 1)
Marco Kurzynski's avatar
Marco Kurzynski committed
529
530
    elif len(record_range.split()) == 1 or len(record_range.split()) == 2:
        records = range(int(record_range.split()[0]), int(record_range.split()[-1]) + 1)
531
    else:
Marco Kurzynski's avatar
Marco Kurzynski committed
532
        logger.error('--record-range has incorrect input: {0}'.format(records))
William Roberts's avatar
William Roberts committed
533
        exit(1)
534

535
536
    logger.info('start through end records: {0} - {1}'.format(records[0], records[-1]))
    if any(record < 1 or record > len(times) for record in records):
William Roberts's avatar
William Roberts committed
537
538
539
        logger.error('Records must be in range 1 to {0}'.format(len(times)))
        exit(1)

Marco Kurzynski's avatar
Marco Kurzynski committed
540
    if excluded_records is None:
541
        excluded_records = []
William Roberts's avatar
William Roberts committed
542
    else:
Marco Kurzynski's avatar
Marco Kurzynski committed
543
        excluded_records = [int(record) for record in excluded_records.split()]
William Roberts's avatar
William Roberts committed
544

Marco Kurzynski's avatar
Marco Kurzynski committed
545
    invalid_records = get_invalid_records(str(c1_file), str(c2_file), excluded_records, records)
546
    if len(invalid_records) > 0:
William Roberts's avatar
William Roberts committed
547
        logger.warning('Records flagged for missing data: {0}. Adding them to excluded records'.format(
548
            ', '.join([str(record) for record in invalid_records])))
549
        excluded_records.extend(invalid_records)
William Roberts's avatar
William Roberts committed
550

551
552
    logger.info('excluded records: {0}'.format(', '.join([str(record) for record in excluded_records])))
    bad_excluded_records = []
553
    excluded_redundant = []
554
    for record in list(set(excluded_records)):
555
556
557
        if record < 1 or record > len(times):
            logger.error('Excluded records must be in range 1 to {0}'.format(len(times)))
            exit(1)
558
559
560
        if record < records[0] or record > records[-1]:
            bad_excluded_records.append(record)
        if excluded_records.count(record) > 1:
561
562
563
564
            excluded_redundant.append(record)
    if len(excluded_redundant) > 0:
        logger.warning('Excluded records that appear more than once: {0}'.format(
            ', '.join([str(record) for record in excluded_redundant])))
565
    if len(bad_excluded_records) > 0:
566
        logger.warning('Excluded records that are outside of record range {0} to {1}: {2}'.format(
567
            records[0], records[-1], ', '.join([str(record) for record in bad_excluded_records])))
568

569
570
    records = [record for record in records if record not in excluded_records]
    if len(records) < 1:
William Roberts's avatar
William Roberts committed
571
572
        logger.error('No records to use. Exiting')
        exit(1)
573

Marco Kurzynski's avatar
Marco Kurzynski committed
574
    if no_spectral is False:
575
576
577
578
579
580
        logger.debug('Skipping spectral graph generation')

    # sceneMirrorPosition = 83
    NBB = []
    # sceneMirrorPosition = 84
    IBB = []
William Roberts's avatar
William Roberts committed
581

Marco Kurzynski's avatar
Marco Kurzynski committed
582
    if separate is True:
Marco Kurzynski's avatar
Marco Kurzynski committed
583
584
585
586
        for record in records:
            logger.info('Record {0} scene mirror position = {1}'.format(record, sceneMirrorPositions[record - 1]))
            # Returns wavenumbers, obs_ice, pred_ice, obs_3body, pred_3body
            logger.debug('Running bbcal')
Marco Kurzynski's avatar
Marco Kurzynski committed
587
            lw,sw = bbcal.process(input_dir, record, output_path=output_path)
Marco Kurzynski's avatar
Marco Kurzynski committed
588
589
590
591
592
593
594
595
596

            if sceneMirrorPositions[record - 1] == 83:
                NBB.append((record, times[record - 1], lw[5], sw[5]))
                if no_spectral is False:
                    make_NBB_plot(str(c1_file), lw, sw, record, output_path)
            elif sceneMirrorPositions[record - 1] == 84:
                IBB.append((record, times[record - 1], lw[6], sw[6]))
                if no_spectral is False:
                    make_IBB_plot(str(c1_file), lw, sw, record, output_path)
597
598
599
600
601
602
603
604
605
606
607
608

        if len(NBB) > 0:
            with open(os.path.join(output_path, 'NBB.avg'), 'w') as file:
                file.write("Record,Time,NBB Chan1 DIFFMean,NBB Chan2 DIFFMean\n")
                for row in NBB:
                    file.write(','.join([str(col) for col in row]) + '\n')

        if len(IBB) > 0:
            with open(os.path.join(output_path, 'IBB.avg'), 'w') as file:
                file.write("Record,Time,IBB Chan1 DIFFMean,IBB Chan2 DIFFMean\n")
                for row in IBB:
                    file.write(','.join([str(col) for col in row]) + '\n')
Marco Kurzynski's avatar
Marco Kurzynski committed
609
    elif no_spectral is False:
610
611
612
613
614
        logger.info('Record {0} scene mirror position = {1}'.format(records[0],
                                                                    sceneMirrorPositions[records[0] - 1]))
        if len(records) > 1:
            logger.info('Record {0} scene mirror position = {1}'.format(records[-1],
                                                                        sceneMirrorPositions[records[-1] - 1]))
615
        # Returns wavenumbers, obs_ice, pred_ice, obs_3body, pred_3body
Marco Kurzynski's avatar
Marco Kurzynski committed
616
        lw,sw = bbcal.process(input_dir, records, output_path=output_path)
617

Marco Kurzynski's avatar
Marco Kurzynski committed
618
619
        make_NBB_plot(str(c1_file), lw, sw, records, output_path)
        make_IBB_plot(str(c1_file), lw, sw, records, output_path)
William Roberts's avatar
William Roberts committed
620

Marco Kurzynski's avatar
Marco Kurzynski committed
621
    make_therm_temps_plot(str(c1_file), NBB, IBB, records, output_path, excluded_records)
622
623

if __name__ == '__main__':
Marco Kurzynski's avatar
Marco Kurzynski committed
624
625
626
627
628
    #
    # Values to change with different data sets in box below
    #
    ######################################################################
    # Argument Parsing
Marco Kurzynski's avatar
Marco Kurzynski committed
629
630
631
    parser = argparse.ArgumentParser(
            description='Generate cal val plots for IBB and NBB. The first record is record 1\n\n' \
            'environment variables:\n' \
632
633
634
635
636
637
638
639
640
641
            f"{'  IDNBB' : <24}NBB Scene ID (defaults to 'S')\n" \
            f"{'  IDIBB' : <24}IBB Scene ID (defaults to 'T')\n" \
            f"{'  NBEMFL' : <24}NBB Emissivity File (defaults to 'nbb-emis.asc')\n" \
            f"{'  IBEMFL' : <24}IBB Emissivity File (defaults to 'ibb-emis.asc')\n" \
            f"{'  CAVFAC2' : <24}BB Cavity Factor (defaults to '39.0')\n" \
            f"{'  WTNB' : <24}NBB Apex, Bottom and Top weightings (defaults to '0.811 0.189 0.0')\n" \
            f"{'  TOSNB' : <24}NBB Apex, Bottom and Top offsets (defaults to '0.0 0.0 0.0')\n" \
            f"{'  WTIB' : <24}IBB Apex, Bottom and Top weightings (defaults to '0.811 0.0945 0.0945')\n" \
            f"{'  TOSIB' : <24}IBB Apex, Bottom and Top offsets (defaults to '0.0 0.0 0.0')\n" \
            f"{'  NRNGS' : <24}Read number of ranges (defaults to '1')\n\n" \
642
643
            f"{'  IBBTEMP' : <24}IBB temperature (defaults to '317.925')\n\n" \
            f"{'  NBBTEMP' : <24}NBB temperature (no default)\n\n" \
Marco Kurzynski's avatar
Marco Kurzynski committed
644
645
646
647
648
649
650
651
            f"{'  C1_MLTFILE' : <24}MLT output filename for channel 1 (defaults to 'bbcal_YYMMDD.mlt1')\n" \
            f"{'  C1_LOGFILE' : <24}Log output filename for channel 1 (defaults to 'bbcal_YYMMDD.lo1')\n" \
            f"{'  C1_SUMFILE' : <24}Sum output filename for channel 1 (defaults to 'bbcal_YYMMDD.sum1')\n" \
            f"{'  C1_WN' : <24}Lower and upper wavenumber for channel 1 (defaults to '900 1100')\n\n" \
            f"{'  C2_MLTFILE' : <24}MLT output filename for channel 2 (defaults to 'bbcal_YYMMDD.mlt2')\n" \
            f"{'  C2_LOGFILE' : <24}Log output filename for channel 2 (defaults to 'bbcal_YYMMDD.lo2')\n" \
            f"{'  C2_SUMFILE' : <24}Sum output filename for channel 2 (defaults to 'bbcal_YYMMDD.sum2')\n" \
            f"{'  C2_WN' : <24}Lower and upper wavenumber for channel 2 (defaults to '2100 2200')",
652
            formatter_class=argparse.RawTextHelpFormatter
653
            )
Marco Kurzynski's avatar
Marco Kurzynski committed
654

Marco Kurzynski's avatar
Marco Kurzynski committed
655
    parser.add_argument("input_dir", help="Directory containing the channel files in HHMMSSCX.RNC format. Realtive or absolute path")
Marco Kurzynski's avatar
Marco Kurzynski committed
656
    parser.add_argument("-o", '--output-dir',
Marco Kurzynski's avatar
Marco Kurzynski committed
657
                        help="Output file directory. It will be created if it does not exist. Defaults to current directory. Relative or absolute path")
Marco Kurzynski's avatar
Marco Kurzynski committed
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
    parser.add_argument("-r", '--record-range',
                        help="Record range. Either a number or a string of two numbers separated by a space. Defaults to all records")
    parser.add_argument('-e', '--excluded-records',
                        help='Records to exclude. Either a number or a string of numbers separated by spaces')
    parser.add_argument('-s', '--separate', action='store_true',
                        help='Plot each record individually. Also generates IBB.avg and NBB.avg files')
    parser.add_argument('-n', '--no-spectral', action='store_true',
                        help='Skip generation of spectral plots and only plot thermistor temperatures and diffmeans.')
    parser.add_argument('-v', '--verbose', action="count", default=0,
                        help='Each occurrence increases verbosity 1 level through ERROR-WARNING-INFO-DEBUG.')

    args = parser.parse_args()

    # Logging setup.
    logger_levels = [logging.ERROR, logging.WARNING, logging.INFO, logging.DEBUG]
    numpy_levels = ['ignore', 'ignore', 'ignore', 'warn']
    logging_format = '[%(levelname)s: %(asctime)s : %(name)s] %(message)s'
    logging.basicConfig(format=logging_format, datefmt='%Y-%m-%d %H:%M:%S', level=logger_levels[min(3, args.verbose)])
    np.seterr(all=numpy_levels[min(3, args.verbose)])

    # File format: HHMMSSC.RNC
Marco Kurzynski's avatar
Marco Kurzynski committed
679
    ae_dirs = file_finder.getDirs(args.input_dir)
Marco Kurzynski's avatar
Marco Kurzynski committed
680
681

    if not ae_dirs:
Marco Kurzynski's avatar
Marco Kurzynski committed
682
        logger.error(f"Invalid input_dir: {args.input_dir} -- No AEYYMMDD "
Marco Kurzynski's avatar
Marco Kurzynski committed
683
                     f"directories found, must pass a folder contaning AEYYMMDD dirs, "
684
                     f"or an AEYYMMDD dir")
Marco Kurzynski's avatar
Marco Kurzynski committed
685
686
687
        exit(1)

    for d in ae_dirs:
Marco Kurzynski's avatar
Marco Kurzynski committed
688
        main(d, output_dir=args.output_dir, record_range=args.record_range,
Marco Kurzynski's avatar
Marco Kurzynski committed
689
690
                excluded_records=args.excluded_records,
                separate=args.separate, no_spectral=args.no_spectral)