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

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

Setting up the new STATION_ASF with sea-ice support

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