New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
plot_station_asf_OCE.py in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/STATION_ASF/EXPREF – NEMO

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

Last change on this file since 13893 was 13893, checked in by laurent, 3 years ago

Improvements and better README for tests/STATION_ASF

  • Property svn:executable set to *
File size: 13.3 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 air-sea fluxes (over open ocean)
6#
7#  L. Brodeau, 2020
8##########################################################################################
9
10import sys
11from os import path, listdir
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.5)
23size_fig0=(12,10)
24
25clr_red = '#AD0000'
26clr_mod = '#008ab8'
27
28rDPI=100.
29
30#l_color = [ '0.85'    ,  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot
31#l_width = [  4        ,     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot
32#l_style = [ '-'       ,    '-'     ,   '-'     , '--'   ] ; # line-style
33
34#ffed00: yellow ON
35#E8A727: ornage
36
37l_color = [ '0.3' , '#E8A727', '0.1'  , '#008ab8'  ] ; # colors to differentiate algos on the plot
38l_width = [   2   ,    2     ,  1.5   ,    2       ] ; # line-width to differentiate algos on the plot
39l_style = [  '-'  ,    '-'   , '--'   ,   '-'      ] ; # line-style
40
41
42# Variables to compare between algorithms
43############################################
44crealm = 'open-ocean'
45L_VNEM = [   'Cd_oce'  ,   'Ce_oce'  ,   'qla_oce' , 'qsb_oce'     ,     'qt_oce' ,   'qlw_oce' ,  'taum'     ,    'dt_skin'         ]
46L_VARO = [     'Cd'    ,     'Ce'    ,   'Qlat'    ,    'Qsen'     ,     'Qnet'   ,   'Qlw'     ,  'Tau'      ,    'dT_skin'         ] ; # name of variable on figure
47L_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
48L_VUNT = [     ''      ,     ''      , r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$N/m^2$'  ,      'K'             ]
49L_BASE = [    0.0005   ,    0.0005   ,      5.     ,      5.       ,       5      ,      5.     ,    0.05      ,       0.05           ]
50L_PREC = [      3      ,      3      ,      0      ,      0        ,       0      ,      0      ,     2       ,        3             ]
51L_ANOM = [   False     ,   False     ,   True      ,    True       ,    True      ,    True     ,   True      ,      False           ]
52L_MAXT = [  10000.   ,     10000.,  10000.   ,  10000.      ,  10000.    ,  10000.      , 10000.     ,      1.5    ]
53L_MINT = [    0.001  ,     0.001 , -10000.   , -10000.      , -10000.    , -10000.      ,-10000.     ,   -10000.   ]
54
55# About STATION_ASF output files to read:
56cpref = 'STATION_ASF-'          ; np = len(cpref)
57csuff = '_gridT.nc'             ; ns = len(csuff)
58cclnd = '_1h_YYYY0101_YYYY1231' ; nc = len(cclnd)
59
60
61################## ARGUMENT PARSING / USAGE ################################################################################################
62parser = ap.ArgumentParser(description='Generate pixel maps of a given scalar.')
63#
64requiredNamed = parser.add_argument_group('required arguments')
65requiredNamed.add_argument('-d', '--dirout' , required=True,                 help='Path to (production) directory where STATION_ASF was run')
66requiredNamed.add_argument('-f', '--forcing', required=True, default="PAPA", help='Name of forcing (ex: PAPA, ERA5_arctic')
67#
68parser.add_argument('-C', '--conf',   default="STATION_ASF",  help='specify NEMO config (ex: STATION_ASF)')
69parser.add_argument('-s', '--ystart', default="2018",         help='specify first year of experiment (ex: 2018)')
70parser.add_argument('-e', '--yend',   default="2018",         help='specify last  year of experiment (ex: 2018)')
71#
72parser.add_argument('-t', '--itype',   default="png",         help='specify the type of image you want to create (ex: png, svg, etc.)')
73#parser.add_argument('-l', '--lev' , type=int, default=0,    help='specify the level to use if 3D field (default: 0 => 2D)')
74#parser.add_argument('-I', '--ice' , action='store_true',    help='draw sea-ice concentration layer onto the field')
75#
76args = parser.parse_args()
77#
78cdir_data = args.dirout
79cforcing  = args.forcing
80#
81CONF      = args.conf
82cy1       = args.ystart
83cy2       = args.yend
84#
85fig_ext   = args.itype
86#jk    = args.lev
87#lshow_ice = args.ice
88#
89#print(''); print(' *** cdir_data = ', cdir_data); print(' *** cforcing  = ', cforcing)
90#print(' *** CONF = ', CONF); print(' *** cy1 = ', cy1); print(' *** cy2 = ', cy2)
91###############################################################################################################################################
92
93
94
95# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
96# Populating and checking existence of files to be read
97# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
98
99dir_out = cdir_data+'/output'
100ldir = listdir(dir_out)
101
102cf_in    = []
103list_exp = []
104list_frc = []
105for fn in ldir:
106    fpn = dir_out+'/'+fn
107    if path.isfile(dir_out+'/'+fn):
108        if fn[:np]==cpref and fn[-ns:]==csuff and cforcing in fn:
109            print('\n file: '+fn)
110            clab = fn[np:-nc-ns]
111            [ cexp, cfrc ] = str.split(clab, '_', 1)
112            print('  ===> Experiment = '+cexp+', Forcing = '+cfrc)
113            list_exp.append(cexp)
114            list_frc.append(cfrc)
115            cf_in.append(fpn)
116nbf = len( set(list_frc) )
117if not nbf == 1:
118    print('PROBLEM: we found files for more that one forcing: ', set(list_frc))
119    sys.exit(0)
120
121nb_exp = len(list_exp)
122
123
124print('\n\nThere are '+str(nb_exp)+' experiments to compare:')
125for ja in range(nb_exp): print('  * '+list_exp[ja]+'\n'+'     ==> '+cf_in[ja]+'\n')
126
127if nb_exp > len(l_color):
128    print('PROBLEM: the max number of experiments for comparison is '+str(len(l_color))+' for now...')
129    sys.exit(0)
130
131
132
133#-----------------------------------------------------------------
134
135
136def round_bounds( x1, x2,  base=5, prec=3 ):
137    rmin =  base * round( floor(float(x1)/base), prec )
138    rmax =  base * round(  ceil(float(x2)/base), prec )
139    return rmin, rmax
140
141
142# Getting time array from the first file:
143id_in = Dataset(cf_in[0])
144vt = id_in.variables['time_counter'][:]
145cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
146id_in.close()
147Nt = len(vt)
148
149vtime = num2date(vt, units=cunit_t) ; # something human!
150vtime = vtime.astype(dtype='datetime64[D]')
151
152ii=Nt/300
153ib=max(ii-ii%10,1)
154xticks_d=int(30*ib)
155
156rat = 100./float(rDPI)
157params = { 'font.family':'Open Sans',
158           'font.size':       int(15.*rat),
159           'legend.fontsize': int(15.*rat),
160           'xtick.labelsize': int(15.*rat),
161           'ytick.labelsize': int(15.*rat),
162           'axes.labelsize':  int(16.*rat)
163}
164mpl.rcParams.update(params)
165font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':18.*rat }
166font_x   = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':15.*rat }
167
168
169
170
171
172# First for each algorithm we compare some input vs out put variables:
173
174# t_skin
175vtemp_in = [ 'sst'           ,               'theta_zu'               ,               'theta_zt'               ,     't_skin'       ]
176vtemp_lb = [ 'SST (bulk SST)', r'$\theta_{zu}$ (pot. air temp. at zu)', r'$\theta_{zt}$ (pot. air temp. at zt)', 'Skin temperature' ]
177vtemp_cl = [  'k'            ,                clr_mod                 ,                'purple'                ,      clr_red       ]
178vtemp_lw = [   3             ,                   1.3                  ,                   1.3                  ,         0.7        ]
179vtemp_ls = [     '-'         ,                   '-'                  ,                   '--'                 ,         '-'        ]
180ntemp = len(vtemp_in)
181
182xxx = nmp.zeros((Nt,ntemp))
183
184for ja in range(nb_exp):
185    #
186    # Temperatures...
187    id_in = Dataset(cf_in[ja])
188    for jv in range(ntemp):
189        xxx[:,jv] = id_in.variables[vtemp_in[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
190    id_in.close()
191
192    fig = plt.figure(num = 1, figsize=size_fig0, facecolor='w', edgecolor='k')
193    ax1 = plt.axes([0.07, 0.2, 0.9, 0.75])
194    ax1.set_xticks(vtime[::xticks_d])
195    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
196    plt.xticks(rotation='60', **font_x)
197
198    for jv in range(ntemp):
199        plt.plot(vtime, xxx[:,jv], '-', color=vtemp_cl[jv], linestyle=vtemp_ls[jv], linewidth=vtemp_lw[jv], label=vtemp_lb[jv], zorder=10)
200
201    idx_okay = nmp.where( nmp.abs(xxx) < 1.e+10 )
202    fmin, fmax = round_bounds( nmp.min(xxx[idx_okay]) , nmp.max(xxx[idx_okay]), base=5, prec=0 )
203    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
204    plt.ylabel(r'Temperature [$^{\circ}$C]')
205
206    ax1.grid(color='k', linestyle='-', linewidth=0.3)
207    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
208    ax1.annotate('Algo: '+list_exp[ja]+', station: '+cforcing, xy=(0.5, 1.), xycoords='axes fraction', ha='center',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
209    plt.savefig('01_temperatures_'+list_exp[ja]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
210    plt.close(1)
211
212
213del xxx
214
215
216
217# Now we compare output variables from bulk algorithms between them:
218
219nb_var = len(L_VNEM)
220
221xF  = nmp.zeros((Nt,nb_exp))
222xFa = nmp.zeros((Nt,nb_exp))
223
224
225for jv in range(nb_var):
226    print('\n *** Treating variable: '+L_VARO[jv]+' !')
227
228    for ja in range(nb_exp):
229        #
230        id_in = Dataset(cf_in[ja])
231        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1]  ; # only the center point of the 3x3 spatial domain!
232        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
233        id_in.close()
234        #
235        id_toolarge, = nmp.where( xF[:,ja] > L_MAXT[jv] ) #
236        xF[id_toolarge,ja] = L_MAXT[jv]
237        id_toosmall, = nmp.where( xF[:,ja] < L_MINT[jv] ) ; print("id_toosmall =", id_toosmall)
238        xF[id_toosmall,ja] = L_MINT[jv]
239
240    idx_okay = nmp.where( nmp.abs(xF) < 1.e+10 )
241
242    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
243    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
244    ax1 = plt.axes([0.083, 0.23, 0.9, 0.7])
245    ax1.set_xticks(vtime[::xticks_d])
246    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
247    plt.xticks(rotation='60', **font_x)
248
249    for ja in range(nb_exp):
250        fplot = nmp.ma.masked_where( xF[:,ja]==0., xF[:,ja] )
251        plt.plot(vtime, fplot, '-', color=l_color[ja], \
252                 linestyle=l_style[ja], linewidth=l_width[ja], label=list_exp[ja], alpha=0.6 ) #zorder=10+ja)
253
254    fmin, fmax = round_bounds( nmp.min(xF[idx_okay]) , nmp.max(xF[idx_okay]), base=L_BASE[jv], prec=L_PREC[jv])
255    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
256    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
257
258    ax1.grid(color='k', linestyle='-', linewidth=0.3)
259    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
260    ax1.annotate(cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
261                 ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
262                 zorder=50, **font_inf)
263    plt.savefig(L_VARO[jv]+'_'+cforcing+'_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
264    plt.close(jv)
265    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
266
267
268
269
270    def symetric_range( pmin, pmax ):
271        # Returns a symetric f-range that makes sense for the anomaly of "f" we're looking at...
272        from math import floor, copysign, log, ceil
273        zmax = max( abs(pmax) , abs(pmin) )
274        romagn = floor(log(zmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
275        rmlt = 10.**(int(romagn)) / 2.
276        frng = copysign( ceil(abs(zmax)/rmlt)*rmlt , zmax)
277        return frng
278
279
280    if L_ANOM[jv]:
281
282        for ja in range(nb_exp): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
283
284        if nmp.sum(nmp.abs(xFa[:,:])) == 0.0:
285            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
286            print('          ==> skipping anomaly plot...')
287
288        else:
289
290            yrng = symetric_range( nmp.min(xFa) , nmp.max(xFa) )
291
292            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
293            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
294            ax1 = plt.axes([0.09, 0.23, 0.9, 0.7])
295            ax1.set_xticks(vtime[::xticks_d])
296            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
297            plt.xticks(rotation='60', **font_x)
298
299            for ja in range(nb_exp):
300                fplot = nmp.ma.masked_where( xF[:,ja]==0., xF[:,ja] )
301                plt.plot(vtime, fplot, '-', color=l_color[ja], \
302                         linewidth=l_width[ja], label=list_exp[ja], alpha=0.6) #, zorder=10+ja)
303
304            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
305            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
306            ax1.grid(color='k', linestyle='-', linewidth=0.3)
307            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
308            ax1.annotate('Anomaly of '+cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
309                         ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
310                         zorder=50, **font_inf)
311            plt.savefig(L_VARO[jv]+'_'+cforcing+'_anomaly_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
312            plt.close(10+jv)
313            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Note: See TracBrowser for help on using the repository browser.