Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / pycama / CAMA_call.py @ 750:ccb55f4d7a99

History | View | Annotate | Download (38.8 KB)

1
#!/usr/bin/env python3
2
# -*- coding: utf-8 -*-
3

    
4
# Copyright 2016-2017 Maarten Sneep, KNMI
5
#
6
# Redistribution and use in source and binary forms, with or without
7
# modification, are permitted provided that the following conditions are met:
8
#
9
# 1. Redistributions of source code must retain the above copyright notice,
10
#    this list of conditions and the following disclaimer.
11
#
12
# 2. Redistributions in binary form must reproduce the above copyright notice,
13
#    this list of conditions and the following disclaimer in the documentation
14
#    and/or other materials provided with the distribution.
15
#
16
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
20
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23
# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26

    
27
## \file CAMA_call.py
28
#  Main entry point for PyCAMA
29
#
30
#  This module parses the command line, and calls the main function.
31

    
32
import datetime
33
import os
34
import sys
35
import logging
36
import plistlib
37
import traceback
38
import math
39
import warnings
40
warnings.filterwarnings("ignore", category=FutureWarning)
41

    
42
import isodate
43
import netCDF4
44
import h5py
45

    
46
import pycama.transform as transform
47
from pycama import Variable, Reader, WorldPlot, ZonalPlot, HistogramPlot, ScatterPlot, AlongTrackPlot, OutlinePlot, EventPlot, IrradiancePlot, Report, JobOrder
48
from pycama.utilities import setup_logging, extract_time, get_variable_list, count_time_steps, CAMAException, JobOrderError
49
from pycama.utilities import total_errors, total_warnings, compare_versions, set_warnings_limit, check_warning_counter, read_version_number_from_file
50

    
51
from pycama import __version__
52

    
53
## \brief The main function for PyCAMA
54
#
55
#  This is the starting point for a run. The function will read the job
56
#  order file and set up the logging system. It also parses the configuration
57
#  file, and creates the variable objects. It will read the input data from the
58
#  L2 files. It will then loop over all actions, and execute the subclass of
59
#  pycama.AnalysisAndPlot.AnalysisAndPlot that belongs to that action.
60
#  In each of these steps data is written to the output netCDF file.
61
#
62
#  If the job order file asks for a report, then the netCDF file is read again,
63
#  and the report is created by the pycama.CAMA_call.cama_plot() function.
64
#
65
#  \param[in] joborder A \c FILE object or path to a job order file.
66
#         The file is treated as read-only.
67
#  \param[in] aggregation An ISO 8601 duration string to indicate the
68
#         aggregation period for the analysis. If this period is shorter
69
#         than the range of the data that is ingested then multiple output
70
#         intervals will be generated. This is an optional parameter.
71
#         Default is \c None, which results in a single output period.
72
#  \param[in] append Append new output to an existing output file (as
73
#         specified in the job order file). This option allows to aggregate
74
#         multiple output intervals in subsequent calls. Equivalent to the
75
#         \c Append_To_Output dynamic processing key in the job order file.
76
#         Optional, default is \c False.
77
#  \param[in] report_file A file path to the location of the output report.
78
#         This is an optional parameter. The default is taken from the
79
#         joborder file. No report is generated if not specified here or in
80
#         the job order file.
81
#  \param[in] report_only This boolean flag changes the meaning of the
82
#         entries in the job order file. Instead of generating the netCDF
83
#         file specified in the output section of the joborder file, this file
84
#         is \em read as input and a report is generated. Thsi construction
85
#         was chosen to allow seamless re-use of the same job order file.
86
#  \param[in] warnings_limit An integer indicating when an error should be
87
#         generated based on the number of warnings. This is to avoid a
88
#         run-away run of PyCAMA. Default is \c None, no limit.
89
#  \author Maarten Sneep
90
#
91
def main(joborder, aggregation=None, append=False,
92
         report_file=None, report_only=False, warnings_limit=None,
93
         actions_argument=None, message_level=None):
94
    rval = 0
95
    logger = None
96
    start_time = datetime.datetime.utcnow()
97
    try:
98
        if h5py.is_hdf5(joborder.name):
99
            raise RuntimeError("The argument must be a joborder (XML) file, not a netCDF4 or HDF5 file.")
100
        jo = JobOrder(joborder, message_level)
101
        logger = jo.logger
102
        if jo.error:
103
            # special error to correctly capture errors while parsing the joborder.
104
            # Error has been passed to the log system already.
105
            raise JobOrderError("Error in JobOrder file")
106
        if not (report_only or jo.report_only):
107
            logger.info("Extracting data from %s in %s mode", jo.product, jo.processing_mode)
108

    
109
            dt = None
110
            if aggregation is not None:
111
                try:
112
                    dt = isodate.parse_duration(aggregation)
113
                    if dt.total_seconds() > 0:
114
                        jo.config['interval_duration'] = dt
115
                except isodate.ISO8601Error:
116
                    logger.warning("Could not parse interval '%s', using value from configuration %s.", aggregation, str(jo.config['interval_duration']))
117
                    dt = None
118

    
119
            if actions_argument is not None:
120
                actions=actions_argument
121
            elif 'actions' in jo.config and jo.processing_mode in jo.config['actions']:
122
                actions = jo.config['actions'][jo.processing_mode]
123
            else:
124
                logger.warning("Configuration error, no actions specified, using default list")
125
                if jo.processing_mode == 'NRTI':
126
                    actions = ['world', 'zonal', 'histogram', 'along_track',            'outline', 'events', 'irradiance']
127
                else:
128
                    actions = ['world', 'zonal', 'histogram', 'along_track', 'scatter', 'outline', 'events', 'irradiance']
129
            logger.info("Performing extraction actions: {0}".format(", ".join(actions)))
130

    
131
            if report_file is None:
132
                report_file = jo.report_file
133

    
134
            if 'spatial_resolution' in jo.config:
135
                spatial_resolution = jo.config['spatial_resolution']
136
            else:
137
                logger.warning("Configuration error, no spatial resolution specified, using default of \u00BC\u00B0")
138
                spatial_resolution = 0.25
139

    
140
            if 'zonal_resolution' in jo.config:
141
                zonal_resolution = jo.config['zonal_resolution']
142
            else:
143
                zonal_resolution = 1.0
144

    
145
            if 'histogram_bincount' in jo.config:
146
                histogram_bin_count = jo.config['histogram_bincount']
147
            else:
148
                logger.warning("Configuration error, no default histogram size specified in configuration, using default of 100")
149
                histogram_bin_count = 100
150

    
151
            if warnings_limit is not None:
152
                set_warnings_limit(warnings_limit)
153
            elif 'warnings_limit' in jo.config:
154
                set_warnings_limit(jo.config['warnings_limit'])
155

    
156
            if 'max_sza' in jo.config:
157
                max_sza = jo.config['max_sza']
158
            else:
159
                max_sza = 88.0
160

    
161
            append = (append or jo.append_to_output)
162

    
163
            if jo.output_file is None:
164
                raise CAMAException("No output file specified")
165

    
166
            if jo.sensing_time[0] is None or jo.sensing_time[1] is None or dt is None:
167
                start = datetime.datetime.utcnow()
168
                data = Reader(joborder=jo, max_sza=max_sza)
169
                data.read()
170
                slice_written = False
171
                slice_written = data.dump(jo.output_file,
172
                                          mode=('a' if (slice_written or append) else 'w'),
173
                                          write_data=jo.config['include_data'])
174
                if not slice_written:
175
                    raise CAMAException("No data in requested interval")
176

    
177
                dtime = datetime.datetime.utcnow() - start
178
                logger.debug("Ingesting took {0} seconds".format(dtime.total_seconds()))
179

    
180
                action_idx = 1
181
                action_count = len(actions)
182

    
183
                for action, AnalysisClass in zip(['world', 'zonal', 'histogram', 'along_track', 'scatter', 'outline', 'events', 'irradiance'],
184
                                                 [WorldPlot, ZonalPlot, HistogramPlot, AlongTrackPlot, ScatterPlot, OutlinePlot, EventPlot, IrradiancePlot]):
185
                    if action in actions:
186
                        check_warning_counter()
187
                        logger.progress("Starting action {0} ({1}/{2})".format(action, action_idx, action_count))
188
                        action_idx += 1
189
                        processor = AnalysisClass(reader_data=data, resolution=spatial_resolution,
190
                                        number_bins=histogram_bin_count, zonal_resolution=zonal_resolution)
191
                        start = datetime.datetime.utcnow()
192
                        processor.process()
193
                        dtime = datetime.datetime.utcnow() - start
194
                        logger.debug("{1} process took {0} seconds".format(dtime.total_seconds(), action))
195
                        start = datetime.datetime.utcnow()
196

    
197
                        processor.dump(jo.output_file)
198

    
199
                        dtime = datetime.datetime.utcnow() - start
200
                        logger.debug("{1} writing took {0} seconds".format(dtime.total_seconds(), action))
201
            else:
202
                start = jo.sensing_time[0]
203
                end = jo.sensing_time[1]
204
                i = 0
205
                slice_written = False
206
                while start + (i+1)*dt <= end:
207
                    data = Reader(joborder=jo, interval=(start+i*dt, start+(i+1)*dt))
208
                    data.read() # add **kwargs to filter on attributes.
209
                    slice_written = data.dump(jo.output_file,
210
                                            mode=('a' if (slice_written or append) else 'w'),
211
                                            write_data=jo.config['include_data'])
212

    
213
                    if not slice_written:
214
                        i += 1
215
                        continue
216
                    action_idx = 1
217
                    action_count = len(actions)
218

    
219
                    for action, AnalysisClass in zip(['world', 'zonal', 'histogram', 'along_track', 'scatter', 'outline', 'events', 'irradiance'],
220
                                                    [WorldPlot, ZonalPlot, HistogramPlot, AlongTrackPlot, ScatterPlot, OutlinePlot, EventPlot, IrradiancePlot]):
221
                        if action in actions:
222
                            check_warning_counter()
223
                            logger.progress("Starting action {0} ({1}/{2})".format(action, action_idx, action_count))
224
                            action_idx += 1
225
                            processor = AnalysisClass(reader_data=data, resolution=spatial_resolution,
226
                                            number_bins=histogram_bin_count, zonal_resolution=zonal_resolution)
227
                            start = datetime.datetime.utcnow()
228
                            processor.process()
229
                            dtime = datetime.datetime.utcnow() - start
230
                            logger.debug("{1} process took {0} seconds".format(dtime.total_seconds(), action))
231
                            start = datetime.datetime.utcnow()
232

    
233
                            processor.dump(jo.output_file)
234

    
235
                            dtime = datetime.datetime.utcnow() - start
236
                            logger.debug("{1} writing took {0} seconds".format(dtime.total_seconds(), action))
237

    
238
                    i += 1
239
        else:
240
            if report_file is None:
241
                report_file = jo.report_file
242
            logger.info("Reading '{0}' and producing report without data extraction".format(jo.output_file))
243
            if 'actions' in jo.config and jo.processing_mode in jo.config['actions']:
244
                actions = jo.config['actions'][jo.processing_mode]
245
            else:
246
                logger.warning("Configuration error, no actions specified, using default list")
247
                if jo.processing_mode == 'NRTI':
248
                    actions = ['world', 'zonal', 'histogram', 'along_track',            'outline', 'events']
249
                else:
250
                    actions = ['world', 'zonal', 'histogram', 'along_track', 'scatter', 'outline', 'events']
251

    
252

    
253
        if report_file is not None:
254
            check_warning_counter()
255
            elapsed_time = datetime.datetime.utcnow() - start_time
256
            if not jo.report_only:
257
                logger.info("Data extraction time: {0:.1f} seconds".format(elapsed_time.total_seconds()))
258
            logger.info("Creating report '%s'", report_file)
259
            dest_directory=os.path.dirname(report_file)
260
            if dest_directory == "":
261
                dest_directory = os.path.curdir
262
            actions = set(actions).union(set(['table_stat', 'table_quantile', 'table_corr', 'table_cov', 'graph_corr', 'matrix_corr']))
263
            if 'events' in actions:
264
                actions.add('event_fig')
265
            if 'outline' in actions:
266
                actions.add('input')
267
            if 'world' in actions:
268
                actions.add('zonal_avg')
269

    
270
            with h5py.File(jo.output_file, 'r') as ref:
271
                time_range = [ref['time_bounds'][0, 0], ref['time_bounds'][-1, 1]]
272
                time_units = ref['time'].attrs['units']
273
                # logger.debug("time range: {0}, type {1}".format(time_range, str(type(time_range[0]))))
274
                times = netCDF4.num2date(time_range, time_units, only_use_cftime_datetimes=False, only_use_python_datetimes=True)
275
            duration_secs = (times[1] - times[0]).total_seconds()
276

    
277
            if math.fabs(duration_secs - 86400.0) < 6060.0:
278
                subject = "MPC report on {0} for {1[0]:%Y-%m-%d}".format(jo.product, times)
279
            elif duration_secs > 3*86400.0:
280
                subject = "MPC report on {0} for {1[0]:%Y-%m-%d} - {1[1]:%Y-%m-%d}".format(jo.product, times)
281
            else:
282
                subject = "MPC report on {0} for {1[0]:%Y-%m-%d (%H:%M:%S)} - {1[1]:%Y-%m-%d (%H:%M:%S)}".format(jo.product, times)
283

    
284
            title = "MPC report on {0}".format(jo.product)
285
            summary = "Processing for mode ``{0}''. ".format(jo.processing_mode)
286
            if 'Summary' in jo.dynamic_processing_parameters:
287
                summary = summary + jo.dynamic_processing_parameters['Summary']
288

    
289

    
290
            cama_plot(jo.output_file,
291
                      dest_directory=os.path.dirname(report_file),
292
                      output_file=os.path.basename(report_file),
293
                      output_type=os.path.splitext(report_file)[-1][1:],
294
                      title=title,
295
                      subject=subject,
296
                      summary=summary,
297
                      keywords=[jo.product, "Sentinel 5 precursor", "PyCAMA", "Mission performance center", "Copernicus"],
298
                      author=jo.station,
299
                      logger=logger,
300
                      add_minmax=False,
301
                      add_mean=True,
302
                      add_stddev=False,
303
                      actions=actions,
304
                      preliminary_unvalidated_data=jo.preliminary_unvalidated_data)
305
    except JobOrderError as err:
306
        rval = 255
307
    except CAMAException as err:
308
        if logger is None:
309
            print("Caught an exception, logging system is not yet set up", file=sys.stderr)
310
            print(err.message, file=sys.stderr)
311
            print(traceback.format_exc(), file=sys.stderr)
312
        else:
313
            logger.error(err.message)
314
        rval = err.return_value
315
    except KeyboardInterrupt:
316
        if logger is None:
317
            print("Keyboard interrupt, logging system is not yet set up", file=sys.stderr)
318
            print(str(err), file=sys.stderr)
319
        else:
320
            logger.error("Terminating because of keyboard interrupt")
321
        rval = 255
322
    except Exception as err:
323
        if logger is None:
324
            print("Caught an unexpected exception, logging system is not yet set up", file=sys.stderr)
325
            print(str(err), file=sys.stderr)
326
        else:
327
            logger.warning("Caught an unexpected exception")
328
            logger.error(str(err))
329
        print(traceback.format_exc(), file=sys.stderr)
330
        rval = 255
331
    else:
332
        if total_errors() > 0:
333
            rval = 255
334
        elif total_warnings() > 0:
335
            rval = 127
336
        else:
337
            rval = 0
338

    
339
    if logger is not None:
340
        elapsed_time = datetime.datetime.utcnow() - start_time
341
        logger.info("Processing time: {0:.1f} seconds".format(elapsed_time.total_seconds()))
342
        logger.info("Errors: {0}, warnings: {1}".format(total_errors(), total_warnings()))
343
        logger.info("Stopping with exit code: %d", rval)
344
        logging.shutdown()
345
    else:
346
        rval = 255
347
    sys.exit(rval)
348

    
349
## Write a default report for PyCAMA.
350
#
351
# The report includes (depending on settings)
352
# * world maps (including separate north- and southpole views),
353
# * histograms of the parameters (including tables with basic statistical parameters)
354
# * along track statistics
355
# * scatter density plots.
356
# * outline plots
357
# * tables with statistical analysis.
358
# * a graphical representation of the correlation matrix
359
#
360
# @param source_filename The netCDF4 file with the gathered input data.
361
# @param dest_directory Directory where the output should be written.
362
# @param output_file The output filename (created within the destination directory).
363
# @param output_type Possible values: 'tex' or 'pdf'. Type of generated report.
364
# @param title Metadata for pdf (also works on tex)
365
# @param subject Metadata for pdf (also works on tex)
366
# @param keywords Metadata for pdf (also works on tex) (list of values)
367
# @param author Metadata for pdf (also works on tex)
368
# @param logger The logger object
369
# @param add_minmax For the along track statistics plot: optionally include minimum & maximum, as
370
#                    these may compress the output to be nearly unreadable. Default is to include.
371
# @param add_mean For the along track statistics plot: optionally include the mean, as
372
#                    this may compress the output to be nearly unreadable. Default is to include.
373
# @param add_stddev For the along track statistics plot: optionally include the standard deviation, as
374
#                    this may compress the output to be nearly unreadable. Default is to include.
375
# @param actions sequence of types to include. Default = ('world', 'histogram', 'along_track', 'scatter', 'table_stat',
376
#                            'table_quantile', 'table_corr', 'table_cov').
377
# @param variable_names Only plot a subset of variables, default is to plot all.
378
# @param start First time slice to plot, index in file.
379
# @param end Last time-slices to plot, index in the file.
380
#               If end is set to -1, then only the last slice will be processed.
381
# @param figure_type Force the output type for the figures from PyCAMA (for use with
382
#               the website or for publications). The default 'dynamic' uses an
383
#               appropriate type for use with a pdf report. Legal values 'dynamic', 'png', 'pdf', 'svg'.
384
# @param add_legend Add legend to along track plots.
385
# @param preliminary_unvalidated_data Add label to worldplot to make clear that
386
#               the data shown in the plot is preliminary.
387
# Note that individual pages from the 'pdf' output can be accessed directly
388
# from LaTeX for publication. For delivery to a publisher you will need
389
# separate figures.
390
#
391
def cama_plot(source_filename, dest_directory=None, output_file=None, output_type='pdf',
392
              title=None, subject=None, summary=None, keywords=None, author=None, logger=None,
393
              add_minmax=True, add_mean=True, add_stddev=True,
394
              actions=None, variable_names=None, start=0, end=None, figure_type='dynamic',
395
              add_legend=False, preliminary_unvalidated_data=False):
396
    import matplotlib
397
    matplotlib.use('svg')
398
    matplotlib.rc('font', family='Helvetica Neue, Helvetica, DejaVu Sans, Bitstream Vera Sans, Lucida Sans, Lucida Grande, Verdana, Geneva, Lucid, Arial, Avant Garde, sans-serif')
399

    
400
    import matplotlib.pyplot as plt
401
    import pycama.cama_plot_time as cama_plot_time
402

    
403
    if logger is None:
404
        logger = setup_logging(loglevel="info", production_logger=True)
405

    
406
    if actions is None:
407
        actions = set('world', 'zonal', 'histogram', 'along_track', 'scatter', 'table_stat', 'events',
408
                   'table_quantile', 'table_corr', 'table_cov', 'graph_corr',
409
                   'matrix_corr', 'outline', 'input')
410

    
411
    if 'outline' in actions:
412
        actions.add('input')
413

    
414
    version_in_netcdf = read_version_number_from_file(source_filename)
415
    comparison_result = compare_versions(version_in_netcdf, level=2)
416
    if comparison_result > 0:
417
        raise CAMAException("The version in the netCDF file is made with a newer version of PyCAMA ({0} in file, {1} in software).".format(version_in_netcdf, __version__))
418
    elif comparison_result < 0:
419
        logger.warning("This version of PyCAMA is newer than the version of PyCAMA used to create the netCDF file ({0}). Proceed with caution".format(version_in_netcdf))
420

    
421
    if variable_names is None:
422
        vname_list = get_variable_list(source_filename)
423
    else:
424
        source_variable_names = get_variable_list(source_filename)
425
        vname_list = [v for v in variable_names if v in source_variable_names]
426
        if len(vname_list) == 0:
427
            logger.error("No valid variable_names given")
428
            return
429

    
430
    report = Report(directory=dest_directory, output_file=output_file, output_type=output_type)
431
    report.add_metadata(title=title, subject=subject, keywords=keywords, author=author, summary=summary)
432

    
433
    numsteps = count_time_steps(source_filename)
434
    if end is not None:
435
        if end < 0:
436
            start = numsteps-1
437
        else:
438
            numsteps = end if (0 < end < numsteps) else numsteps
439
    else:
440
        start = start if (0<= start < numsteps) else 0
441

    
442
    for i in range(start, numsteps):
443
        check_warning_counter()
444
        t, t_start, t_end = extract_time(source_filename, i)
445

    
446
        if 'outline' in actions:
447
            ou = OutlinePlot(filename=source_filename, time_index=i)
448
            if not ou.success:
449
                ou = None
450
        else:
451
            ou = None
452

    
453
        if ou is not None:
454
            name = 'outline_granules'
455
            figure = report.start_figure("outline", name, t_start, t_end, figsize=Report.figure_size(18,10), output_type='png')
456
            if ou.plot(figure, ax=None, resolution='l'):
457
                report.end_figure(figure, dpi=300)
458
            else:
459
                report.clear_figure(figure)
460
            plt.close('all')
461

    
462
        if 'input' in actions:
463
            data = cama_plot_time.outline(logger, input_files=[source_filename])
464
            figure = report.start_figure("inputs", 'inputs', t_start, t_end, figsize=Report.figure_size(18,18), output_type='pdf' if figure_type == 'dynamic' else figure_type)
465
            cama_plot_time.outline_plot(logger, figure, data)
466
            report.end_figure(figure, dpi=300)
467
        else:
468
            logger.warning("Not producing 'inputs' plot")
469

    
470
        if 'events' in actions:
471
            data = cama_plot_time.events(logger, input_files=[source_filename], error=False, filters=False, general=False, warning=False, group=False)
472
            figure = report.start_figure("events", 'events', t_start, t_end, figsize=Report.figure_size(18,18), output_type='pdf' if figure_type == 'dynamic' else figure_type)
473
            cama_plot_time.events_plot(logger, figure, data, error=False, filters=False, general=False, warning=False, group=False, fraction=True)
474
            report.end_figure(figure, dpi=300)
475
        else:
476
            logger.warning("Not producing 'events' plot")
477

    
478
        if 'zonal' in actions:
479
            zp = ZonalPlot(filename=source_filename, time_index=i)
480
            if not zp.success:
481
                zp = None
482
        else:
483
            zp = None
484

    
485
        if 'world' in actions or ('zonal' in actions and zp is None):
486
            wp = WorldPlot(filename=source_filename, time_index=i)
487
            if not wp.success:
488
                wp = None
489
        else:
490
            wp = None
491

    
492
        if wp is not None and 'world' in actions:
493
            for name in vname_list:
494
                check_warning_counter()
495
                if name not in wp.variables_meta:
496
                    continue
497
                caption = "Map of ``{0}'' for {1:%Y-%m-%d} to {2:%Y-%m-%d}".format(wp.variables_meta[name]['title'], t_start, t_end)
498
                figure = report.start_figure("world", name, t_start, t_end,
499
                        figsize=Report.figure_size(18,20),
500
                        output_type='png' if figure_type == 'dynamic' else figure_type,
501
                        caption=caption)
502

    
503
                rval = wp.plot(name, figure, ax_idx=(2,1,1), projection='world',
504
                               resolution='l', add_colorbar=True, linewidth=0.5,
505
                               time=t, preliminary_unvalidated_data=preliminary_unvalidated_data)
506
                rval = rval and wp.plot(name, figure, ax_idx=(2,2,3),
507
                                        projection='north pole', resolution='l',
508
                                        add_colorbar=False, linewidth=0.5,
509
                                        preliminary_unvalidated_data=preliminary_unvalidated_data)
510
                rval = rval and wp.plot(name, figure, ax_idx=(2,2,4),
511
                                        projection='south pole', resolution='l',
512
                                        add_colorbar=False, linewidth=0.5,
513
                                        preliminary_unvalidated_data=preliminary_unvalidated_data)
514

    
515
                plt.tight_layout()
516
                if rval:
517
                    report.end_figure(figure, dpi=300)
518
                else:
519
                    report.clear_figure(figure)
520

    
521
                plt.close('all')
522

    
523
            for name in vname_list:
524
                if name in wp.variables_meta:
525
                    break
526
            check_warning_counter()
527
            caption = "Map of the number of observations for {0:%Y-%m-%d} to {1:%Y-%m-%d}".format(t_start, t_end)
528
            figure = report.start_figure("world", name+'_count', t_start, t_end,
529
                    figsize=Report.figure_size(18,18),
530
                    output_type='png' if figure_type == 'dynamic' else figure_type,
531
                    caption=caption)
532
            rval = wp.plot(name, figure, ax_idx=(2,1,1), projection='world', resolution='l', add_colorbar=True, linewidth=0.5, time=t, count=True)
533
            rval = rval and wp.plot(name, figure, ax_idx=(2,2,3), projection='north pole', resolution='l', add_colorbar=False, linewidth=0.5, count=True)
534
            rval = rval and wp.plot(name, figure, ax_idx=(2,2,4), projection='south pole', resolution='l', add_colorbar=False, linewidth=0.5, count=True)
535
            plt.tight_layout()
536
            if rval:
537
                report.end_figure(figure, dpi=300)
538
            else:
539
                report.clear_figure(figure)
540

    
541
            plt.close('all')
542

    
543
        if 'zonal' in actions and (wp is not None or zp is not None):
544
            if zp is None:
545
                zp = wp
546
            for name in vname_list:
547
                check_warning_counter()
548
                if name not in zp.variables_meta:
549
                    continue
550
                caption = "Zonal average of ``{0}'' for {1:%Y-%m-%d} to {2:%Y-%m-%d}.".format(zp.variables_meta[name]['title'], t_start, t_end)
551
                figure = report.start_figure("zonal_avg", name, t_start, t_end,
552
                        figsize=Report.figure_size(18,18),
553
                        output_type='pdf' if figure_type == 'dynamic' else figure_type,
554
                        caption=caption)
555

    
556
                rval = zp.plot_zonal_mean(name, figure, ax=plt.subplot(1,1,1), autorange=True)
557
                if rval:
558
                    report.end_figure(figure, dpi=300)
559
                else:
560
                    report.clear_figure(figure)
561

    
562
                plt.close('all')
563

    
564
        if 'histogram' in actions:
565
            hp = HistogramPlot(filename=source_filename, time_index=i)
566
            if not hp.success:
567
                hp = None
568
        else:
569
            hp = None
570

    
571
        if hp is not None:
572
            for name in vname_list:
573
                check_warning_counter()
574
                if name not in hp.variables_meta:
575
                    continue
576
                caption = "Histogram of ``{0}'' for {1:%Y-%m-%d} to {2:%Y-%m-%d}".format(hp.variables_meta[name]['title'], t_start, t_end)
577
                figure = report.start_figure("histogram", name, t_start, t_end,
578
                        figsize=Report.figure_size(18,18),
579
                        output_type='pdf' if figure_type == 'dynamic' else figure_type,
580
                        caption=caption)
581

    
582
                rval = hp.plot(name, figure, ax=plt.subplot(211), weighted=True, time=t, clip=True, legend=True, ylog=False, add_xlabel=False)
583
                rval = rval and hp.plot(name, figure, ax=plt.subplot(212), weighted=False, time=t, legend=False, clip=True, ylog=False)
584

    
585
                # plt.tight_layout()
586
                if rval:
587
                    report.end_figure(figure, dpi=120, apply_tight_layout=False)
588
                else:
589
                    report.clear_figure(figure)
590

    
591
            if report.istex:
592
                if 'table_stat' in actions or 'table_quantile' in actions:
593
                    report.add_description(hp.report())
594
            else:
595
                if 'table_stat' in actions:
596
                    figure = report.start_figure("parameters", 'base_stats', t_start, t_end,
597
                            figsize=Report.figure_size(30,10),
598
                            output_type='png' if figure_type == 'dynamic' else figure_type)
599
                    table = hp.table(variant='base')
600
                    report.end_figure(figure, dpi=300, bbox_inches='tight', bbox_extra_artists=table)
601

    
602
                if 'table_quantile' in actions:
603
                    figure = report.start_figure("parameters", 'extra_stats', t_start, t_end,
604
                            figsize=Report.figure_size(30,10),
605
                            output_type='png' if figure_type == 'dynamic' else figure_type)
606
                    table = hp.table(variant='quantile')
607
                    report.end_figure(figure, dpi=300, bbox_inches='tight', bbox_extra_artists=table)
608
                plt.close('all')
609

    
610
        if 'along_track' in actions:
611
            at = AlongTrackPlot(filename=source_filename, time_index=i)
612
            if not at.success:
613
                at = None
614
        else:
615
            at = None
616

    
617
        if at is not None:
618
            for j, name in enumerate(vname_list):
619
                check_warning_counter()
620
                if name not in at.variables_meta:
621
                    continue
622
                caption = "Along track statistics of ``{0}'' for {1:%Y-%m-%d} to {2:%Y-%m-%d}".format(at.variables_meta[name]['title'], t_start, t_end)
623
                figure = report.start_figure("alongtrack", name, t_start, t_end,
624
                        figsize=Report.figure_size(18,18),
625
                        output_type='pdf' if figure_type == 'dynamic' else figure_type,
626
                        caption=caption)
627

    
628
                if at.plot(name, figure, time=t, add_legend=add_legend,
629
                           add_minmax=add_minmax, add_mean=add_mean, add_stddev=add_stddev):
630
                    # do not use tight layout, this will mess up the spacing of the legend
631
                    report.end_figure(figure, dpi=120, apply_tight_layout=(not add_legend))
632
                else:
633
                    report.clear_figure(figure)
634
                plt.close('all')
635

    
636
        if 'scatter' in actions:
637
            sp = ScatterPlot(filename=source_filename, time_index=i)
638
            if not sp.success:
639
                sp = None
640
        else:
641
            sp = None
642

    
643
        if sp is not None:
644
            for name1, name2 in sp.key_pairs:
645
                check_warning_counter()
646
                if variable_names is not None and ((name1 not in vname_list) or (name2 not in vname_list)):
647
                    continue
648
                caption="Scatter density plot of ``{0}'' against ``{1}'' for {2:%Y-%m-%d} to {3:%Y-%m-%d}.".format(sp.variables_meta[name1]['title'], sp.variables_meta[name2]['title'], t_start, t_end)
649
                figure = report.start_figure("scatter", name1, t_start, t_end,
650
                        secondary_variable=name2,
651
                        figsize=Report.figure_size(18,24),
652
                        output_type='png' if figure_type == 'dynamic' else figure_type,
653
                        caption=caption)
654

    
655
                rval = sp.plot((name1, name2), figure, ax=plt.subplot(211), log=False, time=t, colorscale='viridis', hide_R=False)
656
                rval = rval and sp.plot((name1, name2), figure, ax=plt.subplot(212), log=True, time=t, colorscale='viridis', hide_R=True)
657

    
658
                plt.tight_layout()
659
                if rval:
660
                    report.end_figure(figure, dpi=300)
661
                else:
662
                    report.clear_figure(figure)
663
                plt.close('all')
664

    
665
            if report.istex:
666
                if 'table_corr' in actions and len(sp.variable_names) < 45:
667
                    report.add_matrix(sp.report(covariance=False))
668
                if 'table_cov' in actions and len(sp.variable_names) < 45:
669
                    report.add_matrix(sp.report(covariance=True))
670
                if 'graph_corr' in actions:
671
                    figure = report.start_figure("graph", 'correlation_graph', t_start, t_end,
672
                            figsize=Report.figure_size(18,18),
673
                            output_type='pdf' if figure_type == 'dynamic' else figure_type)
674
                    if sp.correlation_matrix_plot():
675
                        report.end_figure(figure, dpi=300, bbox_inches='tight')
676
                    else:
677
                        report.clear_figure(figure)
678
                if 'matrix_corr' in actions:
679
                    figure = report.start_figure("graph", 'correlation_matrix', t_start, t_end,
680
                            figsize=Report.figure_size(18,18),
681
                            output_type='pdf' if figure_type == 'dynamic' else figure_type, )
682
                    if sp.correlation_graph():
683
                        report.end_figure(figure, dpi=300, bbox_inches='tight')
684
                    else:
685
                        report.clear_figure(figure)
686
            else:
687
                if 'table_corr' in actions:
688
                    figure = report.start_figure("matrix", 'correlation', t_start, t_end,
689
                            figsize=Report.figure_size(30,10),
690
                            output_type='pdf' if figure_type == 'dynamic' else figure_type)
691
                    table = sp.table(covariance=False)
692
                    report.end_figure(figure, dpi=300, bbox_inches='tight', bbox_extra_artists=table)
693
                if 'table_cov' in actions:
694
                    figure = report.start_figure("matrix", 'covariance', t_start, t_end,
695
                            figsize=Report.figure_size(30,10),
696
                            output_type='pdf' if figure_type == 'dynamic' else figure_type, )
697
                    table = sp.table(covariance=True)
698
                    report.end_figure(figure, dpi=300, bbox_inches='tight', bbox_extra_artists=table)
699
                if 'graph_corr' in actions:
700
                    figure = report.start_figure("graph", 'correlation_graph', t_start, t_end,
701
                            figsize=Report.figure_size(10,10),
702
                            output_type='pdf' if figure_type == 'dynamic' else figure_type)
703
                    if sp.correlation_matrix_plot():
704
                        report.end_figure(figure, dpi=300, bbox_inches='tight')
705
                    else:
706
                        report.clear_figure(figure)
707
                if 'matrix_corr' in actions:
708
                    figure = report.start_figure("graph", 'correlation_matrix', t_start, t_end,
709
                            figsize=Report.figure_size(10,10),
710
                            output_type='pdf' if figure_type == 'dynamic' else figure_type)
711
                    if sp.correlation_graph():
712
                        report.end_figure(figure, dpi=300, bbox_inches='tight')
713
                    else:
714
                        report.clear_figure(figure)
715
            plt.close('all')
716
    report.close()
717

    
718
## command line interface for the data extractor.
719
if __name__ == "__main__":
720
    import argparse
721

    
722
    descr = """Extract data from Sentinel 5 precursor level 2 files and produce a
723
summary of the contents.
724

725
    This is the data extractor for the mission performance centre (MPC) for Sentinel 5 precursor. In the PDGS the interface follows the thin layer interface.
726

727
    A description of the job order file and the configuration file can be found in the software user manual [MPC-KNMI-CC-0014-MA]. """
728
    parser = argparse.ArgumentParser(description=descr,
729
                                     prog='PyCAMA')
730
    parser.add_argument('--version', action='version',
731
                        version='%(prog)s {0}'.format(__version__))
732
    parser.add_argument('joborder', metavar='FILE',
733
                        type=argparse.FileType('r'),
734
                        help='The joborder file')
735
    parser.add_argument('--aggregation', metavar='PERIOD',
736
                        help='Aggregation interval, specified as an ISO 8601 duration.')
737
    parser.add_argument('--actions', nargs="+", default=None,
738
                        choices=('world', 'zonal', 'histogram', 'along_track', 'scatter',
739
                                 'outline', 'events', 'irradiance', 'input'),
740
                        metavar="ACTION",
741
                        help="Explicitly specify the atcions to perform. Defaults to actions in configuration file.")
742
    parser.add_argument('--append', action='store_true',
743
                        help='Append to existing file instead of creating a new file.')
744
    parser.add_argument('--report', metavar='REPORT',
745
                        help='Where to write a pdf or latex output report.')
746
    parser.add_argument('--report-only', action='store_true',
747
                        help='Create a report based on a previous data extraction run.')
748
    parser.add_argument('--limit-warnings', type=int, default=None,
749
                        help="Limit the number of warnings before bailing out with an error.")
750
    parser.add_argument('--msg', choices=('DEBUG', 'INFO', 'WARNING', 'ERROR'), help='Logging message level.')
751

    
752
    args = parser.parse_args()
753
    main(args.joborder, aggregation=args.aggregation,
754
         append=args.append, report_file=args.report,
755
         report_only=args.report_only, warnings_limit=args.limit_warnings,
756
         actions_argument=args.actions, message_level=args.msg)