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 @ 13723

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

Keep up with trunk r13718 + figure generation for STATION_ASF (ocean & ice).
  • Property svn:executable set to *
File size: 12.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 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
30l_color = [ '0.85'    ,  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot
31l_width = [  4        ,     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot
32l_style = [ '-'       ,    '-'     ,   '-'     , '--'   ] ; # line-style
33
34
35# Variables to compare between algorithms
36############################################
37crealm = 'open-ocean'
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.0005   ,    0.0005   ,      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# About STATION_ASF output files to read:
47cpref = 'STATION_ASF-'          ; np = len(cpref)
48csuff = '_gridT.nc'             ; ns = len(csuff)
49cclnd = '_1h_YYYY0101_YYYY1231' ; nc = len(cclnd)
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# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
89
90dir_out = cdir_data+'/output'
91ldir = listdir(dir_out)
92
93cf_in    = []
94list_exp = []
95list_frc = []
96for fn in ldir:
97    fpn = dir_out+'/'+fn
98    if path.isfile(dir_out+'/'+fn):
99        if fn[:np]==cpref and fn[-ns:]==csuff and cforcing in fn:
100            print('\n file: '+fn)
101            clab = fn[np:-nc-ns]
102            [ cexp, cfrc ] = str.split(clab, '_', 1)
103            print('  ===> Experiment = '+cexp+', Forcing = '+cfrc)
104            list_exp.append(cexp)
105            list_frc.append(cfrc)
106            cf_in.append(fpn)
107nbf = len( set(list_frc) )
108if not nbf == 1:
109    print('PROBLEM: we found files for more that one forcing: ', set(list_frc))
110    sys.exit(0)
111
112nb_exp = len(list_exp)
113
114
115print('\n\nThere are '+str(nb_exp)+' experiments to compare:')
116for ja in range(nb_exp): print('  * '+list_exp[ja]+'\n'+'     ==> '+cf_in[ja]+'\n')
117
118if nb_exp > len(l_color):
119    print('PROBLEM: the max number of experiments for comparison is '+str(len(l_color))+' for now...')
120    sys.exit(0)
121
122
123
124#-----------------------------------------------------------------
125
126
127def round_bounds( x1, x2,  base=5, prec=3 ):
128    rmin =  base * round( floor(float(x1)/base), prec )
129    rmax =  base * round(  ceil(float(x2)/base), prec )
130    return rmin, rmax
131
132
133# Getting time array from the first file:
134id_in = Dataset(cf_in[0])
135vt = id_in.variables['time_counter'][:]
136cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
137id_in.close()
138Nt = len(vt)
139
140vtime = num2date(vt, units=cunit_t) ; # something human!
141vtime = vtime.astype(dtype='datetime64[D]')
142
143ii=Nt/300
144ib=max(ii-ii%10,1)
145xticks_d=int(30*ib)
146
147rat = 100./float(rDPI)
148params = { 'font.family':'Open Sans',
149           'font.size':       int(15.*rat),
150           'legend.fontsize': int(15.*rat),
151           'xtick.labelsize': int(15.*rat),
152           'ytick.labelsize': int(15.*rat),
153           'axes.labelsize':  int(16.*rat)
154}
155mpl.rcParams.update(params)
156font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':18.*rat }
157font_x   = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':15.*rat }
158
159
160
161
162
163# First for each algorithm we compare some input vs out put variables:
164
165# t_skin
166vtemp_in = [ 'sst'           ,               'theta_zu'               ,               'theta_zt'               ,     't_skin'       ]
167vtemp_lb = [ 'SST (bulk SST)', r'$\theta_{zu}$ (pot. air temp. at zu)', r'$\theta_{zt}$ (pot. air temp. at zt)', 'Skin temperature' ]
168vtemp_cl = [  'k'            ,                clr_mod                 ,                'purple'                ,      clr_red       ]
169vtemp_lw = [   3             ,                   1.3                  ,                   1.3                  ,         0.7        ]
170vtemp_ls = [     '-'         ,                   '-'                  ,                   '--'                 ,         '-'        ]
171ntemp = len(vtemp_in)
172
173xxx = nmp.zeros((Nt,ntemp))
174
175for ja in range(nb_exp):
176    #
177    # Temperatures...
178    id_in = Dataset(cf_in[ja])
179    for jv in range(ntemp):
180        xxx[:,jv] = id_in.variables[vtemp_in[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
181    id_in.close()
182
183    fig = plt.figure(num = 1, figsize=size_fig0, facecolor='w', edgecolor='k')
184    ax1 = plt.axes([0.07, 0.2, 0.9, 0.75])
185    ax1.set_xticks(vtime[::xticks_d])
186    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
187    plt.xticks(rotation='60', **font_x)
188
189    for jv in range(ntemp):
190        plt.plot(vtime, xxx[:,jv], '-', color=vtemp_cl[jv], linestyle=vtemp_ls[jv], linewidth=vtemp_lw[jv], label=vtemp_lb[jv], zorder=10)
191
192    idx_okay = nmp.where( nmp.abs(xxx) < 1.e+10 )
193    fmin, fmax = round_bounds( nmp.min(xxx[idx_okay]) , nmp.max(xxx[idx_okay]), base=5, prec=0 )
194    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
195    plt.ylabel(r'Temperature [$^{\circ}$C]')
196
197    ax1.grid(color='k', linestyle='-', linewidth=0.3)
198    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
199    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)
200    plt.savefig('01_temperatures_'+list_exp[ja]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
201    plt.close(1)
202
203
204del xxx
205
206
207
208# Now we compare output variables from bulk algorithms between them:
209
210nb_var = len(L_VNEM)
211
212xF  = nmp.zeros((Nt,nb_exp))
213xFa = nmp.zeros((Nt,nb_exp))
214
215
216for jv in range(nb_var):
217    print('\n *** Treating variable: '+L_VARO[jv]+' !')
218
219    for ja in range(nb_exp):
220        #
221        id_in = Dataset(cf_in[ja])
222        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
223        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
224        id_in.close()
225
226    idx_okay = nmp.where( nmp.abs(xF) < 1.e+10 )
227
228    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
229    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
230    ax1 = plt.axes([0.083, 0.23, 0.9, 0.7])
231    ax1.set_xticks(vtime[::xticks_d])
232    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
233    plt.xticks(rotation='60', **font_x)
234
235    for ja in range(nb_exp):
236        fplot = nmp.ma.masked_where( xF[:,ja]==0., xF[:,ja] )
237        plt.plot(vtime, fplot, '-', color=l_color[ja], \
238                 linestyle=l_style[ja], linewidth=l_width[ja], label=list_exp[ja], zorder=10+ja)
239
240    fmin, fmax = round_bounds( nmp.min(xF[idx_okay]) , nmp.max(xF[idx_okay]), base=L_BASE[jv], prec=L_PREC[jv])
241    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
242    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
243
244    ax1.grid(color='k', linestyle='-', linewidth=0.3)
245    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
246    ax1.annotate(cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
247                 ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
248                 zorder=50, **font_inf)
249    plt.savefig(L_VARO[jv]+'_'+cforcing+'_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
250    plt.close(jv)
251    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
252
253
254
255
256    def symetric_range( pmin, pmax ):
257        # Returns a symetric f-range that makes sense for the anomaly of "f" we're looking at...
258        from math import floor, copysign, log, ceil
259        zmax = max( abs(pmax) , abs(pmin) )
260        romagn = floor(log(zmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
261        rmlt = 10.**(int(romagn)) / 2.
262        frng = copysign( ceil(abs(zmax)/rmlt)*rmlt , zmax)
263        return frng
264
265
266    if L_ANOM[jv]:
267
268        for ja in range(nb_exp): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
269
270        if nmp.sum(nmp.abs(xFa[:,:])) == 0.0:
271            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
272            print('          ==> skipping anomaly plot...')
273
274        else:
275
276            yrng = symetric_range( nmp.min(xFa) , nmp.max(xFa) )
277
278            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
279            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
280            ax1 = plt.axes([0.09, 0.23, 0.9, 0.7])
281            ax1.set_xticks(vtime[::xticks_d])
282            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
283            plt.xticks(rotation='60', **font_x)
284
285            for ja in range(nb_exp):
286                fplot = nmp.ma.masked_where( xF[:,ja]==0., xF[:,ja] )
287                plt.plot(vtime, fplot, '-', color=l_color[ja], \
288                         linewidth=l_width[ja], label=list_exp[ja], zorder=10+ja)
289
290            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
291            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
292            ax1.grid(color='k', linestyle='-', linewidth=0.3)
293            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
294            ax1.annotate('Anomaly of '+cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
295                         ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
296                         zorder=50, **font_inf)
297            plt.savefig(L_VARO[jv]+'_'+cforcing+'_anomaly_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
298            plt.close(10+jv)
299            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Note: See TracBrowser for help on using the repository browser.