source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/STATION_ASF/EXPREF/plot_station_asf_OCE.py @ 13690

Last change on this file since 13690 was 13690, checked in by laurent, 6 months ago

Improvements for "oce+ice capable" STATION_ASF!

  • Property svn:executable set to *
File size: 11.5 KB
Line 
1#!/usr/bin/env python3
2# -*- Mode: Python; coding: utf-8; indent-tabs-mode: nil; tab-width: 4 -*-
3#
4##########################################################################################
5# Post-diagnostic of STATION_ASF for ocean only  (no sea-ice support)
6#
7#  L. Brodeau, 2020
8##########################################################################################
9
10import sys
11from os import path as path
12import argparse as ap
13from math import floor, ceil, copysign, log
14import numpy as nmp
15from netCDF4 import Dataset,num2date
16import matplotlib as mpl
17mpl.use('Agg')
18import matplotlib.pyplot as plt
19import matplotlib.dates as mdates
20
21dir_figs='.'
22size_fig=(13,8)
23size_fig0=(12,10)
24
25clr_red = '#AD0000'
26clr_mod = '#008ab8'
27
28rDPI=100.
29
30L_ALGOS = [ 'ANDREAS' , 'COARE3p6' , 'ECMWF'   , 'NCAR' ]
31l_color = [ '0.85'    ,  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot
32l_width = [  4        ,     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot
33l_style = [ '-'       ,    '-'     ,   '-'     , '--'   ] ; # line-style
34
35
36# Variables to compare between algorithms
37############################################
38L_VNEM  = [   'Cd_oce'  ,   'Ce_oce'  ,   'qla_oce' , 'qsb_oce'     ,     'qt_oce' ,   'qlw_oce' ,  'taum'     ,    'dt_skin'         ]
39L_VARO  = [     'Cd'    ,     'Ce'    ,   'Qlat'    ,    'Qsen'     ,     'Qnet'   ,   'Qlw'     ,  'Tau'      ,    'dT_skin'         ] ; # name of variable on figure
40L_VARL  = [ r'$C_{D}$'  , r'$C_{E}$'  , r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$|\tau|$' , r'$\Delta T_{skin}$' ] ; # name of variable in latex mode
41L_VUNT  = [     ''      ,     ''      , r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$N/m^2$'  ,      'K'             ]
42L_BASE  = [    0.005    ,    0.005    ,      5.     ,      5.       ,       5      ,      5.     ,    0.05      ,       0.05           ]
43L_PREC  = [      3      ,      3      ,      0      ,      0        ,       0      ,      0      ,     2       ,        3             ]
44L_ANOM  = [   False     ,   False     ,   True      ,    True       ,    True      ,    True     ,   True      ,      False           ]
45
46
47nb_algos = len(L_ALGOS) ; print(nb_algos)
48
49
50
51
52################## ARGUMENT PARSING / USAGE ################################################################################################
53parser = ap.ArgumentParser(description='Generate pixel maps of a given scalar.')
54#
55requiredNamed = parser.add_argument_group('required arguments')
56requiredNamed.add_argument('-d', '--dirout' , required=True,                 help='Path to (production) directory where STATION_ASF was run')
57requiredNamed.add_argument('-f', '--forcing', required=True, default="PAPA", help='Name of forcing (ex: PAPA, ERA5_arctic')
58#
59parser.add_argument('-C', '--conf',   default="STATION_ASF",  help='specify NEMO config (ex: STATION_ASF)')
60parser.add_argument('-s', '--ystart', default="2018",         help='specify first year of experiment (ex: 2018)')
61parser.add_argument('-e', '--yend',   default="2018",         help='specify last  year of experiment (ex: 2018)')
62#
63parser.add_argument('-t', '--itype',   default="png",         help='specify the type of image you want to create (ex: png, svg, etc.)')
64#parser.add_argument('-l', '--lev' , type=int, default=0,    help='specify the level to use if 3D field (default: 0 => 2D)')
65#parser.add_argument('-I', '--ice' , action='store_true',    help='draw sea-ice concentration layer onto the field')
66#
67args = parser.parse_args()
68#
69cdir_data = args.dirout
70cforcing  = args.forcing
71#
72CONF      = args.conf
73cy1       = args.ystart
74cy2       = args.yend
75#
76fig_ext   = args.itype
77#jk    = args.lev
78#lshow_ice = args.ice
79#
80#print(''); print(' *** cdir_data = ', cdir_data); print(' *** cforcing  = ', cforcing)
81#print(' *** CONF = ', CONF); print(' *** cy1 = ', cy1); print(' *** cy2 = ', cy2)
82###############################################################################################################################################
83
84
85
86# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
87# Populating and checking existence of files to be read
88# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
89def chck4f(cf):
90    cmesg = 'ERROR: File '+cf+' does not exist !!!'
91    if not path.exists(cf): print(cmesg) ; sys.exit(0)
92
93###cf_in = nmp.empty((), dtype="S10")
94cf_in = []
95for ja in range(nb_algos):
96    cfi = cdir_data+'/output/'+CONF+'-'+L_ALGOS[ja]+'_'+cforcing+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc'
97    chck4f(cfi)
98    cf_in.append(cfi)
99print('Files we are goin to use:')
100for ja in range(nb_algos): print(cf_in[ja])
101#-----------------------------------------------------------------
102
103
104def round_bounds( x1, x2,  base=5, prec=3 ):
105    rmin =  base * round( floor(float(x1)/base), prec )
106    rmax =  base * round(  ceil(float(x2)/base), prec )
107    return rmin, rmax
108
109
110
111# Getting time array from the first file:
112id_in = Dataset(cf_in[0])
113vt = id_in.variables['time_counter'][:]
114cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
115id_in.close()
116nbr = len(vt)
117
118vtime = num2date(vt, units=cunit_t) ; # something understandable!                                                                 
119vtime = vtime.astype(dtype='datetime64[D]')
120
121ii=nbr/300
122ib=max(ii-ii%10,1)
123xticks_d=int(30*ib)
124
125
126rat = 100./float(rDPI)
127params = { 'font.family':'Open Sans',
128           'font.size':       int(15.*rat),
129           'legend.fontsize': int(15.*rat),
130           'xtick.labelsize': int(15.*rat),
131           'ytick.labelsize': int(15.*rat),
132           'axes.labelsize':  int(16.*rat)
133}
134mpl.rcParams.update(params)
135font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':18.*rat }
136font_x   = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':15.*rat }
137
138
139
140
141
142# First for each algorithm we compare some input vs out put variables:                                                                       
143
144# t_skin
145vtemp_in = [ 'sst'            ,               'theta_zu'                 ,               't_skin'                   ] ; #,               'theta_zt'               
146vtemp_lb = [ 'SST (bulk SST)' , r'$\theta_{zu}$ (pot. air temp. at 10m)' ,           'Skin temperature'             ] ; #, r'$\theta_{zt}$ (pot. air temp. at 2m)'
147vtemp_cl = [  'k'             ,                clr_mod                   ,                clr_red                   ] ; #,                clr_mod                 
148vtemp_lw = [   3              ,                   1.3                    ,                   0.7                    ] ; #,                 2                       
149vtemp_ls = [     '-'          ,                   '-'                    ,                   '-'                    ] ; #,                 '-'                     
150
151ntemp = len(vtemp_in)
152
153xxx = nmp.zeros((nbr,ntemp))
154
155for ja in range(nb_algos):
156    #
157    # Temperatures...
158    id_in = Dataset(cf_in[ja])
159    for jv in range(ntemp):
160        xxx[:,jv] = id_in.variables[vtemp_in[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
161    id_in.close()
162
163    fig = plt.figure(num = 1, figsize=size_fig0, facecolor='w', edgecolor='k')
164    ax1 = plt.axes([0.07, 0.2, 0.9, 0.75])
165    ax1.set_xticks(vtime[::xticks_d])
166    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
167    plt.xticks(rotation='60', **font_x)
168
169    for jv in range(ntemp):
170        plt.plot(vtime, xxx[:,jv], '-', color=vtemp_cl[jv], linestyle=vtemp_ls[jv], linewidth=vtemp_lw[jv], label=vtemp_lb[jv], zorder=10)
171
172    idx_okay = nmp.where( nmp.abs(xxx) < 1.e+10 )
173    fmin, fmax = round_bounds( nmp.min(xxx[idx_okay]) , nmp.max(xxx[idx_okay]), base=5, prec=0 )
174    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
175    plt.ylabel(r'Temperature [$^{\circ}$C]')
176
177    ax1.grid(color='k', linestyle='-', linewidth=0.3)
178    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
179    ax1.annotate('Algo: '+L_ALGOS[ja]+', station: '+cforcing, xy=(0.4, 1.), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
180    plt.savefig('01_temperatures_'+L_ALGOS[ja]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
181    plt.close(1)
182
183
184del xxx
185
186
187
188# Now we compare output variables from bulk algorithms between them:
189
190nb_var = len(L_VNEM)
191
192xF  = nmp.zeros((nbr,nb_algos))
193xFa = nmp.zeros((nbr,nb_algos))
194
195
196for jv in range(nb_var):
197    print('\n *** Treating variable: '+L_VARO[jv]+' !')
198
199    for ja in range(nb_algos):
200        #
201        id_in = Dataset(cf_in[ja])
202        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
203        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
204        id_in.close()
205
206    idx_okay = nmp.where( nmp.abs(xF) < 1.e+10 )
207       
208    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
209    ax1 = plt.axes([0.08, 0.25, 0.9, 0.7])
210    ax1.set_xticks(vtime[::xticks_d])
211    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
212    plt.xticks(rotation='60', **font_x)
213
214    for ja in range(nb_algos):
215        plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
216       
217    fmin, fmax = round_bounds( nmp.min(xF[idx_okay]) , nmp.max(xF[idx_okay]), base=L_BASE[jv], prec=L_PREC[jv])
218    #print("LOLO: fmin, fmax =",nmp.min(xF[idx_okay]), nmp.max(xF[idx_okay]) ); print("LOLO: fmin, fmax =",fmin, fmax) ; #sys.exit(0)
219   
220    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
221    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
222
223    ax1.grid(color='k', linestyle='-', linewidth=0.3)
224    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
225    ax1.annotate(cvar_lnm+', station: '+cforcing, xy=(0.3, 1.), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
226                 zorder=50, **font_inf)
227    plt.savefig(L_VARO[jv]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
228    plt.close(jv)
229
230    if L_ANOM[jv]:
231
232        for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
233
234        if nmp.sum(xFa[:,:]) == 0.0:
235            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
236            print('          ==> skipping anomaly plot...')
237
238        else:
239
240            # Want a symetric y-range that makes sense for the anomaly we're looking at:
241            rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
242            rmax = max( abs(rmax) , abs(rmin) )
243            romagn = floor(log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
244            rmlt = 10.**(int(romagn)) / 2.
245            yrng = copysign( ceil(abs(rmax)/rmlt)*rmlt , rmax)
246
247            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
248            ax1 = plt.axes([0.08, 0.25, 0.9, 0.7])
249            ax1.set_xticks(vtime[::xticks_d])
250            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
251            plt.xticks(rotation='60', **font_x)
252
253            for ja in range(nb_algos):
254                plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
255
256            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
257            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
258            ax1.grid(color='k', linestyle='-', linewidth=0.3)
259            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
260            ax1.annotate('Anomaly of '+cvar_lnm, xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
261            plt.savefig(L_VARO[jv]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
262            plt.close(10+jv)
263
264
265
266
Note: See TracBrowser for help on using the repository browser.