Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / pycama / cama_plot_time.py @ 833:354ddca04053

History | View | Annotate | Download (58.9 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
        ylabel = []
162
        yunits = set()
163
        for var_name in data['variables']:
164
            variable = data[var_name]
165
            time_ax = time_transform(variable['time'])
166
            if parameters is not None:
167
                actions = parameters
168
            else:
169
                actions = ['mean', 'median', 'mode', 'inter_quartile_range', 'standard_deviation', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99']
170
                if not disable_min_max:
171
                    actions = actions.append('max')
172
                    actions = actions.append('min')
173
            for k in actions:
174
                pltlabel="{0} ({1})".format(variable['meta']['title'], k)
175
                plt.plot_date(time_ax,
176
                              variable[k], label=pltlabel,
177
                              linestyle=('-' if not k.startswith('q') else ':'),
178
                              xdate=True, ydate=False, fmt='-',
179
                              linewidth=(1.5 if k in ('mean', 'median', 'inter_quartile_range', 'standard_deviation') else 0.8))
180
            ylabel.append(variable['meta']['title'])
181
            yunits.add(variable['meta']['units'])
182
        plt.xlabel('Date')
183
        if 'count' in actions:
184
            plt.ylabel("Number of retrievals per day")
185
        else:
186
            plt.ylabel("{0} [{1}]".format(", ".join(ylabel), ", ".join(list(yunits))))
187
        plt.legend(loc='best', ncol=2, fontsize='xx-small', frameon=False)
188
        ax = fig.gca()
189
        date_format = mdates.DateFormatter('%Y-%m-%d')
190
        date_loc = mdates.AutoDateLocator(minticks=6, maxticks=14, interval_multiples=True)
191
        ax.xaxis.set_major_formatter(date_format)
192
        ax.xaxis.set_major_locator(date_loc)
193
        fig.autofmt_xdate()
194

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

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

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

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

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

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

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

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

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

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

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

    
447
    rval['orbits'] = np.asarray(x_list)
448

    
449
    return rval
450

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

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

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

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

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

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

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

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

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

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

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

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

    
563
    return rval
564

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

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

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

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

    
607
    rval = {}
608
    rval['times'] = []
609

    
610
    files_n_granules = {}
611
    files = set()
612
    n_steps = 0
613

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

    
635
    input_files = sorted(list(files_n_granules.keys()))
636

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

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

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

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

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

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

    
676
    return rval
677

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

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

    
689
    t_end.append(stop[-1])
690

    
691
    return t_start, t_end, rval
692

    
693
def outline_table(logger, data, **kwargs):
694
    pass
695

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
1034

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
1143
def cmdline():
1144
    descr = """Plot time-dependent data from PyCAMA."""
1145
    parser = argparse.ArgumentParser(description=descr,
1146
                                     prog=os.path.basename(sys.argv[0]))
1147
    parser.add_argument('--version', action='version',
1148
                        version='%(prog)s {0}'.format(__version__))
1149
    parser.add_argument('-l', '--log', metavar='LEVEL', choices=('debug', 'info', 'warning', 'error'),
1150
                        help="Set the log level ('debug', 'info', 'warning', 'error')", default="info")
1151
    parser.add_argument('-o', '--output', choices=('screen', 'png', 'svg', 'pdf'),
1152
                        help="Determine the destination ('screen', 'png', 'svg', 'pdf')",
1153
                        default="screen")
1154
    parser.add_argument('-f', '--file', dest='outfile', help="Output file", default=None)
1155
    parser.add_argument('--time-start', dest='time_start', type=timeParser, default=None,
1156
                        help="Start of the time-range to plot. From start of provided data if not given.")
1157
    parser.add_argument('--time-end', dest='time_end', type=timeParser, default=None,
1158
                        help="End of the time-range to plot. Up to end of provided data if not given.")
1159

    
1160
    subparsers = parser.add_subparsers(help='Extraction class')
1161

    
1162
    parser_l = subparsers.add_parser('list', help='List variables')
1163
    parser_l.add_argument('input_files', metavar='FILE', nargs="+",
1164
                          help='The netCDF files extracted by PyCAMA')
1165
    parser_l.set_defaults(func=list_variables)
1166

    
1167
    parser_a = subparsers.add_parser('alongtrack', help='Along track data', aliases=['along'])
1168
    parser_a.add_argument('-V', '--variable', required=True,
1169
                          help='Name of the variable to extract.')
1170
    parser_a.add_argument('-R', '--set-range', type=float, nargs=2, default=None, help="Set data range for plot.")
1171
    parser_a.add_argument('-P', '--parameter', help='Parameter to extract for all rows.',
1172
                          choices=('iqr', 'min', 'max', 'mean', 'median', 'mode', 'stddev', 'count', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99'))
1173
    parser_a.add_argument('-r', '--row', type=int, help="Row to display as a time-line with all parameters.")
1174
    parser_a.add_argument('input_files', metavar='FILE', nargs="+",
1175
                          help='The netCDF files extracted by PyCAMA')
1176
    parser_a.set_defaults(func=alongtrack)
1177

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

    
1193
    parser_h = subparsers.add_parser('histogram', help='Histogram data',
1194
                                     aliases=['hist', 'histo', 'histomap'])
1195
    parser_h.add_argument('-V', '--variable', required=True, nargs="+",
1196
                          help='Name of the variable(s) to extract. Multiple values are only allowed for a parameter plot.')
1197
    parser_h.add_argument('-H', '--histomap', action='store_true', dest='show_hist',
1198
                          help="time-series of histograms")
1199
    parser_h.add_argument('-P', '--parameter', action='store_const', dest='show_hist',
1200
                          default=False, const=False,
1201
                          help="time-series of variable (min, max, mean, median, ...)")
1202
    parser_h.add_argument('-e', '--equalized', action='store_true',
1203
                          help='When showing a histomap, use the equalized histogram')
1204
    parser_h.add_argument('-M', '--disable_min_max', action='store_true',
1205
                          help="Do not include the min/max values in the plot.")
1206
    parser_h.add_argument('-p', '--parameters', nargs="+", choices=('iqr', 'min', 'max', 'mean', 'median', 'mode', 'stddev', 'count', 'q01', 'q05', 'q10', 'q25', 'q75', 'q90', 'q95', 'q99'),
1207
                          help="Select which parameter(s) to plot. ")
1208
    parser_h.add_argument('input_files', metavar='FILE', nargs="+",
1209
                          help='The netCDF files extracted by PyCAMA')
1210
    parser_h.set_defaults(func=histomap)
1211

    
1212
    parser_i = subparsers.add_parser('irradiance', help='Irradiance data', aliases=['irrad', 'irr'])
1213
    parser_i.add_argument('-V', '--variable', help='The name of the variable to extract.')
1214
    parser_i.add_argument('-2', '--E2-only', action='store_true', help="Include only values from E2")
1215
    parser_i.add_argument('input_files', metavar='FILE', nargs="+",
1216
                          help='The netCDF files extracted by PyCAMA')
1217
    parser_i.set_defaults(func=irradiance)
1218

    
1219
    parser_o = subparsers.add_parser('outline', help='Outline data')
1220
    parser_o.add_argument('--processing-time', action='store_true', dest="processing_time", help="Make a plot of the processing time per granule")
1221
    parser_o.add_argument('--orbit-range', nargs=2, type=int, dest='orbit_range', help='Orbit range to display when showing processing time')
1222
    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")
1223
    parser_o.add_argument('input_files', metavar='FILE', nargs="+",
1224
                          help='The netCDF files extracted by PyCAMA')
1225
    parser_o.set_defaults(func=outline)
1226

    
1227
    parser_s = subparsers.add_parser('scatter', help='Scatter data (not supported)')
1228
    parser_s.add_argument('input_files', metavar='FILE', nargs="+",
1229
                          help='The netCDF files extracted by PyCAMA')
1230
    parser_s.set_defaults(func=None)
1231

    
1232
    parser_w = subparsers.add_parser('world', help='Level 3 data (not supported)',
1233
                                     aliases=['l3', 'level3'])
1234
    parser_w.add_argument('-V', '--variable', help='The name of the variable to extract.')
1235
    parser_w.add_argument('-m', '--movie', default=None,
1236
                          help='Where to write the movie (a single slice will be witten as normal)')
1237
    parser_w.add_argument('-p', '--projection', choices=('north_pole', 'south_pole', 'world', 'magic'),
1238
                          help='What to plot')
1239
    parser_w.add_argument('input_files', metavar='FILE', nargs="+",
1240
                          help='The netCDF files extracted by PyCAMA')
1241
    parser_w.set_defaults(func=None)
1242

    
1243
    parser_z = subparsers.add_parser('zonal', help='Zonal mean',
1244
                                     aliases=['zone'])
1245
    parser_z.add_argument('-V', '--variable', help='The name of the variable to extract.')
1246
    z_group = parser_z.add_mutually_exclusive_group(required=False)
1247
    z_group.add_argument('-a', '--autorange', action='store_true', help="Auto-scale the output")
1248
    z_group.add_argument('-R', '--range', dest='fixedrange', nargs=2, type=float, help="Scale the zonal average to this range")
1249
    parser_z.add_argument('input_files', metavar='FILE', nargs="+",
1250
                          help='The netCDF files extracted by PyCAMA')
1251
    parser_z.set_defaults(func=zonal)
1252

    
1253

    
1254
    args = parser.parse_args()
1255

    
1256
    return args, parser
1257

    
1258

    
1259

    
1260
if __name__ == "__main__":
1261
    main()