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

Last change on this file since 13691 was 13691, checked in by laurent, 4 years ago

Improvements for "oce+ice capable" STATION_ASF!
  • Property svn:executable set to *
File size: 11.4 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.5)
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'               ,               'theta_zt'               ,     't_skin'       ]
146vtemp_lb = [ 'SST (bulk SST)', r'$\theta_{zu}$ (pot. air temp. at zu)', r'$\theta_{zt}$ (pot. air temp. at zt)', 'Skin temperature' ]
147vtemp_cl = [  'k'            ,                clr_mod                 ,                'purple'                ,      clr_red       ]
148vtemp_lw = [   3             ,                   1.3                  ,                   1.3                  ,         0.7        ]
149vtemp_ls = [     '-'         ,                   '-'                  ,                   '--'                 ,         '-'        ]
150ntemp = len(vtemp_in)
151
152xxx = nmp.zeros((nbr,ntemp))
153
154for ja in range(nb_algos):
155    #
156    # Temperatures...
157    id_in = Dataset(cf_in[ja])
158    for jv in range(ntemp):
159        xxx[:,jv] = id_in.variables[vtemp_in[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
160    id_in.close()
161
162    fig = plt.figure(num = 1, figsize=size_fig0, facecolor='w', edgecolor='k')
163    ax1 = plt.axes([0.07, 0.2, 0.9, 0.75])
164    ax1.set_xticks(vtime[::xticks_d])
165    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
166    plt.xticks(rotation='60', **font_x)
167
168    for jv in range(ntemp):
169        plt.plot(vtime, xxx[:,jv], '-', color=vtemp_cl[jv], linestyle=vtemp_ls[jv], linewidth=vtemp_lw[jv], label=vtemp_lb[jv], zorder=10)
170
171    idx_okay = nmp.where( nmp.abs(xxx) < 1.e+10 )
172    fmin, fmax = round_bounds( nmp.min(xxx[idx_okay]) , nmp.max(xxx[idx_okay]), base=5, prec=0 )
173    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
174    plt.ylabel(r'Temperature [$^{\circ}$C]')
175
176    ax1.grid(color='k', linestyle='-', linewidth=0.3)
177    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
178    ax1.annotate('Algo: '+L_ALGOS[ja]+', station: '+cforcing, xy=(0.5, 1.), xycoords='axes fraction', ha='center',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
179    plt.savefig('01_temperatures_'+L_ALGOS[ja]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
180    plt.close(1)
181
182
183del xxx
184
185
186
187# Now we compare output variables from bulk algorithms between them:
188
189nb_var = len(L_VNEM)
190
191xF  = nmp.zeros((nbr,nb_algos))
192xFa = nmp.zeros((nbr,nb_algos))
193
194
195for jv in range(nb_var):
196    print('\n *** Treating variable: '+L_VARO[jv]+' !')
197
198    for ja in range(nb_algos):
199        #
200        id_in = Dataset(cf_in[ja])
201        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
202        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
203        id_in.close()
204
205    idx_okay = nmp.where( nmp.abs(xF) < 1.e+10 )
206       
207    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
208    ax1 = plt.axes([0.08, 0.23, 0.9, 0.7])
209    ax1.set_xticks(vtime[::xticks_d])
210    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
211    plt.xticks(rotation='60', **font_x)
212
213    for ja in range(nb_algos):
214        plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
215       
216    fmin, fmax = round_bounds( nmp.min(xF[idx_okay]) , nmp.max(xF[idx_okay]), base=L_BASE[jv], prec=L_PREC[jv])
217    #print("LOLO: fmin, fmax =",nmp.min(xF[idx_okay]), nmp.max(xF[idx_okay]) ); print("LOLO: fmin, fmax =",fmin, fmax) ; #sys.exit(0)
218   
219    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
220    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
221
222    ax1.grid(color='k', linestyle='-', linewidth=0.3)
223    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
224    ax1.annotate(cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', ha='center',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
225                 zorder=50, **font_inf)
226    plt.savefig(L_VARO[jv]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
227    plt.close(jv)
228
229    if L_ANOM[jv]:
230
231        for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
232
233        if nmp.sum(xFa[:,:]) == 0.0:
234            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
235            print('          ==> skipping anomaly plot...')
236
237        else:
238
239            # Want a symetric y-range that makes sense for the anomaly we're looking at:
240            rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
241            rmax = max( abs(rmax) , abs(rmin) )
242            romagn = floor(log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
243            rmlt = 10.**(int(romagn)) / 2.
244            yrng = copysign( ceil(abs(rmax)/rmlt)*rmlt , rmax)
245
246            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
247            ax1 = plt.axes([0.08, 0.25, 0.9, 0.7])
248            ax1.set_xticks(vtime[::xticks_d])
249            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
250            plt.xticks(rotation='60', **font_x)
251
252            for ja in range(nb_algos):
253                plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
254
255            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
256            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
257            ax1.grid(color='k', linestyle='-', linewidth=0.3)
258            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
259            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)
260            plt.savefig(L_VARO[jv]+'_'+cforcing+'.'+fig_ext, dpi=int(rDPI), transparent=False)
261            plt.close(10+jv)
262
263
264
265
Note: See TracBrowser for help on using the repository browser.