Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / pycama / cama_plot_time.py @ 850:8e96afa2665a

History | View | Annotate | Download (59.4 KB)

1
#!/usr/bin/env python3
2

    
3
import argparse
4
import sys
5
import re
6
import os.path
7
import datetime
8
import traceback
9
import logging
10
import getpass
11
import pprint
12
import math
13
import itertools
14
from collections import OrderedDict
15
import warnings
16
warnings.filterwarnings("ignore", category=FutureWarning)
17

    
18
import numpy as np
19
import h5py, netCDF4
20
import matplotlib as mpl
21
mpl.use('Agg') # plot without X
22
from matplotlib import pyplot as plt
23
import matplotlib.cm as cm
24
from matplotlib import ticker
25
import matplotlib.dates as mdates
26
from matplotlib import colors as mcolors
27
from matplotlib import gridspec
28
from matplotlib.colors import Normalize, LogNorm
29

    
30
from scipy import polyfit
31

    
32
import pycama.transform as transform
33
from pycama import Variable, Reader, WorldPlot, ZonalPlot, HistogramPlot, ScatterPlot, AlongTrackPlot, OutlinePlot, EventPlot, IrradiancePlot, Report, JobOrder
34
from pycama.utilities import setup_logging, extract_time, extract_input_count, extract_time_step_count, get_variable_list, count_time_steps, CAMAException, JobOrderError
35
from pycama.utilities import total_errors, total_warnings, compare_versions, set_warnings_limit, check_warning_counter, read_version_number_from_file, extract_granule_count
36

    
37
from pycama import __version__
38

    
39
def main():
40
    args, parser = cmdline()
41

    
42
    logger = setup_logging(production_logger=True, loglevel=args.log, errlevel="error")
43
    start_time = datetime.datetime.utcnow()
44

    
45
    logger.info("Starting PyCAMA (cama_plot_time.py) version {0}".format(__version__))
46
    logger.info("PyCAMA copyright \u00A9 2016-2020, KNMI (Maarten Sneep). See LICENSE file for details.")
47
    logger.info("Running under Python version {0[0]}.{0[1]}.{0[2]} ({0[3]}), numpy {1}, netCDF4 {2} ({4}), h5py {3} ({5})".format(sys.version_info,
48
                np.__version__, netCDF4.__version__, h5py.version.version, netCDF4.getlibversion().split()[0], h5py.version.hdf5_version))
49

    
50
    try:
51
        args.input_files = sorted(args.input_files)
52
    except AttributeError:
53
        logger.error("No input files supplied.")
54
        parser.print_help(file=sys.stderr)
55
        sys.exit(1)
56

    
57
    if args.func is None:
58
        logger.error("No valid command given.")
59
        parser.print_help(file=sys.stderr)
60
        sys.exit(1)
61

    
62
    try:
63
        data = args.func(logger, **vars(args))
64
    except (CAMAException, NotImplementedError) as err:
65
        logger.error(str(err))
66
        sys.exit(1)
67

    
68
    if data is None:
69
        sys.exit(0)
70

    
71
    fig = plt.figure(figsize=(25/2.54, 16/2.54))
72
    try:
73
        plot_func = globals()[args.func.__name__ + '_plot']
74
    except KeyError:
75
        logger.error("Plotting function '{0}' not found.".format(args.func.__name__ + '_plot'))
76
        sys.exit(1)
77

    
78
    try:
79
        plot_func(logger, fig, data, **vars(args))
80
    except (CAMAException, NotImplementedError) as err:
81
        logger.error(str(err))
82
        sys.exit(1)
83

    
84
    if args.output == "screen":
85
        fig.show()
86
        s = input('Hit enter\n\n')
87
    else:
88
        logger.info("Saving figure to '{0}'".format(args.outfile))
89
        fig.savefig(args.outfile, format=args.output)
90
    if total_errors() > 0:
91
        rval = 255
92
    elif total_warnings() > 0:
93
        rval = 127
94
    else:
95
        rval = 0
96

    
97
    elapsed_time = datetime.datetime.utcnow() - start_time
98
    logger.info("Processing time: {0:.1f} seconds".format(elapsed_time.total_seconds()))
99
    logger.info("Errors: {0}, warnings: {1}".format(total_errors(), total_warnings()))
100
    logger.info("Stopping with exit code: %d", rval)
101
    logging.shutdown()
102
    sys.exit(rval)
103

    
104

    
105
def time_transform(times):
106
    """Convert time to days sine 0001-01-01 for matplotlib plotting.
107

108
    Datetime objects are rounded down to midnight before the time in times.
109
    """
110
    return mpl.dates.date2num(times)
111

    
112
def select_time_range(input_files, time_start, time_end, logger):
113
    n_steps = 0
114
    times = []
115
    files = {}
116
    for src in input_files:
117
        steps_in_file = extract_time_step_count(src)
118
        for idx in range(steps_in_file):
119
            logger.debug("processing file '{0}'".format(os.path.basename(src)))
120
            t, tstart, tend = extract_time(src, idx)
121
            if (time_start is None or t >= time_start) and (time_end is None or t <= time_end):
122
                times.append(t)
123
                n_steps += 1
124
                if src in files:
125
                    files[src].append(idx)
126
                else:
127
                    files[src] = [idx]
128
    return files, times, n_steps
129

    
130

    
131
def histomap_plot(logger, fig, data, variable=None, show_hist=True, equalized=False, disable_min_max=False, parameters=None, **kwargs):
132
    if show_hist:
133
        data = data[data['variables'][0]]
134
        x_lims = mdates.date2num([data['time'][0], data['time'][-1]])
135
        x_lims[-1]+=1
136
        y_lims = data['meta']['data_range']
137
        ax = fig.gca()
138
        # Using ax.imshow we set two keyword arguments. The first is extent.
139
        # We give extent the values from x_lims and y_lims above.
140
        # We also set the aspect to "auto" which should set the plot up nicely.
141
        im_data = data['histomap'][:, ::-1]
142
        im = ax.imshow(im_data.T,
143
                       extent=[x_lims[0], x_lims[1],  y_lims[0], y_lims[1]],
144
                       aspect='auto', interpolation='nearest', cmap=cm.nipy_spectral)
145

    
146
        plt.xlabel('Date')
147
        plt.ylabel("{0} [{1}]".format(data['meta']['title'], data['meta']['units']))
148
        cbar = plt.colorbar(im, ax=ax, orientation='vertical')
149
        cbar.set_label("Density" if equalized else "Count")
150

    
151
        # We tell Matplotlib that the x-axis is filled with datetime data,
152
        # this converts it from a float (which is the output of date2num)
153
        # into a nice datetime string.
154
        ax.xaxis_date()
155
        date_format = mdates.DateFormatter('%Y-%m-%d')
156
        date_loc = mdates.AutoDateLocator(minticks=6, maxticks=14, interval_multiples=True)
157
        ax.xaxis.set_major_formatter(date_format)
158
        ax.xaxis.set_major_locator(date_loc)
159
        fig.autofmt_xdate()
160
    else:
161
        translation = {'stddev':'standard_deviation',
162
                       'iqr':'inter_quartile_range'}
163
        ylabel = []
164
        yunits = set()
165
        for var_name in data['variables']:
166
            variable = data[var_name]
167
            time_ax = time_transform(variable['time'])
168
            if parameters is not None:
169
                actions = [translation.get(p, p) for p in parameters]
170
            else:
171
                actions = ['mean', 'median', 'mode', 'inter_quartile_range', 'standard_deviation', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99']
172
                if not disable_min_max:
173
                    actions = actions.append('max')
174
                    actions = actions.append('min')
175
            for k in actions:
176
                pltlabel="{0} ({1})".format(variable['meta']['title'], k.replace('_', ' '))
177
                plt.plot_date(time_ax,
178
                              variable[k], label=pltlabel,
179
                              linestyle=('-' if not k.startswith('q') else ':'),
180
                              xdate=True, ydate=False, fmt='-',
181
                              linewidth=(1.5 if k in ('mean', 'median', 'inter_quartile_range', 'standard_deviation') else 0.8))
182
            ylabel.append(variable['meta']['title'])
183
            yunits.add(variable['meta']['units'])
184
        plt.xlabel('Date')
185
        if 'count' in actions:
186
            plt.ylabel("Number of retrievals per day")
187
        else:
188
            plt.ylabel("{0} [{1}]".format(", ".join(ylabel), ", ".join(list(yunits))))
189
        plt.legend(loc='best', ncol=2, fontsize='xx-small', frameon=False)
190
        ax = fig.gca()
191
        date_format = mdates.DateFormatter('%Y-%m-%d')
192
        date_loc = mdates.AutoDateLocator(minticks=6, maxticks=14, interval_multiples=True)
193
        ax.xaxis.set_major_formatter(date_format)
194
        ax.xaxis.set_major_locator(date_loc)
195
        fig.autofmt_xdate()
196

    
197
def histomap(logger, input_files=None, variable=None, show_hist=True, equalized=False, disable_min_max=False, **kwargs):
198
    logger.debug("Extracting histogram data for plotting")
199
    if input_files is None:
200
        raise CAMAException("Missing input files")
201
    if variable is None:
202
        raise CAMAException("Missing variable to extract.")
203
    
204
    if show_hist:
205
        variables = [variable[0]]
206
    else:
207
        variables = variable
208
    
209
    files_n_indices, times, n_steps = select_time_range(input_files, kwargs['time_start'], kwargs['time_end'], logger)
210
    files = sorted(files_n_indices.keys())
211
    
212
    logger.info("Extracting metadata for variable '{0}'".format(", ".join(variables)))
213
    d = HistogramPlot(filename=files[0], time_index=files_n_indices[files[0]][0])
214
    if d is None or not d.success:
215
        raise CAMAException("Histogram data not recorded")
216
    
217
    rvals = {'variables': variables}
218
    
219
    for variable in variables:
220
        rval = {'time': times}
221
        rval['meta'] = d.variables_meta[variable]
222
        v_name = variable + ('.histogram_eq' if equalized else '.histogram')
223
        n_hist = d.variables[v_name].shape[0]
224
        rval['v_name'] = v_name
225
        rval['n_hist'] = n_hist
226
        rval['hist_axis'] = d.variables[variable+'.bins']
227
        for k in rval['meta'].keys():
228
            if k not in ('noscanline', 'title', 'units', 'color_scale', 'log_range', 'data_range'):
229
                rval[k] = np.zeros((n_steps,), dtype=np.float32)
230

    
231
        logger.info("Extracting {2}histogram of variable '{0}' for {1} time steps".format(variable, n_steps, "equalized " if equalized else ""))
232
        rval['histomap'] = np.zeros((n_steps, n_hist), dtype=(np.float32 if equalized else np.int32))
233
        step=0
234
        rvals[variable] = rval
235
        
236
    for src in files:
237
        logger.info("Processing file '{0}'".format(os.path.basename(src)))
238
        steps_in_file = files_n_indices[src]
239
        for i, idx in enumerate(steps_in_file):
240
            d = HistogramPlot(filename=src, time_index=idx)
241
            for variable in variables:
242
                v_name = rvals[variable]['v_name']
243
                n_hist = rvals[variable]['n_hist']
244
                try:
245
                    v = d.variables[v_name][:]
246
                except KeyError:
247
                    continue
248
                rvals[variable]['histomap'][step+i, :] = v
249
                for k in rvals[variable]['meta'].keys():
250
                    if k not in ('noscanline', 'title', 'units', 'color_scale', 'log_range', 'data_range', 'v_name', 'n_hist'):
251
                        try:
252
                            rvals[variable][k][step+i] = d.variables_meta[variable][k]
253
                        except KeyError:
254
                            pass
255
        step += len(steps_in_file)
256
    return rvals
257

    
258
def events_plot(logger, fig, data, error=False, filters=False,
259
                general=False, warning=False, group=False, fraction=False, nolog=False, **kwargs):
260
    marker = itertools.cycle(('.', '+', 'o', '*', '^', '<', '>', 's', 'x', 'd', '1', '2', '3', 'v'))
261
    color = itertools.cycle(('k', 'r', 'b', 'g', 'm', 'c', 'grey', 'darkviolet'))
262
    items = [key for key in data.keys() if key not in ('orbits', 'number_of_groundpixels')]
263
    if len(items) == 0:
264
        return False
265
    for k in items:
266
        if fraction:
267
            if k not in ('orbits', 'number_of_groundpixels') and np.any(data[k] > 0):
268
                label = k.replace('_', ' ').replace('number of ', '')
269
                if nolog:
270
                    plt.plot(data['orbits'], data[k]/data['number_of_groundpixels'], '-',
271
                                 marker=next(marker), color=next(color),
272
                                 label=label)
273
                else:
274
                    plt.semilogy(data['orbits'], data[k]/data['number_of_groundpixels'], '-',
275
                                 marker=next(marker), color=next(color),
276
                                 label=label)
277
        else:
278
            if k not in ('orbits',) and np.any(data[k] > 0):
279
                label = k.replace('_', ' ').replace('number of ', '')
280
                if nolog:
281
                    plt.plot(data['orbits'], data[k], '-',
282
                             marker=next(marker), color=next(color),
283
                             label=label)
284
                else:
285
                    plt.semilogy(data['orbits'], data[k], '-',
286
                                 marker=next(marker), color=next(color),
287
                                 label=label)
288

    
289
    plt.xlabel('Orbit')
290
    plt.ylabel('Fraction of total pixels' if fraction else "Count")
291
    x_formatter = ticker.ScalarFormatter(useOffset=False)
292
    #y_formatter = ticker.ScalarFormatter(useOffset=False)
293
    ax = fig.gca()
294
    ax.xaxis.set_major_formatter(x_formatter)
295
    #ax.yaxis.set_major_formatter(y_formatter)
296
    plt.legend(loc='best', ncol=3, fontsize='xx-small', frameon=False)
297
    plt.grid(True)
298

    
299
def events(logger, input_files=None, error=False, filters=False,
300
           general=False, warning=False, group=False, specific=None, nolog=False, **kwargs):
301
    logger.info("Extracting event data for plotting")
302
    if input_files is None:
303
        raise CAMAException("Missing input files")
304

    
305
    counters = {}
306
    files = set()
307
    for src in input_files:
308
        with h5py.File(src, 'r') as ref:
309
            try:
310
                product = list(ref['eventplot_data'].keys())[0]
311
            except KeyError:
312
                raise CAMAException("Events not recorded")
313

    
314
            for granule in ref['eventplot_data'][product].keys():
315
                logger.debug("Found granule '{0}' in file {1}.".format(granule, os.path.basename(src)))
316
                t = timeParser(granule.split('_')[2])
317
                if not (('time_start' not in kwargs or kwargs['time_start'] is None or t >= kwargs['time_start'])
318
                        and ('time_end' not in kwargs or kwargs['time_end'] is None or t <= kwargs['time_end'])):
319
                    continue
320
                orbit = int(granule.split('_')[1])
321
                files.add(src)
322
                if orbit in counters:
323
                    counters[orbit] += 1
324
                else:
325
                    counters[orbit] = 1
326
    if group:
327
        granule_count = len(counters.keys())
328
    else:
329
        granule_count = sum([n for n in counters.values()])
330

    
331
    input_files = sorted(list(files))
332
    logger.debug("Granule count: {0} in {1} files".format(granule_count, len(input_files)))
333

    
334
    rval = {}
335
    with h5py.File(input_files[0], 'r') as ref:
336
        product = list(ref['eventplot_data'].keys())[0]
337
        if 'time_start' not in kwargs or kwargs['time_start'] is None:
338
            granule = sorted(list(ref['eventplot_data'][product].keys()))[0]
339
        else:
340
            granule = sorted([s for s in list(ref['eventplot_data'][product].keys()) if timeParser(s.split('_')[2]) >= kwargs['time_start']])[0]
341
        start_orbit = min(counters.keys())
342
        keys = [k for k in list(ref['eventplot_data'][product][granule].attrs.keys()) if k != 'comment']
343

    
344
    
345
    err_list = ('number_of_generic_exception_occurrences',
346
            'number_of_input_spectrum_missing_occurrences',
347
            'number_of_radiance_missing_occurrences',
348
            'number_of_rejected_pixels_not_enough_spectrum',
349
            'number_of_irradiance_missing_occurrences',
350
            'number_of_groundpixels')
351
    error_keys = [k for k in keys if 'error' in k or k in err_list]
352
    filter_keys = [k for k in keys if 'filter' in k and k != 'number_of_cloud_filter_convergence_error_occurrences']
353
    filter_keys.append('number_of_groundpixels')
354
    warn_list = ('number_of_sun_glint_correction_occurrences',
355
                 'number_of_groundpixels',
356
                 'number_of_pixel_level_input_data_missing_occurrences')
357
    warning_keys = [k for k in keys if 'warning' in k or k in warn_list]
358
    general_keys = ['number_of_groundpixels',
359
            'number_of_processed_pixels',
360
            'number_of_successfully_processed_pixels',
361
            'number_of_failed_retrievals',
362
            'number_of_errors',
363
            'number_of_warnings',
364
            'number_of_filtered_pixels',
365
            'number_of_missing_scanlines']
366
            
367
    if specific is not None:
368
        specific_keys = [k for k in keys if specific in k]
369
        if not specific_keys:
370
            raise CAMAException("Specific key '{0}' not found in {1}".format(specific, ", ".join(keys)))
371
        
372
        specific_keys.extend(['number_of_groundpixels'])
373
        
374
    if specific is not None and len(specific_keys) > 1:
375
        logger.info("Extracting a specific key only '{0}'".format(specific_keys[0]))
376
        for k in specific_keys:
377
            rval[k] = np.zeros((granule_count,), dtype=np.int32)
378
    elif error:
379
        logger.info("Extracting error keys only")
380
        for k in error_keys:
381
            rval[k] = np.zeros((granule_count,), dtype=np.int32)
382
    elif warning:
383
        logger.info("Extracting warning keys only")
384
        for k in warning_keys:
385
            rval[k] = np.zeros((granule_count,), dtype=np.int32)
386
    elif filters:
387
        logger.info("Extracting filter keys only")
388
        for k in filter_keys:
389
            rval[k] = np.zeros((granule_count,), dtype=np.int32)
390
    elif general:
391
        logger.info("Extracting general keys only")
392
        for k in general_keys:
393
            rval[k] = np.zeros((granule_count,), dtype=np.int32)
394
    else:
395
        logger.info("Extracting all keys")
396
        for k in keys:
397
            rval[k] = np.zeros((granule_count,), dtype=np.int32)
398
            
399
    granule_idx = 0
400
    dorbit = 0
401
    x_list = [float(start_orbit)] if group else []
402
    for src in input_files:
403
        logger.info("Processing file '{0}'".format(os.path.basename(src)))
404
        with h5py.File(src, 'r') as ref:
405
            product = list(ref['eventplot_data'].keys())[0]
406

    
407
            if 'time_start' not in kwargs or kwargs['time_start'] is None:
408
                if 'time_end' not in kwargs or kwargs['time_end'] is None:
409
                    sorted_granules = sorted(list(ref['eventplot_data'][product].keys()))
410
                else:
411
                    sorted_granules = sorted([s for s in list(ref['eventplot_data'][product].keys()) if timeParser(s.split('_')[2]) <= kwargs['time_end']])
412
            else:
413
                if 'time_end' not in kwargs or kwargs['time_end'] is None:
414
                    sorted_granules = sorted([s for s in list(ref['eventplot_data'][product].keys()) if timeParser(s.split('_')[2]) >= kwargs['time_start']])
415
                else:
416
                    sorted_granules = sorted([s for s in list(ref['eventplot_data'][product].keys()) if kwargs['time_start'] <= timeParser(s.split('_')[2]) <= kwargs['time_end']])
417
            
418
            if len(sorted_granules) == 0:
419
                continue
420
            
421
            for granule in sorted_granules:
422
                orbit = int(granule.split('_')[1])
423
                grp = ref['eventplot_data'][product][granule]
424
                for att_name in grp.attrs.keys():
425
                    if att_name in rval:
426
                        try:
427
                            rval[att_name][granule_idx] += grp.attrs[att_name]
428
                        except TypeError:
429
                            continue
430
                    if general and att_name not in ('number_of_groundpixels', ):
431
                        if att_name in error_keys:
432
                            rval['number_of_errors'][granule_idx] += grp.attrs[att_name]
433
                        elif att_name in filter_keys:
434
                            rval['number_of_filtered_pixels'][granule_idx] += grp.attrs[att_name]
435
                        elif att_name in warning_keys:
436
                            rval['number_of_warnings'][granule_idx] += grp.attrs[att_name]
437
                if not group:
438
                    granule_idx += 1
439
                    if orbit != start_orbit:
440
                        dorbit = 0
441
                        start_orbit = orbit
442
                    x_list.append(float(orbit) + dorbit/counters[orbit])
443
                    dorbit += 1
444
                if group and orbit != start_orbit:
445
                    granule_idx += 1
446
                    start_orbit = orbit
447
                    x_list.append(float(orbit))
448

    
449
    rval['orbits'] = np.asarray(x_list)
450

    
451
    return rval
452

    
453
def alongtrack_plot(logger, fig, data, variable=None, parameter=None, row=None, **kwargs):
454
    if parameter is not None:
455
        x_lims = mdates.date2num([data['time'][0], data['time'][-1]])
456
        x_lims[-1]+=1
457
        y_lims = [0, data[parameter].shape[1]-1]
458
        ax = fig.gca()
459
        if 'set_range' in kwargs and kwargs['set_range'] is not None:
460
            data_range = sorted(kwargs['set_range'])
461
        else:
462
            # get rid of fill values
463
            data[parameter] = np.where(data[parameter] > 1e36, np.nan, data[parameter])
464
            data_range = [np.nanmin(data[parameter]), np.nanmax(data[parameter])]
465
        logger.info("Data range for {1}({2}): [{0[0]}, {0[1]}]".format(data_range, variable, parameter))
466
        normalizer = Normalize(vmin=data_range[0], vmax=data_range[1], clip=True)
467

    
468
        # Using ax.imshow we set two keyword arguments. The first is extent.
469
        # We give extent the values from x_lims and y_lims above.
470
        # We also set the aspect to "auto" which should set the plot up nicely.
471
        im = ax.imshow(data[parameter].T, norm=normalizer,
472
                       extent=[x_lims[0], x_lims[1],  y_lims[0], y_lims[1]],
473
                       aspect='auto', interpolation='nearest', cmap=cm.nipy_spectral)
474

    
475
        plt.xlabel('Date')
476
        plt.ylabel('Row')
477
        cbar = plt.colorbar(im, ax=ax, orientation='vertical')
478
        cbar.set_label("{2} of {0} [{1}]".format(data['meta']['title'], data['meta']['units'], parameter))
479

    
480
        # We tell Matplotlib that the x-axis is filled with datetime data,
481
        # this converts it from a float (which is the output of date2num)
482
        # into a nice datetime string.
483
        ax.xaxis_date()
484
        date_format = mdates.DateFormatter('%Y-%m-%d')
485
        date_loc = mdates.AutoDateLocator(minticks=6, maxticks=8, interval_multiples=True)
486
        ax.xaxis.set_major_formatter(date_format)
487
        ax.xaxis.set_major_locator(date_loc)
488
        fig.autofmt_xdate()
489
    elif row is not None:
490
        time_ax = time_transform(data['time'])
491
        for k in ('mean', 'median', 'iqr', 'stddev', 'min', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99', 'max'):
492
            plt.plot_date(time_ax,
493
                          data[k], label=k, linestyle=('-' if not k.startswith('q') else ':'),
494
                          xdate=True, ydate=False, fmt='-',
495
                          linewidth=(2.0 if k in ('mean', 'median', 'iqr', 'stddev') else 1.0))
496
        plt.xlabel('Date')
497
        plt.ylabel("{0} [{1}], row {2}".format(data['meta']['title'], data['meta']['units'], row))
498
        plt.legend(loc='best', ncol=4, fontsize='x-small', frameon=False)
499
        ax = fig.gca()
500
        date_format = mdates.DateFormatter('%Y-%m-%d')
501
        date_loc = mdates.AutoDateLocator(minticks=6, maxticks=14, interval_multiples=True)
502
        ax.xaxis.set_major_formatter(date_format)
503
        ax.xaxis.set_major_locator(date_loc)
504
        fig.autofmt_xdate()
505

    
506
def alongtrack(logger, input_files=None, variable=None, parameter=None, row=None, **kwargs):
507
    logger.info("Extracting alongtrack data for plotting")
508
    translation = {'stddev':'standard_deviation',
509
                   'iqr':'inter_quartile_range'}
510

    
511
    if input_files is None:
512
        raise CAMAException("Missing input files")
513
    if variable is None:
514
        raise CAMAException("Missing variable to extract.")
515

    
516
    files_n_indices, times, n_steps = select_time_range(input_files, kwargs['time_start'], kwargs['time_end'], logger)
517
    files = sorted(files_n_indices.keys())
518

    
519
    with h5py.File(files[0], 'r') as ref:
520
        ground_pixel = ref['ground_pixel'].shape[0]
521

    
522
    rval = {'time': times}
523
    logger.info("Extracting metadata for variable '{0}'".format(variable))
524
    d = AlongTrackPlot(filename=files[0], time_index=0)
525
    if d is None or not d.success:
526
        raise CAMAException("Along track data not recorded")
527
    rval['meta'] = d.variables_meta[variable]
528

    
529
    if parameter is not None:
530
        logger.info("Extracting parameter '{0}' for variable '{1}' in {2} time steps for {3} ground pixels".format(parameter, variable, n_steps, ground_pixel))
531
        rval[parameter] = np.zeros((n_steps, ground_pixel), dtype=np.float32)
532

    
533
        step=0
534
        for src in files:
535
            logger.info("Processing file '{0}'".format(os.path.basename(src)))
536
            steps_in_file = files_n_indices[src]
537
            for i, idx in enumerate(steps_in_file):
538
                d = AlongTrackPlot(filename=src, time_index=idx)
539
                try:
540
                    v = d.variables[variable][translation.get(parameter, parameter)]
541
                except KeyError:
542
                    raise CAMAException("Parameter '{0}' for variable '{1}' not found".format(parameter, variable))
543
                rval[parameter][step+i, :] = v
544
            step += len(steps_in_file)
545
    elif row is not None:
546
        logger.info("Extracting variable '{0}' for row {1} in {2} time steps".format(variable, row, n_steps))
547
        if row >= ground_pixel or row < 0:
548
            raise CAMAException("Row must be larger than 0 and smaller than the ground_pixel dimension ({0}). Found {1}.".format(ground_pixel, row))
549

    
550
        for k in ('iqr', 'min', 'max', 'mean', 'median', 'mode', 'stddev', 'count', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99'):
551
            rval[k] = np.zeros((n_steps,), dtype=np.float32)
552
        step=0
553
        for src in files:
554
            steps_in_file = files_n_indices[src]
555
            for i,idx in enumerate(steps_in_file):
556
                d = AlongTrackPlot(filename=src, time_index=idx)
557
                for k in rval.keys():
558
                    if k in ('time', 'meta'):
559
                        continue
560
                    rval[k][step+i] = d.variables[variable][translation.get(k, k)][row]
561
            step += steps_in_file
562
    else:
563
        raise CAMAException("Either row or parameter must be specified.")
564

    
565
    return rval
566

    
567
def irradiance_plot(logger, fig, data, variable=None, **kwargs):
568
    #pp = pprint.PrettyPrinter(indent=2, width=120)
569
    #pp.pprint(data)
570
    x_lims = mdates.date2num([data['times'][0], data['times'][-1]])
571
    y_lims = [0, data[variable].shape[-1]-1]
572
    ax = fig.gca()
573
    # Using ax.imshow we set two keyword arguments. The first is extent.
574
    # We give extent the values from x_lims and y_lims above.
575
    # We also set the aspect to "auto" which should set the plot up nicely.
576
    drange = data[variable].min(), data[variable].max()
577
    norm = mcolors.Normalize(vmin=drange[0], vmax=drange[1], clip=True)
578
    im = ax.imshow(data[variable].T,
579
                   extent=[x_lims[0], x_lims[1],  y_lims[0], y_lims[1]],
580
                   aspect='auto', interpolation='nearest',
581
                   cmap=getattr(cm, data['color_scale'], "seismic"),
582
                   norm=norm)
583

    
584
    plt.xlabel('Date')
585
    plt.ylabel("Ground pixel index")
586
    cbar = plt.colorbar(im, ax=ax, orientation='vertical')
587
    if data['units'] != "1":
588
        cbar.set_label("{0} [{1}]".format(data['title'], data['units']))
589
    else:
590
        cbar.set_label(data['title'])
591

    
592
    # We tell Matplotlib that the x-axis is filled with datetime data,
593
    # this converts it from a float (which is the output of date2num)
594
    # into a nice datetime string.
595
    ax.xaxis_date()
596
    date_format = mdates.DateFormatter('%Y-%m-%d')
597
    date_loc = mdates.AutoDateLocator(minticks=6, maxticks=14, interval_multiples=True)
598
    ax.xaxis.set_major_formatter(date_format)
599
    ax.xaxis.set_major_locator(date_loc)
600
    fig.autofmt_xdate()
601

    
602
def irradiance(logger, input_files=None, variable=None, **kwargs):
603
    logger.debug("Extracting irradiance data for plotting")
604
    if input_files is None:
605
        raise CAMAException("Missing input files")
606
    if variable is None:
607
        raise CAMAException("Missing variable to extract.")
608

    
609
    rval = {}
610
    rval['times'] = []
611

    
612
    files_n_granules = {}
613
    files = set()
614
    n_steps = 0
615

    
616
    for src in input_files:
617
        logger.debug("processing file '{0}'".format(os.path.basename(src)))
618
        with h5py.File(src, 'r') as ref:
619
            try:
620
                irr_granules = sorted(list(ref['irradianceplot_data'].keys()))
621
            except KeyError:
622
                raise CAMAException("Events not recorded")
623
            for i, irr_granule in enumerate(irr_granules):
624
                t = timeParser(irr_granule[11:26])
625
                if not ((kwargs['time_start'] is None or t >= kwargs['time_start'])
626
                        and (kwargs['time_end'] is None or t <= kwargs['time_end'])):
627
                    continue
628
                rval['times'].append(t)
629
                files.add(src)
630
                if src in files_n_granules:
631
                    files_n_granules[src].append((i,irr_granule))
632
                else:
633
                    files_n_granules[src] = [(i,irr_granule)]
634
            if src in files_n_granules:
635
                n_steps += len(files_n_granules[src])
636

    
637
    input_files = sorted(list(files_n_granules.keys()))
638

    
639
    d = IrradiancePlot(filename=input_files[0], time_index=0)
640
    if d is None or not d.success:
641
        raise CAMAException("Irradiance data not recorded")
642

    
643
    if variable not in d.variables[list(d.variables.keys())[0]]:
644
        raise CAMAException("Variable '{0}' not in input".format(variable))
645
    else:
646
        n_groundpixel = d.variables[list(d.variables.keys())[0]][variable]['data'].shape[0]
647

    
648
    for k in ('units', 'title', 'data_range', 'color_scale', 'log_range'):
649
        rval[k] = d.variables_meta[variable][k]
650

    
651
    logger.info("Extracting irradiance data for variable '{0}' in {1} steps".format(variable, n_steps))
652

    
653
    rval[variable] = np.zeros((n_steps, n_groundpixel), dtype=np.float32)
654
    rval['orbits'] = np.zeros((n_steps,), dtype=np.int32)
655
    step=0
656
    for src in input_files:
657
        logger.info("Processing file '{0}'".format(os.path.basename(src)))
658
        steps_in_file = [s[0] for s in files_n_granules[src]]
659
        granules = sorted([s[1] for s in files_n_granules[src]])
660
        for i, idx in enumerate(steps_in_file):
661
            granule = granules[i]
662
            d = IrradiancePlot(filename=src, time_index=idx)
663
            if granule not in d.variables:
664
                continue
665
            v = d.variables[granule][variable]['data']
666
            rval[variable][step+i, :] = v
667
            rval['orbits'][step+i] = d.variables[granule][variable]['orbit']
668
        step+=len(steps_in_file)
669

    
670
    if 'E2_only' in kwargs and kwargs['E2_only']:
671
        idx = rval['orbits'] >= 2818
672
        logger.debug("Filtering for E2 only: total {0} orbits, {1} >= 2818. Found {2} times.".format(n_steps, np.sum(idx), len(rval['times'])))
673
        rval['orbits'] = rval['orbits'][idx]
674
        rval[variable] = rval[variable][idx, :]
675
        if np.any(np.logical_not(idx)):
676
            rval['times'][0] = datetime.datetime.strptime("2018-04-30T00:19:00", "%Y-%m-%dT%H:%M:%S")
677

    
678
    return rval
679

    
680
def stable_ranges(start, stop, d):
681
    t_start = [start[0]]
682
    t_end = []
683
    rval = [d[0]]
684

    
685
    for i in range(1, len(d)):
686
        if d[i] != d[i-1]:
687
            t_end.append(stop[i-1])
688
            t_start.append(start[i])
689
            rval.append(d[i])
690

    
691
    t_end.append(stop[-1])
692

    
693
    return t_start, t_end, rval
694

    
695
def outline_table(logger, data, **kwargs):
696
    pass
697

    
698
def outline_plot(logger, fig, data, **kwargs):
699
    # set height ratios for sublots
700
    n = len(data['items'])
701
    logger.debug("Plotting {0} data traces".format(n))
702

    
703
    if 'processing_time' in kwargs and kwargs['processing_time']:
704
        gs = gridspec.GridSpec(1, 1, hspace=0.0, bottom=0.08, left=0.10, right=0.97, top=0.97)
705
        axes = [plt.subplot(gs[0])]
706
        max_proc_time = np.max(data['time_for_processing'])
707
        if 'max_proc_time_s' in kwargs and kwargs['max_proc_time_s']:
708
            max_proc_time = kwargs['max_proc_time_s']
709
        if max_proc_time < 300:
710
            proc_time = data['time_for_processing']
711
            unit = "s"
712
        elif 300 <= max_proc_time < 7200:
713
            proc_time = data['time_for_processing']/60
714
            max_proc_time /= 60
715
            unit = "min"
716
        else:
717
            proc_time = data['time_for_processing']/3600
718
            max_proc_time /= 3600
719
            unit = "hours"
720

    
721
        if 'orbit_range' in kwargs and kwargs['orbit_range']:
722
            axes[0].set_xlim([kwargs['orbit_range'][0], kwargs['orbit_range'][-1]])
723
        else:
724
            axes[0].set_xlim([data['orbit'][0], data['orbit'][-1]])
725

    
726
        if 'max_proc_time_s' in kwargs and kwargs['max_proc_time_s']:
727
            axes[0].set_ylim([0, max_proc_time])
728

    
729
        axes[0].plot(data['orbit'], proc_time, 'k.')
730
        axes[0].set_xlabel("Orbit")
731
        axes[0].set_ylabel("Processing time [{0}]".format(unit))
732
        return
733

    
734
    # reorder to put overal processing status on top.
735
    try:
736
        idx = data['items'].index('processing_status')
737
        del data['items'][idx]
738
        data['items'].insert(0, 'processing_status')
739
    except ValueError:
740
        pass
741

    
742
    logger.debug("Items: {0}".format(", ".join(data['items'])))
743

    
744
    gs = gridspec.GridSpec(n, 1, hspace=0.0, bottom=0.08, left=0.20, right=0.97, top=0.95)
745

    
746
    # build axes array, one for each 'item'.
747
    # axes 1...n share the x-axis with the first instance.
748
    axes = [plt.subplot(gs[0])]
749
    for i in range(1,n):
750
        axes.append(plt.subplot(gs[i], sharex=axes[0]))
751

    
752
    # time-scale (middle of coverage)
753
    t = [t0 + (t1-t0)/2 for t0, t1 in zip(data['time_coverage_start'], data['time_coverage_end'])]
754
    axes[0].set_xlim([data['time_coverage_start'][0], data['time_coverage_end'][-1]])
755

    
756
    # loop over all items for plotting.
757
    for i, ax in enumerate(axes):
758
        name = data['items'][i]
759
        logger.debug("Plotting '{0}'".format(name))
760
        # plot data
761
        if name.startswith('Status_') or name == 'processing_status':
762
            for j, (t_start, t_end) in enumerate(zip(data['time_coverage_start'],
763
                                                     data['time_coverage_end'])):
764
                mplt_start, mplt_end = mdates.date2num([t_start, t_end])
765
                # red (bad)/blue (good) -> colorblind should be able to distinguish this.
766
                # TBA: NO2 status flags.
767
                try:
768
                    low_val = data[name][j].lower()
769
                except AttributeError:
770
                    color = 'red'
771
                    low_val = 'unavailable'
772
                if 'nominal' in low_val or 'nrti' in low_val:
773
                    color = 'royalblue'
774
                elif 'unneeded' in low_val:
775
                    color='lightgrey'
776
                elif 'earth' in low_val or 'backup' in low_val:
777
                    color = "darkorange"
778
                elif 'solar' in low_val or 'external' in low_val:
779
                    color = "green"
780
                else:
781
                    color = 'red'
782
                ax.axvspan(mplt_start, mplt_end, edgecolor=color, facecolor=color)
783
        elif name in ('time_standard_deviation_per_pixel', 'time_for_processing',
784
                    'time_per_pixel', 'time_for_algorithm_initialization'):
785
            ax.plot(mdates.date2num(t), data[name], 'b.')
786
        elif name in ('orbit', ):
787
            t_start, t_end, val = stable_ranges(data['time_coverage_start'],
788
                                                data['time_coverage_end'],
789
                                                data[name])
790
            nval = len(val)
791
            ax.set_ylim(-1,1)
792
            for j, (s, e, v) in enumerate(zip(t_start, t_end, val)):
793
                ax.text(mdates.date2num(s+(e-s)/2), 0, "{0}".format(v), horizontalalignment='center',
794
                      verticalalignment='center', stretch='ultra-condensed', size='xx-small',
795
                      rotation=0 if nval < 20 else 60)
796
        else:
797
            ax.set_ylim(-1,1)
798
            t_start, t_end, val = stable_ranges(data['time_coverage_start'],
799
                                                data['time_coverage_end'],
800
                                                data[name])
801
            nval = len(val)
802
            colors = ['navy', 'darkviolet', 'maroon', 'slategrey', 'darkgreen', 'orangered']
803
            for j, (s, e, v) in enumerate(zip(t_start, t_end, val)):
804
                if not v:
805
                    continue
806
                logger.debug("Start: {0}; End: {1}; val: {2}".format(s, e, v))
807
                mplt_start, mplt_end = mdates.date2num([s, e])
808
                ax.axvspan(mplt_start, mplt_end,
809
                    edgecolor=colors[j % len(colors)],
810
                    facecolor=colors[j % len(colors)])
811
                if len(name) == 10 and name.upper() == name:
812
                    if name == "AUX_NISE__":
813
                        if v == "None":
814
                            continue
815
                        try:
816
                            v = "{0:%Y-%m-%d}".format(datetime.datetime.strptime(v[14:22], "%Y%m%d"))
817
                        except ValueError:
818
                            logger.warning("Datetime parse error in '{0}' ('{1}' does not match '%Y%m%d')".format(v, v[14:22]))
819
                            continue
820
                    elif 'L1B' in name:
821
                        v = v[52:57]
822
                    else:
823
                        vvv = []
824
                        for vv in v.split(','):
825
                            vv = vv.strip()
826
                            if vv == "None":
827
                                vvv.append(vv)
828
                            elif vv[20:35] == "00000000T000000":
829
                                try:
830
                                    dt = datetime.datetime.strptime(vv[52:67], "%Y%m%dT%H%M%S")
831
                                except ValueError:
832
                                    logger.warning("Datetime parse error in '{0}' ('{1} does not match '%Y%m%dT%H%M%S')".format(vv, vv[52:67]))
833
                                    continue
834
                                vvv.append("{0:%Y-%m-%d %H:%M}".format(dt))
835
                            else:
836
                                fmt = "{0:%Y-%m-%d %H:%M}" if 'MET' in name else "{0:%Y-%m-%d}"
837
                                try:
838
                                    dt = datetime.datetime.strptime(vv[20:35], "%Y%m%dT%H%M%S")
839
                                except ValueError:
840
                                    logger.warning("Datetime parse error in '{0}' ('{1}' does not match '{2}')".format(vv, vv[20:35], fmt))
841
                                    continue
842
                                vvv.append(fmt.format(dt))
843
                        v = ", ".join(sorted(vvv))
844
                if nval < 20:
845
                    ax.text(mdates.date2num(s+(e-s)/2), 0, "{0}".format(v),
846
                            horizontalalignment='center',
847
                            verticalalignment='center', color='azure', stretch='condensed', size='x-small')
848
                else:
849
                    ax.text(mdates.date2num(s+(e-s)/2), 0, "{0}".format(v),
850
                            horizontalalignment='center',
851
                            verticalalignment='center', color='azure', stretch='ultra-condensed', size='xx-small',
852
                            rotation=60 if 'L1B' in name else 0)
853

    
854
    # x-axis labeling (in time).
855
    maj_loc = mdates.AutoDateLocator()
856
    maj_format = mdates.DateFormatter('%Y-%m-%d\n%H:%M' if (t[-1] - t[0]) > datetime.timedelta(seconds=26*3600) else '%H:%M')
857
    axes[0].xaxis.set_major_locator(maj_loc)
858
    axes[0].xaxis.set_major_formatter(maj_format)
859

    
860
    translation = {'time_standard_deviation_per_pixel': r'$\sigma$ time per pixel',
861
                   'time_for_processing': 'processing (s)',
862
                   'time_for_algorithm_initialization': 'initialization (s)',
863
                   'revision_control_identifier': 'revision',
864
                   'ProcessingMode': 'processing mode',
865
                   'Status_CTMFCT_CTMANA': 'status TM5 data',
866
                   'Status_reference_spectrum': 'reference spectrum',
867
                   'Status_L2__O3____': 'ozone reference'}
868
    # loop over all items.
869
    for i, ax in enumerate(axes):
870
        name = data['items'][i]
871
        # disable y-ticks and tick marks on y-axis.
872
        plt.setp(ax.get_yticklabels(), visible=False)
873
        ax.yaxis.set_major_locator(plt.NullLocator())
874
        ax.yaxis.set_minor_locator(plt.NullLocator())
875

    
876
        # disable x-tick labels on axes 0...n-2
877
        try:
878
            tt = ax.get_xticklabels()
879
            logger.debug(str(tt))
880
            if i < n-1:
881
                plt.setp(tt, visible=False)
882
        except ValueError:
883
            logger.debug("No valid x-axis for {0}".format(name))
884

    
885
        # label with attribute name
886
        ax.set_ylabel(translation.get(name, name).replace('_', ' ').strip(), rotation=0,
887
                      horizontalalignment='right',
888
                      verticalalignment='center', stretch='condensed', size='small')
889

    
890
    if (t[-1] - t[0]) < datetime.timedelta(seconds=26*3600):
891
        axes[-1].set_xlabel("{0:%Y-%m-%d}".format(data['time_coverage_start'][0]))
892
    axes[0].set_title(data['product'])
893
    axes[-1].tick_params(axis='x', which='major', labelsize='x-small')
894

    
895
def outline(logger, input_files=None, **kwargs):
896
    logger.info("Extracting outline data for plotting the input file names")
897
    if input_files is None:
898
        raise CAMAException("Missing input files")
899
    counters = {}
900
    rval = OrderedDict()
901
    n_granules = 0
902
    keys = set()
903
    status = set()
904
    input_pointer = set()
905
    dates = set()
906
    granule_list = []
907

    
908
    for src in input_files:
909
        n_granules += extract_granule_count(src)
910
        logger.info("Processing file '{0}'.".format(os.path.basename(src)))
911
        with h5py.File(src, 'r') as ref:
912
            try:
913
                product = list(ref['outlineplot_data'].keys())[0]
914
            except KeyError:
915
                raise CAMAException("Outlines not recorded in {0}".format(src))
916

    
917
            for granule in sorted(list(ref['outlineplot_data'][product].keys())):
918
                logger.debug("Found granule '{0}' in file '{1}'.".format(granule, os.path.basename(src)))
919
                rval[granule] = {}
920
                granule_list.append(granule)
921
                ggrp_attrs = dict(ref['outlineplot_data'][product][granule].attrs)
922
                attnames = sorted(list(ggrp_attrs.keys()))
923
                try:
924
                    orbit = int(ggrp_attrs['orbit'])
925
                except ValueError as err:
926
                    logger.warning("The orbit number could not be converted to an integer value")
927
                    orbit = 0
928
                if orbit in counters:
929
                    counters[orbit] += 1
930
                else:
931
                    counters[orbit] = 1
932
                logger.debug("orbit: {0}".format(orbit))
933
                # status for dynamic aux data
934
                for attname in [a for a in attnames if a.startswith('Status_') or a == 'processing_status']:
935
                    v = ggrp_attrs[attname]
936
                    try:
937
                        rval[granule][attname] = v.decode('utf-8', 'ignore')
938
                    except AttributeError:
939
                        rval[granule][attname] = v
940
                    status.add(attname)
941
                    logger.debug("{0}: '{1}'".format(attname, rval[granule][attname]))
942
                # input pointer
943
                for attname in [a for a in attnames if len(a) == 10 and a.upper() == a]:
944
                    if attname.startswith('L1B_RA') or attname.startswith('L2__'):
945
                        continue
946
                    v = ggrp_attrs[attname]
947
                    try:
948
                        rval[granule][attname] = v.decode('utf-8', 'ignore')
949
                    except AttributeError:
950
                        rval[granule][attname] = v
951
                    input_pointer.add(attname)
952
                    logger.debug("{0}: '{1}'".format(attname, rval[granule][attname]))
953
                # generic metadata expressed as float values.
954
                for attname in ('time_standard_deviation_per_pixel', 'time_for_processing',
955
                                'time_per_pixel', 'time_for_algorithm_initialization', ):
956
                    if attname in ggrp_attrs:
957
                        v = ggrp_attrs[attname]
958
                        try:
959
                            rval[granule][attname] = float(v[0])
960
                        except IndexError:
961
                            rval[granule][attname] = float(v)
962
                        keys.add(attname)
963
                        logger.debug("{0}: '{1}'".format(attname, rval[granule][attname]))
964
                # generic metadata expressed as strings
965
                for attname in ('algorithm_version', 'revision_control_identifier',
966
                                'product_version', 'processor_version', 'ProcessingMode'):
967
                    if attname in ggrp_attrs:
968
                        v = ggrp_attrs[attname]
969
                        try:
970
                            rval[granule][attname] = v.decode('utf-8').strip()
971
                        except AttributeError:
972
                            rval[granule][attname] = v
973
                        keys.add(attname)
974
                        logger.debug("{0}: '{1}'".format(attname, rval[granule][attname]))
975
                # date-time objects
976
                for attname in ('time_coverage_start', 'time_coverage_end', 'date_created'):
977
                    if attname in ggrp_attrs:
978
                        v = ggrp_attrs[attname]
979
                        try:
980
                            v = v.decode('utf-8').strip()
981
                        except AttributeError:
982
                            pass
983
                        # use the KNMI/DLR common part of the date/time string.
984
                        v = datetime.datetime.strptime(v[0:19], '%Y-%m-%dT%H:%M:%S')
985
                        rval[granule][attname] = v
986
                        dates.add(attname)
987
                        logger.debug("{0}: '{1:%Y-%m-%dT%H:%M:%S}'".format(attname, rval[granule][attname]))
988
                # optional (only if not fill value)
989
                attname = 'LongitudeOfDaysideNadirEquatorCrossing'
990
                if attname in ggrp_attrs:
991
                    v = ggrp_attrs[attname]
992
                    if math.fabs(float(v[0])) < 400.0:
993
                        rval[granule][attname] = float(v[0])
994
                        keys.add(attname)
995
                        logger.debug("{0}: '{1}'".format(attname, rval[granule][attname]))
996
                # generic metadata expressed as int values
997
                attname = 'orbit'
998
                rval[granule][attname] = int(orbit)
999
                keys.add(attname)
1000

    
1001
    x_list = []
1002
    for orbit in sorted(list(counters.keys())):
1003
        x_list.extend(list(np.linspace(orbit, orbit+1,
1004
                num=counters[orbit], endpoint=False, dtype=np.float32)))
1005

    
1006
    # restructure return value
1007
    sval = OrderedDict()
1008
    sval['fractional_orbits'] = np.asarray(x_list)
1009
    sval['granule_list'] = granule_list
1010
    sval['items'] = []
1011
    sval['product'] = product
1012

    
1013
    for keylist in (status, keys, input_pointer):
1014
        for key in sorted(list(keylist)):
1015
            l = []
1016
            for gr in granule_list:
1017
                try:
1018
                    l.append(rval[gr][key])
1019
                except KeyError:
1020
                    try:
1021
                        l.append("None" if isinstance(l[0], str) else np.nan)
1022
                    except IndexError:
1023
                        l.append("None")
1024
            sval[key] = np.asarray(l)
1025
            sval['items'].append(key)
1026
    for key in sorted(list(dates)):
1027
        l = []
1028
        for gr in granule_list:
1029
            try:
1030
                l.append(rval[gr][key])
1031
            except KeyError:
1032
                l.append(datetime.datetime.utcnow())
1033
        sval[key] = l
1034
    return sval
1035

    
1036

    
1037
def zonal_plot(logger, fig, data, variable=None, **kwargs):
1038
    x_lims = mdates.date2num([data['time'][0], data['time'][-1]])
1039
    x_lims[-1]+=1
1040
    y_lims = [-90.0, 90.0]
1041

    
1042
    if 'autorange' in kwargs and kwargs['autorange']:
1043
        data['range'] = (np.min(data['zonalmap']), np.max(data['zonalmap']))
1044

    
1045
    try:
1046
        if 'fixedrange' in kwargs and len(kwargs['fixedrange']) == 2:
1047
            data['range'] = kwargs['fixedrange']
1048
    except:
1049
        pass
1050

    
1051
    if data['log_range']:
1052
        normalizer = LogNorm(vmin=data['range'][0], vmax=data['range'][1])
1053
    else:
1054
        normalizer = Normalize(vmin=data['range'][0], vmax=data['range'][1])
1055

    
1056
    ax = fig.gca()
1057
    # Using ax.imshow we set two keyword arguments. The first is extent.
1058
    # We give extent the values from x_lims and y_lims above.
1059
    # We also set the aspect to "auto" which should set the plot up nicely.
1060
    im = ax.imshow((data['zonalmap'][:, ::-1]).T, norm=normalizer,
1061
                   extent=[x_lims[0], x_lims[1],  y_lims[0], y_lims[1]],
1062
                   aspect='auto', interpolation='nearest',
1063
                   cmap=getattr(cm, data['colorscale'], 'nipy_spectral'))
1064

    
1065
    plt.xlabel('Date')
1066
    plt.ylabel('Latitude [degrees]')
1067
    cbar = plt.colorbar(im, ax=ax, orientation='vertical')
1068
    cbar.set_label("{0} [{1}]".format(data['title'], data['units']) if data['units'] != "1" else data['title'])
1069

    
1070
    # We tell Matplotlib that the x-axis is filled with datetime data,
1071
    # this converts it from a float (which is the output of date2num)
1072
    # into a nice datetime string.
1073
    ax.xaxis_date()
1074
    date_format = mdates.DateFormatter('%Y-%m-%d')
1075
    date_loc = mdates.AutoDateLocator(minticks=6, maxticks=14, interval_multiples=True)
1076
    ax.xaxis.set_major_formatter(date_format)
1077
    ax.xaxis.set_major_locator(date_loc)
1078
    fig.autofmt_xdate()
1079

    
1080
def zonal(logger, input_files=None, variable=None, **kwargs):
1081
    logger.debug("Extracting zonal mean data for plotting")
1082
    if input_files is None:
1083
        raise CAMAException("Missing input files")
1084
    if variable is None:
1085
        raise CAMAException("Missing variable to extract.")
1086

    
1087
    files_n_indices, times, n_steps = select_time_range(input_files, kwargs['time_start'], kwargs['time_end'], logger)
1088
    files = sorted(files_n_indices.keys())
1089

    
1090
    d = ZonalPlot(filename=files[0], time_index=0)
1091
    if d is None or not d.success:
1092
        d = WorldPlot(filename=input_files[0], time_index=0)
1093
        if d is None or not d.success:
1094
            raise CAMAException("level 3 or zonal data not recorded")
1095

    
1096
    zmean = d.zonal_mean(variable)
1097
    n_zones = zmean['zonal_mean'].shape[0]
1098
    for k in ('title', 'units', 'range', 'log_range', 'colorscale', 'lat_bounds', 'lat'):
1099
        rval[k] = zmean[k]
1100

    
1101
    logger.info("Extracting zonal mean of variable '{0}' for {1} time steps".format(variable, n_steps))
1102
    rval['zonalmap'] = np.ma.zeros((n_steps, n_zones), dtype=(np.float64))
1103
    step=0
1104
    for src in input_files:
1105
        logger.info("Processing file '{0}'".format(os.path.basename(src)))
1106
        steps_in_file = files_n_indices[src]
1107
        for i, idx in enumerate(steps_in_file):
1108
            d = ZonalPlot(filename=src, time_index=idx)
1109
            if d is None or not d.success:
1110
                d = WorldPlot(filename=src, time_index=idx)
1111
                if d is None or not d.success:
1112
                    raise CAMAException("level 3 or zonal data not recorded in {0} for time-step {1}".format(src, idx))
1113
            v = d.zonal_mean(variable)
1114
            rval['zonalmap'][step+i, :] = v['zonal_mean']
1115
            rval['zonalmap'].mask[step+i, :] = v['zonal_mean'].mask
1116
        step += steps_in_file
1117
    if isinstance(rval['range'], bool) and not data_range:
1118
        rval['range'] = (np.min(rval['zonalmap']), np.max(rval['zonalmap']))
1119
    return rval
1120

    
1121
def list_variables_plot(logger, fig, data, **kwargs):
1122
    """Dummy function, shall do nothing."""
1123
    pass
1124

    
1125
def list_variables(logger, input_files=None, **kwargs):
1126
    logger.info("Printing variable names to standard output.")
1127
    with h5py.File(input_files[0], 'r') as ref:
1128
        names = ref['variable_names'][:]
1129
    for name in names:
1130
        print(str(name, encoding='utf-8'))
1131
    return None
1132

    
1133
def timeParser(string):
1134
    formats = ("%Y-%m-%dT%H:%M:%S", "%Y%m%dT%H%M%S", "%Y-%m-%d", "%Y%m%d")
1135
    for format in formats:
1136
        try:
1137
            value = datetime.datetime.strptime(string, format)
1138
        except ValueError:
1139
            continue
1140
        else:
1141
            return value
1142

    
1143
    raise argparse.ArgumentTypeError("The time '{0}' could not be parsed".format(string))
1144

    
1145
def cmdline():
1146
    descr = """Plot time-dependent data from PyCAMA."""
1147
    parser = argparse.ArgumentParser(description=descr,
1148
                                     prog=os.path.basename(sys.argv[0]))
1149
    parser.add_argument('--version', action='version',
1150
                        version='%(prog)s {0}'.format(__version__))
1151
    parser.add_argument('-l', '--log', metavar='LEVEL', choices=('debug', 'info', 'warning', 'error'),
1152
                        help="Set the log level ('debug', 'info', 'warning', 'error')", default="info")
1153
    parser.add_argument('-o', '--output', choices=('screen', 'png', 'svg', 'pdf'),
1154
                        help="Determine the destination ('screen', 'png', 'svg', 'pdf')",
1155
                        default="screen")
1156
    parser.add_argument('-f', '--file', dest='outfile', help="Output file", default=None)
1157
    parser.add_argument('--time-start', dest='time_start', type=timeParser, default=None,
1158
                        help="Start of the time-range to plot. From start of provided data if not given.")
1159
    parser.add_argument('--time-end', dest='time_end', type=timeParser, default=None,
1160
                        help="End of the time-range to plot. Up to end of provided data if not given.")
1161
    
1162
    parameter_choices = ['iqr', 'min', 'max', 'mean', 'median', 'mode', 'stddev', 'count', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99']
1163
    
1164
    histogram_parameter_choices = []
1165
    for parameter_choice in parameter_choices:
1166
        histogram_parameter_choices.append(parameter_choice)
1167
        if parameter_choice != 'count':
1168
            for sub_choice in ('land', 'sea', 'north', 'south'):
1169
                histogram_parameter_choices.append(f"{parameter_choice}_{sub_choice}")
1170
    
1171
    subparsers = parser.add_subparsers(help='Extraction class')
1172

    
1173
    parser_l = subparsers.add_parser('list', help='List variables')
1174
    parser_l.add_argument('input_files', metavar='FILE', nargs="+",
1175
                          help='The netCDF files extracted by PyCAMA')
1176
    parser_l.set_defaults(func=list_variables)
1177

    
1178
    parser_a = subparsers.add_parser('alongtrack', help='Along track data', aliases=['along'])
1179
    parser_a.add_argument('-V', '--variable', required=True,
1180
                          help='Name of the variable to extract.')
1181
    parser_a.add_argument('-R', '--set-range', type=float, nargs=2, default=None, help="Set data range for plot.")
1182
    parser_a.add_argument('-P', '--parameter', help='Parameter to extract for all rows.',
1183
                          choices=parameter_choices)
1184
    parser_a.add_argument('-r', '--row', type=int, help="Row to display as a time-line with all parameters.")
1185
    parser_a.add_argument('input_files', metavar='FILE', nargs="+",
1186
                          help='The netCDF files extracted by PyCAMA')
1187
    parser_a.set_defaults(func=alongtrack)
1188

    
1189
    parser_e = subparsers.add_parser('event', help='Event counting', aliases=['events'])
1190
    parser_e.add_argument('input_files', metavar='FILE', nargs="+",
1191
                          help='The netCDF files extracted by PyCAMA')
1192
    group_e = parser_e.add_mutually_exclusive_group()
1193
    group_e.add_argument("-S", "--specific", help="Track a specific error or warning over time")
1194
    group_e.add_argument('-E', '--error', action="store_true", help="Only display errors")
1195
    group_e.add_argument('-W', '--warning', action="store_true", help="Only display warnings types")
1196
    group_e.add_argument('-F', '--filters', action="store_true", help="Only display filters")
1197
    group_e.add_argument('-G', '--general', action="store_true", help="Display general counts")
1198
    parser_e.add_argument('-g', '--group', action="store_true", help="group NRT granules into orbits")
1199
    parser_e.add_argument('-f', '--fraction', action="store_true", help="Display as a fraction relative to the total number of pixels.")
1200
    parser_e.add_argument('-L', '--log', action="store_false", dest="nolog", default=False, help="Use a logarithmic vertical scale (default)")
1201
    parser_e.add_argument('-l', '--linear', action="store_true", dest="nolog", default=False, help="Use a linear vertical scale")
1202
    parser_e.set_defaults(func=events)
1203

    
1204
    parser_h = subparsers.add_parser('histogram', help='Histogram data',
1205
                                     aliases=['hist', 'histo', 'histomap'])
1206
    parser_h.add_argument('-V', '--variable', required=True, nargs="+",
1207
                          help='Name of the variable(s) to extract. Multiple values are only allowed for a parameter plot.')
1208
    parser_h.add_argument('-H', '--histomap', action='store_true', dest='show_hist',
1209
                          help="time-series of histograms")
1210
    parser_h.add_argument('-P', '--parameter', action='store_const', dest='show_hist',
1211
                          default=False, const=False,
1212
                          help="time-series of variable (min, max, mean, median, ...)")
1213
    parser_h.add_argument('-e', '--equalized', action='store_true',
1214
                          help='When showing a histomap, use the equalized histogram')
1215
    parser_h.add_argument('-M', '--disable_min_max', action='store_true',
1216
                          help="Do not include the min/max values in the plot.")
1217
    parser_h.add_argument('-p', '--parameters', nargs="+", choices=histogram_parameter_choices,
1218
                          help="Select which parameter(s) to plot. ")
1219
    parser_h.add_argument('input_files', metavar='FILE', nargs="+",
1220
                          help='The netCDF files extracted by PyCAMA')
1221
    parser_h.set_defaults(func=histomap)
1222

    
1223
    parser_i = subparsers.add_parser('irradiance', help='Irradiance data', aliases=['irrad', 'irr'])
1224
    parser_i.add_argument('-V', '--variable', help='The name of the variable to extract.')
1225
    parser_i.add_argument('-2', '--E2-only', action='store_true', help="Include only values from E2")
1226
    parser_i.add_argument('input_files', metavar='FILE', nargs="+",
1227
                          help='The netCDF files extracted by PyCAMA')
1228
    parser_i.set_defaults(func=irradiance)
1229

    
1230
    parser_o = subparsers.add_parser('outline', help='Outline data')
1231
    parser_o.add_argument('--processing-time', action='store_true', dest="processing_time", help="Make a plot of the processing time per granule")
1232
    parser_o.add_argument('--orbit-range', nargs=2, type=int, dest='orbit_range', help='Orbit range to display when showing processing time')
1233
    parser_o.add_argument('--max-proc-time', type=float, dest='max_proc_time_s', help="top of the scale in seconds when showing processing time")
1234
    parser_o.add_argument('input_files', metavar='FILE', nargs="+",
1235
                          help='The netCDF files extracted by PyCAMA')
1236
    parser_o.set_defaults(func=outline)
1237

    
1238
    parser_s = subparsers.add_parser('scatter', help='Scatter data (not supported)')
1239
    parser_s.add_argument('input_files', metavar='FILE', nargs="+",
1240
                          help='The netCDF files extracted by PyCAMA')
1241
    parser_s.set_defaults(func=None)
1242

    
1243
    parser_w = subparsers.add_parser('world', help='Level 3 data (not supported)',
1244
                                     aliases=['l3', 'level3'])
1245
    parser_w.add_argument('-V', '--variable', help='The name of the variable to extract.')
1246
    parser_w.add_argument('-m', '--movie', default=None,
1247
                          help='Where to write the movie (a single slice will be witten as normal)')
1248
    parser_w.add_argument('-p', '--projection', choices=('north_pole', 'south_pole', 'world', 'magic'),
1249
                          help='What to plot')
1250
    parser_w.add_argument('input_files', metavar='FILE', nargs="+",
1251
                          help='The netCDF files extracted by PyCAMA')
1252
    parser_w.set_defaults(func=None)
1253

    
1254
    parser_z = subparsers.add_parser('zonal', help='Zonal mean',
1255
                                     aliases=['zone'])
1256
    parser_z.add_argument('-V', '--variable', help='The name of the variable to extract.')
1257
    z_group = parser_z.add_mutually_exclusive_group(required=False)
1258
    z_group.add_argument('-a', '--autorange', action='store_true', help="Auto-scale the output")
1259
    z_group.add_argument('-R', '--range', dest='fixedrange', nargs=2, type=float, help="Scale the zonal average to this range")
1260
    parser_z.add_argument('input_files', metavar='FILE', nargs="+",
1261
                          help='The netCDF files extracted by PyCAMA')
1262
    parser_z.set_defaults(func=zonal)
1263

    
1264

    
1265
    args = parser.parse_args()
1266

    
1267
    return args, parser
1268

    
1269

    
1270

    
1271
if __name__ == "__main__":
1272
    main()