Project

General

Profile

Pycama-4-O3_TCL.py

Maarten Sneep, 10/30/2019 04:52 PM

Download (19.5 KB)

 
1
# -*- coding: utf-8 -*-
2
"""
3
genertae overview plots for O3_TCL
4
and histogramms as used in Pycma
5

6
Call as "Pycama-4-O3_TCL.py -i S5P_OFFL_L2__O3_TCl_*.nc -o "xy.pdf"
7

8
@author: heue_kl
9
"""
10

    
11
# reorder the loading of modules: first teh standard library
12
import os
13
import sys
14
import argparse
15
from datetime import datetime
16

    
17
# external packages
18
import numpy as np
19
from matplotlib.backends.backend_pdf import PdfPages
20
import h5py
21

    
22
# own modules
23
import tropoO3ResultsCSA as tropoO3Results
24

    
25
def parseCmdLine():
26
    parser = argparse.ArgumentParser( description="""Generator to create tropospheric O3 products""",
27
            formatter_class=argparse.ArgumentDefaultsHelpFormatter )
28

    
29
    parser.add_argument( "-o", "--output", required=True, help="PDF output for Pycama_report" )
30
    parser.add_argument( "-i", "--input", dest='input_file', required=True, help="O3_TCL Level2C file")
31
    args = parser.parse_args()
32

    
33
    return args
34

    
35

    
36

    
37

    
38
def readO3_TCL(filename):
39
    #from netCDF4 import Dataset ## Do not use NetCDF4 tries to multiply with scale_factor ='1e-9' which is a string
40

    
41
    #MS: This is a bug in the O3_TCL processor. Following the standards, this constant should have the same data-type
42
    #MS: as the target type for a "packed" variable.
43
    #MS: see http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#packed-data
44
    #MS: Scaling can be turned off using nc_file.set_auto_mask(False) if the netCDF4 module is used to access the file.
45

    
46
    import h5py
47
    nc_file=h5py.File(filename,'r')
48

    
49
    with nc_file:
50
        #lat=nc_file.groups['PRODUCT'].variables['latitude'][:]
51
        start = str(nc_file.attrs.get('time_coverage_troposphere_start'))
52
        end   = str(nc_file.attrs.get('time_coverage_troposphere_end'))
53

    
54
        lat = nc_file['PRODUCT/latitude_ccd'][:]
55
        lon = nc_file['PRODUCT/longitude_ccd'][:]
56

    
57
        latcoarse = nc_file['PRODUCT/latitude_csa'][:]
58
        loncoarse = nc_file['PRODUCT/longitude_csa'][:]
59

    
60
        O3_TCL = nc_file['PRODUCT/ozone_tropospheric_vertical_column'][0,:,:]
61
        O3_TCL_sig = nc_file['PRODUCT/ozone_tropospheric_vertical_column_precision'][0,:,:]
62
        O3_TCL_fv = nc_file['PRODUCT/ozone_tropospheric_vertical_column'].attrs.get('_FillValue')[0]
63

    
64
        O3_TCL_strat_ref = nc_file['PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/ozone_stratospheric_vertical_column_reference'][0,:]
65
        O3_TCL_strat_flag = nc_file['PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/ozone_stratospheric_vertical_column_reference_flag'][0,:]
66

    
67
        O3_TCL_nr = nc_file['PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/number_of_observations_ozone_tropospheric_vertical_column'][0,:,:]
68
        O3_TCL_nr_fv = nc_file['PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/number_of_observations_ozone_tropospheric_vertical_column'].attrs.get('_FillValue')[0]
69

    
70
        # recalcualte rather than reading
71
        #    O3_TCL_hist=nc_file['METADATA/QA_STATISTICS/ozone_tropospheric_vertical_column_histogram'][:]
72
        #    O3_TCL_hist_axis=nc_file['METADATA/QA_STATISTICS/histogram_axis_tropospheric_ozone'][:]
73

    
74
        O3_MRu = nc_file['PRODUCT/ozone_upper_tropospheric_mixing_ratio'][0,:,:]
75
        O3_MRu_sig = nc_file['PRODUCT/ozone_upper_tropospheric_mixing_ratio_precision'][0,:,:]
76
        O3_MRu_fv = nc_file['PRODUCT/ozone_upper_tropospheric_mixing_ratio'].attrs.get('_FillValue')[0]
77
        O3_MRu_nr = nc_file['PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/number_of_observations_ozone_upper_tropospheric_mixing_ratio'][0,:,:]
78
        O3_MRu_flag = nc_file['PRODUCT/ozone_upper_tropospheric_mixing_ratio_flag'][0,:,:]
79

    
80
        # additional data to be read?
81
        # ozone_stratospheric_vertical_column_reference?
82
        # ozone_stratospheric_vertical_column_reference_flag?
83
        # ozone_stratospheric_vertical_column_reference_precision?
84

    
85
    f_dlat = abs(lat[1]-lat[0])
86
    f_latmin = np.min(lat)-0.5*f_dlat
87
    f_latmax = np.max(lat)+0.5*f_dlat
88

    
89
    f_dlong = abs( lon[1]-lon[0])
90
    f_longmin = np.min(lon)-0.5*f_dlong
91
    f_longmax = np.max(lon)+0.5*f_dlong
92

    
93
    f_dlatc = abs(latcoarse[1]-latcoarse[0])
94
    f_dlongc = abs(loncoarse[1]-loncoarse[0])
95

    
96
    results = tropoO3Results.TropoO3ResultsCSA(f_latmin, f_dlat, f_dlatc, f_latmax, f_longmin,f_dlong, f_dlongc, f_longmax)
97

    
98
    results._tropoO3 = np.ma.masked_values(O3_TCL, O3_TCL_fv)
99
    results._tropoO3_std = np.ma.masked_values(O3_TCL_sig, O3_TCL_fv)
100
    results._tropoO3_count = np.ma.masked_values(O3_TCL_nr, O3_TCL_nr_fv)
101

    
102
    results._stratO3Reference = np.ma.masked_values(O3_TCL_strat_ref, O3_TCL_fv)
103
    results._stratO3Reference_flag = O3_TCL_strat_flag
104

    
105
    results._upper_tropoO3_mixingratio = np.ma.masked_values(O3_MRu, O3_MRu_fv)
106
    results._upper_tropoO3_mixingratio_std = np.ma.masked_values(O3_MRu_sig, O3_MRu_fv)
107
    results._upper_tropoO3_count = np.ma.masked_values(O3_MRu_nr, O3_TCL_nr_fv)
108
    results._upper_tropoO3_flag = O3_MRu_flag
109

    
110
    start = start.replace('b','')
111
    end   = end.replace('b','')
112
    results._centralStart = start.split('T')[0]
113
    results._centralEnd = end.split('T')[0]
114

    
115
    return(results)
116

    
117

    
118
def worldPlots(results, out=None):
119
    """plot the data and the deviation in the netCDF on a worldmap and store the figue as postscript if filename is given"""
120

    
121
    import matplotlib.pyplot as plt
122
    from matplotlib.ticker import MaxNLocator
123
    from matplotlib.colors import BoundaryNorm
124
    import cartopy.crs as ccrs
125

    
126
    fa_lat = results._latCenters
127
    fa_long = results._lonCenters
128

    
129
    fa_latC = results._latCentersCoarse
130
    fa_longC = results._lonCentersCoarse
131

    
132
    f_dlat = abs(fa_lat[1]-fa_lat[0])
133
    f_dlong = abs(fa_long[1]-fa_long[0])
134

    
135
    f_dlatC = abs(fa_latC[1]-fa_latC[0])
136
    f_dlongC = abs(fa_longC[1]-fa_longC[0])
137

    
138
    fa_VCDtrop = results._tropoO3
139
    ia_VCDtropnr = results._tropoO3_count
140
    #    fa_UtropMR = results._upper_tropoO3_mixingratio
141

    
142
    fa_VCDtrop = np.ma.masked_invalid(fa_VCDtrop)
143
    fa_VCDtrop = np.ma.masked_values(fa_VCDtrop, 0)
144

    
145
    #    ia_UtropMRnr=results._upper_tropoO3_count
146
    #    fa_UtropMR=np.ma.masked_invalid(fa_UtropMR)
147
    #    fa_UtropMR=np.ma.masked_values(fa_UtropMR,0)
148

    
149
    ## for the meshgrid not the centres but the corners should be used
150

    
151
    fa_lat = np.append(fa_lat-0.5*f_dlat, fa_lat[-1]+0.5*f_dlat)
152
    fa_long = np.append(fa_long-0.5*f_dlong, fa_long[-1]+0.5*f_dlong)
153

    
154
    fa_latC = np.append(fa_latC-0.5*f_dlatC, fa_latC[-1]+0.5*f_dlatC)
155
    fa_longC = np.append(fa_longC-0.5*f_dlongC, fa_longC[-1]+0.5*f_dlongC)
156

    
157

    
158
    f_latmin = min(fa_lat)
159
    f_latmax = max(fa_lat)
160

    
161
    f_longmin = min(fa_long)
162
    f_longmax = max(fa_long)
163

    
164
    plt.close('all')
165

    
166
    fig,(trop, MR) = plt.subplots(2,1) #,sharex='all',sharey=True) ,strat,cldhgt,cldfr
167

    
168
    fig.set_size_inches(20, 13)
169
    fig.subplots_adjust(bottom=0.2, right=0.9,hspace=0.0)
170

    
171
    x, y = np.meshgrid(fa_long, fa_lat)
172
    xc, yc = np.meshgrid(fa_longC, fa_latC)
173

    
174
    cm = plt.get_cmap('jet')       #gist_rainbow_r
175
    cm.set_bad('white',1)
176
    cm.set_under('indigo',1)
177

    
178
    # tropopospheric VCD
179
    trop = plt.subplot(2, 1, 1, projection=ccrs.PlateCarree())
180
    trop.set_extent([-180, 180, -22, 22]) #MS: shouldn't that come from the f_latmin and f_latmax variables?
181

    
182
    trop.coastlines()
183
    trop.gridlines(ylocs=([-30, -20, -10, 0, 10, 20, 30]), xlocs=([-180, -90, 0, 90, 180]))
184

    
185
    levels = MaxNLocator(nbins=30).tick_values(0,30)
186
    norm = BoundaryNorm(levels, ncolors=cm.N, clip=True)
187
    cs = trop.pcolor(x, y, fa_VCDtrop*1000, cmap=cm, norm=norm, transform=ccrs.PlateCarree())
188
    cbar = plt.colorbar(cs, ax=trop, shrink=0.6, pad=0.11, aspect=35, orientation="horizontal", extendfrac='auto', extend='both')
189
    cbar.set_label("mmol/m²", fontsize=16, labelpad=10)
190
    cbar.set_ticks(np.arange(0, 31,5))
191
    cbar.ax.tick_params(labelsize=16)
192
    s_title='Tropospheric Ozone Column Density'
193
    plt.title(s_title, fontsize=16, horizontalalignment='left', x=0.05)
194

    
195
    MR = plt.subplot(2, 1, 2, projection=ccrs.PlateCarree())
196
    MR.set_extent([-180, 180, -22, 22]) #MS: shouldn't that come from the f_latmin and f_latmax variables?
197

    
198
    # draw coastlines
199
    MR.coastlines()
200
    MR.gridlines(ylocs=(np.arange(-30, 30, 10)), xlocs=(-180, -90, 0, 90, 180))
201

    
202
    #    levels=MaxNLocator(nbins=50).tick_values(0,100)
203
    #    norm=BoundaryNorm(levels,ncolors=cm.N, clip=True)
204
    ##    cs=MR.pcolor(xc,yc,fa_UtropMR,cmap=cm,norm=norm,transform=ccrs.PlateCarree())
205
    #    cbar=plt.colorbar(cs,ax=MR,shrink=0.6,pad=0.11,aspect=35, orientation="horizontal",extendfrac='auto',extend='both')
206
    #    cbar.set_label("ppb",fontsize=16,labelpad=10)
207
    #    cbar.set_ticks(np.arange(0,101,10))
208
    #    cbar.ax.tick_params(labelsize=16)
209
    s_title='Upper Tropospheric Ozone Mixing Ratio'
210
    plt.title(s_title, fontsize=16, horizontalalignment='left', x=0.05)
211

    
212
    # print(out)
213

    
214
    out.savefig()
215
    plt.close()
216

    
217
    fig2,(tropnr,MRnr) = plt.subplots(2, 1) #,sharex='all',sharey=True) ,strat,cldhgt,cldfr
218

    
219
    fig2.set_size_inches(20, 13)
220
    fig2.subplots_adjust(bottom=0.2, right=0.9, hspace=0.0)
221

    
222
    x, y = np.meshgrid(fa_long, fa_lat)
223
    xc, yc = np.meshgrid(fa_longC, fa_latC)
224

    
225
    cm = plt.get_cmap('jet')       #gist_rainbow_r
226
    cm.set_bad('white',1)
227
    cm.set_under('indigo',1)
228

    
229
    # tropopospheric number
230
    trop=plt.subplot(2, 1, 1, projection=ccrs.PlateCarree())
231
    trop.set_extent([-180, 180, -22, 22]) #MS: shouldn't that come from the f_latmin and f_latmax variables?
232

    
233
    trop.coastlines()
234
    trop.gridlines(ylocs=(np.arange(-30, 30, 10)), xlocs=(-180, -90, 0, 90, 180))
235

    
236
    levels = MaxNLocator(nbins=20).tick_values(0, 600)
237
    norm = BoundaryNorm(levels, ncolors=cm.N, clip=True)
238
    cs = trop.pcolor(x, y, ia_VCDtropnr, cmap=cm, norm=norm, transform=ccrs.PlateCarree())
239
    cbar = plt.colorbar(cs, ax=trop, shrink=0.6, pad=0.11, aspect=35, orientation="horizontal", extendfrac='auto', extend='both')
240
    cbar.set_label("", fontsize=16, labelpad=-5)
241
    cbar.set_ticks(np.arange(0, 601, 100))
242
    cbar.ax.tick_params(labelsize=16)
243
    s_title='number of observations in Tropospheric Ozone Column Density'
244
    plt.title(s_title, fontsize=16, horizontalalignment='left', x=0.05)
245

    
246

    
247
    MR=plt.subplot(2, 1, 2, projection=ccrs.PlateCarree())
248
    MR.set_extent([-180, 180, -22, 22])
249

    
250
    # draw coastlines
251
    MR.coastlines()
252
    MR.gridlines(ylocs=(np.arange(-30, 30, 10)),xlocs=(-180, -90, 0, 90, 180))
253

    
254
    #    levels = MaxNLocator(nbins=50).tick_values(0, 10000)
255
    #    norm = BoundaryNorm(levels, ncolors=cm.N, clip=True)
256
    #    cs = MR.pcolor(x, y, ia_UtropMRnr, cmap=cm, norm=norm, transform=ccrs.PlateCarree())
257
    #    cbar = plt.colorbar(cs, ax=MR, shrink=0.6, pad=0.11, aspect=35, orientation="horizontal", extendfrac='auto', extend='both')
258
    #    cbar.set_label("", fontsize=16, labelpad=-5)
259
    #    cbar.set_ticks(np.arange(0, 10000, 1000))
260
    #    cbar.ax.tick_params(labelsize=16)
261
    s_title='number of observations in Upper Tropospheric Ozone Mixing Ratio'
262
    plt.title(s_title, fontsize=16, horizontalalignment='left', x=0.05) #MS: shouldn't that come from the f_latmin and f_latmax variables?
263

    
264
    if out is None:
265
        plt.show()
266
    else:
267
        out.savefig()#,bbox_inches='tight')
268
    plt.close()
269

    
270
    return
271

    
272
def Histogramplots(results, out=None):
273

    
274
    import matplotlib.pyplot as plt
275
    f_max_vcdtrop = 3e-2
276
    f_max_umrtrop = 100
277
    i_vcdbins = 100
278
    i_umrbins = 50
279

    
280
    fa_VCDtrop = results._tropoO3
281
    fa_UtropMR = results._upper_tropoO3_mixingratio
282

    
283
    fa_VCDtrop = np.ma.masked_invalid(fa_VCDtrop)
284
    fa_VCDtrop = np.ma.masked_values(fa_VCDtrop,0)
285

    
286
    #    fa_UtropMR=np.ma.masked_invalid(fa_UtropMR)
287
    #    fa_UtropMR=np.ma.masked_values(fa_UtropMR,0)
288

    
289
    fa_histo3 = np.ma.compressed(fa_VCDtrop)
290
    i_higher = np.ma.count(np.ma.masked_less(fa_histo3,f_max_vcdtrop))
291

    
292
    #    fa_histuo3=np.ma.compressed(fa_UtropMR)
293
    #    i_uhigher=np.ma.count(np.ma.masked_less(fa_histuo3,f_max_umrtrop))
294

    
295
    plt.close('all')
296

    
297
    fig, (histVCD, histUMR) = plt.subplots(2, 1)
298
    fig.set_size_inches(20, 13)
299
    plt.subplots_adjust()
300

    
301
    histVCD = plt.subplot(2, 1, 1)
302
    histVCD = plt.hist(fa_histo3, bins=i_vcdbins, range=(0, f_max_vcdtrop))
303
    vcdtext = 'number of data higher than {0}: {1}'.format(f_max_vcdtrop, i_higher)
304
    plt.text(0.028, 800, vcdtext, horizontalalignment='right', fontsize=16)
305
    plt.axis([0, f_max_vcdtrop, 0, 1200])
306
    plt.xticks(fontsize=16)
307
    plt.yticks(fontsize=16)
308
    plt.xlabel('tropospheric ozone column (mol/m²)', fontsize=16)
309

    
310
    histUMR = plt.subplot(2, 1, 2)
311
    #    histUMR = plt.hist(fa_histuo3, bins=i_umrbins, range=(0,f_max_umrtrop))
312
    #    utext='number of data higher than {0}: {1}'.format(f_max_umrtrop, i_uhigher)
313
    #    plt.text(90, 12, utext, horizontalalignment='right')
314
    plt.axis([0, f_max_umrtrop, 0, 15])
315
    plt.xticks(fontsize=16)
316
    plt.xlabel('upper tropospheric mixing ratio (ppb)',fontsize=16)
317

    
318
    if out is None:
319
        plt.show()
320
    else:
321
        out.savefig()#,bbox_inches='tight')
322

    
323
    plt.close()
324

    
325
    return
326

    
327
def statistics(results, out=None):
328
    import matplotlib.pyplot as plt
329
    #from pylatex import Documentv  # geht nicht
330
    fa_VCDtrop  =results._tropoO3*1000
331
    fa_VCDtrop_std = results._tropoO3_std*1000
332

    
333
    fa_UtropMR = results._upper_tropoO3_mixingratio
334
    fa_UtropMR_std = results._upper_tropoO3_mixingratio_std
335

    
336
    f_meanVCD = np.ma.mean(fa_VCDtrop)
337
    f_medianVCD = np.ma.median(fa_VCDtrop)
338
    f_stdVCD = np.ma.std(fa_VCDtrop)
339
    f_minVCD = np.ma.min(fa_VCDtrop)
340
    f_maxVCD = np.ma.max(fa_VCDtrop)
341

    
342
    f_meansigVCD = np.ma.mean(fa_VCDtrop_std)
343

    
344
    f_meanUMR = np.ma.mean(fa_UtropMR)
345
    f_medianUMR = np.ma.median(fa_UtropMR)
346
    f_stdUMR = np.ma.std(fa_UtropMR)
347
    f_minUMR = np.ma.min(fa_UtropMR)
348
    f_maxUMR = np.ma.max(fa_UtropMR)
349

    
350
    f_meansigUMR = np.ma.mean(fa_UtropMR_std)
351

    
352
    StartDate = str(results._centralStart)[1:]
353
    EndDate   = str(results._centralEnd)[1:]
354

    
355
    fig = plt.figure(figsize=(20, 13))
356
    fig.clf()
357
    txt = 'Tropospheric Ozone Column monitoring analysis'
358
    fig.text(0.5, 0.85, txt, transform=fig.transFigure, size=24, ha="center")
359
    txt = 'Date range: {0} - {1}'.format(StartDate, EndDate)
360
    fig.text(0.5, 0.80, txt, transform=fig.transFigure, size=14, ha="center")
361
    txt = 'Statistics: tropospheric ozone column'
362
    fig.text(0.2, 0.72, txt, transform=fig.transFigure, size=14, ha="left")
363

    
364
    slist = ['Mean tropospheric ozone colum', 'Median tropospheric ozone column', 'Minimum tropospheric ozone column',
365
             'Maximum tropospheric ozone column', 'Standard deviation in tropospheric ozone columns',
366
             'Mean standard deviation of tropospheric ozone columns']
367
    vlist = [f_meanVCD, f_medianVCD, f_minVCD, f_maxVCD, f_stdVCD, f_meansigVCD]
368

    
369
    # print HTML table to stdout, as well as to pdf.
370
    print('<table class="table table-bordered">\n<tbody>')
371
    for i, (s, v) in enumerate(zip(slist, vlist)):
372
        print('<tr><td>{0}</td><td>{1:.2f} mmol/m²</td></tr>'.format(s, v))
373
        fig.text(0.2, 0.68-i*0.025, s, transform=fig.transFigure, size=14, ha="left")
374
        fig.text(0.52, 0.68-i*0.025, "{0:.2f} mmol/m²".format(v), transform=fig.transFigure, size=14, ha="left")
375
    print('</tbody>\n</table>')
376

    
377
    if out is None:
378
        plt.show()
379
    else:
380
        out.savefig()#,bbox_inches='tight')
381

    
382
    plt.close()
383

    
384
    # print(out, 'Statistics:')
385
    # print(out, 'tropospheric ozone column')
386
    # print(out, 'Mean tropospheric ozone column {0:04.2f} mmol/m²'.format(f_meanVCD))
387
    # print(out, 'Median tropospheric ozone column {0:04.2f} mmol/m²'.format(f_medianVCD))
388
    # print(out, 'Minimum tropospheric ozone column {0:04.2f} mmol/m²'.format(f_minVCD))
389
    # print(out, 'Maximum tropospheric ozone column {0:04.2f} mmol/m²'.format(f_maxVCD))
390
    # print(out, 'standard deviation in tropospheric ozone columns {0:04.2f} mmol/m²'.format(f_stdVCD))
391
    # print(out, 'mean standard deviation of tropospheric ozone columns {0:04.2f} mmol/m²'.format(f_meansigVCD))
392
    #
393
    #    print(out, 'upper tropospheric ozone mixing ratio')
394
    #    print(out, 'Mean upper tropospheric ozone mixing ratio {0:04.2f} ppb'.format(f_meanUMR))
395
    #    print(out, 'Median upper tropospheric ozone mixing ratio {0:04.2f} ppb'.format(f_medianUMR))
396
    #    print(out, 'Minimum upper tropospheric ozone mixing ratio {0:04.2f} ppb'.format(f_minUMR))
397
    #    print(out, 'Maximum upper tropospheric ozone mixing ratio {0:04.2f} ppb'.format(f_maxUMR))
398
    #    print(out, 'standard deviation in upper tropospheric ozone mixing ratios {0:04.2f} ppb'.format(f_stdUMR))
399
    #    print(out, 'mean standard deviation of upper tropospheric ozone mixing ratios {0:04.2f} ppb'.format(f_meansigUMR))
400

    
401
    return
402

    
403
def flags(results, out=None):
404
    import matplotlib.pyplot as plt
405
    from matplotlib.ticker import MaxNLocator
406
    from matplotlib.colors import BoundaryNorm
407
    import cartopy.crs as ccrs
408

    
409
    fa_VCD_strat = results._stratO3Reference
410
    ia_VCD_strat_flag = results._stratO3Reference_flag
411
    fa_latitude = results._latCenters
412

    
413
    fa_latC = results._latCentersCoarse
414
    fa_longC = results._lonCentersCoarse
415

    
416
    f_dlatC = abs(fa_latC[1]-fa_latC[0])
417
    f_dlongC = abs(fa_longC[1]-fa_longC[0])
418

    
419
    ia_UtropMR = results._upper_tropoO3_flag
420

    
421
    ## for the meshgrid not the centres but the corners should be used
422

    
423
    fa_latC = np.append(fa_latC-0.5*f_dlatC, fa_latC[-1]+0.5*f_dlatC)
424
    fa_longC = np.append(fa_longC-0.5*f_dlongC, fa_longC[-1]+0.5*f_dlongC)
425
    f_latmin = min(fa_latC)
426
    f_latmax = max(fa_latC)
427

    
428
    f_longmin = min(fa_longC)
429
    f_longmax = max(fa_longC)
430

    
431
    xc, yc = np.meshgrid(fa_longC, fa_latC)
432

    
433
    fig, (VCD, uMR) = plt.subplots(2, 1)
434
    fig.set_size_inches(20, 13)
435
    # plt.subplots_adjust(hspace=0.3)
436
    VCD = plt.subplot(2, 1, 1)
437
    VCD.plot(fa_latitude, fa_VCD_strat, 'b-')
438
    plt.axis([-22, 22, 0.1, 0.12])
439
    plt.ylabel('stratospheric reference column', color='b', fontsize=16)
440
    plt.yticks(fontsize=16)
441
    plt.xticks(fontsize=16)
442
    plt.xlabel('Latitude', fontsize=16)
443

    
444
    flag=VCD.twinx()
445
    flag.step(fa_latitude, ia_VCD_strat_flag, 'r')
446
    plt.axis([-22, 22, 0, 15])
447
    plt.ylabel('stratospheric reference column flag', color='r', fontsize=16)
448
    plt.yticks(fontsize=16)
449

    
450

    
451
    ## worldplot for upper tropospheric mixing ratio
452
    cm = plt.get_cmap('jet')       #gist_rainbow_r
453
    cm.set_under('white', 1)
454

    
455
    MR=plt.subplot(2, 1, 2, projection=ccrs.PlateCarree())
456
    MR.set_extent([-180, 180, -22, 22])
457

    
458
    # draw coastlines
459
    MR.coastlines()
460
    MR.gridlines(ylocs=(np.arange(-30, 30, 10)), xlocs=(-180, -90, 0, 90, 180))
461

    
462

    
463
    #    cs=MR.pcolormesh(x,y,ia_UtropMRnr,cmap=cm,norm=norm,transform=ccrs.PlateCarree())
464
    #
465
    levels = MaxNLocator(nbins=14).tick_values(1, 15)
466
    #    norm=BoundaryNorm(levels,ncolors=cm.N, clip=False)
467
    #    cs=uMR.pcolor(xc,yc,ia_UtropMR,cmap=cm,norm=norm,transform=ccrs.PlateCarree())
468
    #    cbar=plt.colorbar(cs,ax=uMR.ax,shrink=0.6,pad=0.11,aspect=35, orientation="horizontal",extendfrac='auto',extend='both')
469
    # cbar.set_label("",fontsize=16,labelpad=-5)
470
    #    cbar.set_ticks(np.arange(0,16,1))
471
    #    cbar.ax.tick_params(labelsize=16)
472
    s_title = 'Upper Tropospheric Ozone Mixing Ratio flag'
473
    plt.title(s_title,fontsize=16,horizontalalignment='left',x=0.05)
474

    
475
    if out is None:
476
        plt.show()
477
    else:
478
        out.savefig()#,bbox_inches='tight')
479
    plt.close()
480

    
481
    return()
482

    
483

    
484
##main()
485
def main():
486
    args = parseCmdLine()
487

    
488
    O3_TCLdata = readO3_TCL(args.input_file)
489

    
490
    #date as strings
491
    StartDate = str(O3_TCLdata._centralStart)
492
    EndDate   = str(O3_TCLdata._centralEnd)
493

    
494
    with PdfPages(args.output) as pdf:
495
        d = pdf.infodict()
496
        d['Title'] = 'PYCAMA report'
497
        d['Author'] = 'Klaus-Peter Heue'
498
        d['Subject'] = 'Tropospheric Ozone ({0}-{1})'.format(StartDate, EndDate)
499
        d['CreationDate'] = datetime.now()
500
        # print(pdf,d)
501
        statistics(O3_TCLdata, out=pdf)
502
        worldPlots(O3_TCLdata, out=pdf)
503
        Histogramplots(O3_TCLdata, out=pdf)
504
        flags(O3_TCLdata, out=pdf)
505

    
506
if __name__ == "__main__":
507
    main()