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/trunk/tests/STATION_ASF/EXPREF – NEMO

source: NEMO/trunk/tests/STATION_ASF/EXPREF/plot_station_asf_OCE.py @ 14072

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

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
  • 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
135def round_bounds( x1, x2,  base=5, prec=3 ):
136    rmin =  base * round( floor(float(x1)/base), prec )
137    rmax =  base * round(  ceil(float(x2)/base), prec )
138    return rmin, rmax
139
140
141# Getting time array from the first file:
142id_in = Dataset(cf_in[0])
143vt = id_in.variables['time_counter'][:]
144cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
145id_in.close()
146Nt = len(vt)
147
148vtime = num2date(vt, units=cunit_t) ; # something human!
149vtime = vtime.astype(dtype='datetime64[D]')
150
151ii=Nt/300
152ib=max(ii-ii%10,1)
153xticks_d=int(30*ib)
154
155rat = 100./float(rDPI)
156params = { 'font.family':'Open Sans',
157           'font.size':       int(15.*rat),
158           'legend.fontsize': int(15.*rat),
159           'xtick.labelsize': int(15.*rat),
160           'ytick.labelsize': int(15.*rat),
161           'axes.labelsize':  int(16.*rat)
162}
163mpl.rcParams.update(params)
164font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':18.*rat }
165font_x   = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':15.*rat }
166
167
168
169
170
171# First for each algorithm we compare some input vs out put variables:
172
173# t_skin
174vtemp_in = [ 'sst'           ,               'theta_zu'               ,               'theta_zt'               ,     't_skin'       ]
175vtemp_lb = [ 'SST (bulk SST)', r'$\theta_{zu}$ (pot. air temp. at zu)', r'$\theta_{zt}$ (pot. air temp. at zt)', 'Skin temperature' ]
176vtemp_cl = [  'k'            ,                clr_mod                 ,                'purple'                ,      clr_red       ]
177vtemp_lw = [   3             ,                   1.3                  ,                   1.3                  ,         0.7        ]
178vtemp_ls = [     '-'         ,                   '-'                  ,                   '--'                 ,         '-'        ]
179ntemp = len(vtemp_in)
180
181xxx = nmp.zeros((Nt,ntemp))
182
183for ja in range(nb_exp):
184    #
185    # Temperatures...
186    id_in = Dataset(cf_in[ja])
187    for jv in range(ntemp):
188        xxx[:,jv] = id_in.variables[vtemp_in[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
189    id_in.close()
190
191    fig = plt.figure(num = 1, figsize=size_fig0, facecolor='w', edgecolor='k')
192    ax1 = plt.axes([0.07, 0.2, 0.9, 0.75])
193    ax1.set_xticks(vtime[::xticks_d])
194    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
195    plt.xticks(rotation='60', **font_x)
196
197    for jv in range(ntemp):
198        plt.plot(vtime, xxx[:,jv], '-', color=vtemp_cl[jv], linestyle=vtemp_ls[jv], linewidth=vtemp_lw[jv], label=vtemp_lb[jv], zorder=10)
199
200    idx_okay = nmp.where( nmp.abs(xxx) < 1.e+10 )
201    fmin, fmax = round_bounds( nmp.min(xxx[idx_okay]) , nmp.max(xxx[idx_okay]), base=5, prec=0 )
202    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
203    plt.ylabel(r'Temperature [$^{\circ}$C]')
204
205    ax1.grid(color='k', linestyle='-', linewidth=0.3)
206    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
207    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)
208    plt.savefig('01_temperatures_'+list_exp[ja]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
209    plt.close(1)
210
211
212del xxx
213
214
215
216# Now we compare output variables from bulk algorithms between them:
217
218nb_var = len(L_VNEM)
219
220xF  = nmp.zeros((Nt,nb_exp))
221xFa = nmp.zeros((Nt,nb_exp))
222
223
224for jv in range(nb_var):
225    print('\n *** Treating variable: '+L_VARO[jv]+' !')
226
227    for ja in range(nb_exp):
228        #
229        id_in = Dataset(cf_in[ja])
230        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1]  ; # only the center point of the 3x3 spatial domain!
231        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
232        id_in.close()
233        #
234        id_toolarge, = nmp.where( xF[:,ja] > L_MAXT[jv] ) #
235        xF[id_toolarge,ja] = L_MAXT[jv]
236        id_toosmall, = nmp.where( xF[:,ja] < L_MINT[jv] ) ; #print("id_toosmall =", id_toosmall)
237        xF[id_toosmall,ja] = L_MINT[jv]
238
239    idx_okay = nmp.where( nmp.abs(xF) < 1.e+10 )
240
241    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
242    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
243    ax1 = plt.axes([0.083, 0.23, 0.9, 0.7])
244    ax1.set_xticks(vtime[::xticks_d])
245    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
246    plt.xticks(rotation='60', **font_x)
247
248    for ja in range(nb_exp):
249        fplot = nmp.ma.masked_where( xF[:,ja]==0., xF[:,ja] )
250        plt.plot(vtime, fplot, '-', color=l_color[ja], \
251                 linestyle=l_style[ja], linewidth=l_width[ja], label=list_exp[ja], alpha=0.6 ) #zorder=10+ja)
252
253    fmin, fmax = round_bounds( nmp.min(xF[idx_okay]) , nmp.max(xF[idx_okay]), base=L_BASE[jv], prec=L_PREC[jv])
254    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
255    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
256
257    ax1.grid(color='k', linestyle='-', linewidth=0.3)
258    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
259    ax1.annotate(cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
260                 ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
261                 zorder=50, **font_inf)
262    plt.savefig(L_VARO[jv]+'_'+cforcing+'_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
263    plt.close(jv)
264    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
265
266
267
268
269    def symetric_range( pmin, pmax ):
270        # Returns a symetric f-range that makes sense for the anomaly of "f" we're looking at...
271        from math import floor, copysign, log, ceil
272        zmax = max( abs(pmax) , abs(pmin) )
273        romagn = floor(log(zmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
274        rmlt = 10.**(int(romagn)) / 2.
275        frng = copysign( ceil(abs(zmax)/rmlt)*rmlt , zmax)
276        return frng
277
278
279    if L_ANOM[jv]:
280
281        for ja in range(nb_exp): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
282
283        if nmp.sum(nmp.abs(xFa[:,:])) == 0.0:
284            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
285            print('          ==> skipping anomaly plot...')
286
287        else:
288
289            yrng = symetric_range( nmp.min(xFa) , nmp.max(xFa) )
290
291            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
292            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
293            ax1 = plt.axes([0.09, 0.23, 0.9, 0.7])
294            ax1.set_xticks(vtime[::xticks_d])
295            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
296            plt.xticks(rotation='60', **font_x)
297
298            for ja in range(nb_exp):
299                fplot = nmp.ma.masked_where( xFa[:,ja]==0., xFa[:,ja] )
300                plt.plot(vtime, fplot, '-', color=l_color[ja], \
301                         linewidth=l_width[ja], label=list_exp[ja], alpha=0.6) #, zorder=10+ja)
302
303            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
304            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
305            ax1.grid(color='k', linestyle='-', linewidth=0.3)
306            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
307            ax1.annotate('Anomaly of '+cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
308                         ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
309                         zorder=50, **font_inf)
310            plt.savefig(L_VARO[jv]+'_'+cforcing+'_anomaly_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
311            plt.close(10+jv)
312            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Note: See TracBrowser for help on using the repository browser.