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

Last change on this file since 13675 was 13675, checked in by laurent, 7 months ago

Setting up the new STATION_ASF with sea-ice support

  • Property svn:executable set to *
File size: 7.4 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
16CONFIG='STATION_ASF'
17FORCING='ERA5_arctic'
18# (files are: output/3x3/'+CONFIG+'-'+algo+'_1h_'+year+'0101_'+year+'1231_gridT_'+FORCING+'.nc' )
19cstation = 'ERA5 81N, 36.75E'
20
21cy1     = '2018' ; # First year
22cy2     = '2018' ; # Last year
23
24jt0 = 0
25
26dir_figs='.'
27size_fig=(13,8)
28size_fig0=(12,10)
29fig_ext='png'
30
31clr_red = '#AD0000'
32clr_sat = '#ffed00'
33clr_mod = '#008ab8'
34
35rDPI=100.
36
37L_ALGOS = [ 'ECMWF-LG15', 'ECMWF-LU12', 'ECMWF-CSTC' ]
38l_color = [  '#ffed00'  , '#008ab8'   ,     '0.4'    ] ; # colors to differentiate algos on the plot
39l_width = [     3       ,      2      ,       1      ] ; # line-width to differentiate algos on the plot
40l_style = [    '-'      ,     '-'     ,      '--'    ] ; # line-style
41
42
43# Variables to compare for A GIVEN algorithm
44###############################################
45#L_VNEM0 = [
46
47
48
49# Variables to compare between algorithms
50############################################
51L_VNEM = [   'Cd_ice',   'Ce_ice',   'qla_ice' ,   'qsb_ice'   ,  'qt_ice'    ,  'qlw_ice'  ,  'qsr_ice'  , 'taum_ai'   ]
52L_VARO = [     'Cd'  ,     'Ce'  ,   'Qlat'    ,    'Qsen'     ,   'Qnet'     ,   'Qlw'     ,    'Qsw'    ,  'Tau'      ]
53L_VARL = [ r'$C_{D}$', r'$C_{E}$', r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$Q_{sw}$' , r'$|\tau|$' ]
54L_VUNT = [     ''    ,     ''    , r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$W/m^2$'  , r'$N/m^2$'  ]
55L_VMAX = [    0.003  ,    0.003  ,     75.     ,     75.       ,    800.      ,    200.     ,    200.     ,    1.2      ]
56L_VMIN = [    0.     ,    0.     ,   -250.     ,   -125.       ,   -400.      ,   -200.     ,      0.     ,    0.       ]
57L_ANOM = [   False   ,   False   ,   True      ,    True       ,    True      ,    True     ,    True     ,   True      ]
58
59
60nb_algos = len(L_ALGOS)
61
62# Getting arguments:
63narg = len(sys.argv)
64if narg != 2:
65    print('Usage: '+sys.argv[0]+' <DIR_OUT_SASF>'); sys.exit(0)
66cdir_data = sys.argv[1]
67
68
69
70# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
71# Populating and checking existence of files to be read
72# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
73def chck4f(cf):
74    cmesg = 'ERROR: File '+cf+' does not exist !!!'
75    if not path.exists(cf): print(cmesg) ; sys.exit(0)
76
77cf_in = []
78for ja in range(nb_algos):
79    cfi = cdir_data+'/output/3x3/'+CONFIG+'-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_icemod_'+FORCING+'.nc'
80    chck4f(cfi)
81    cf_in.append(cfi)
82print('Files we are goin to use:')
83for ja in range(nb_algos): print(cf_in[ja])
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()
91Nt = len(vt)
92
93vtime = num2date(vt, units=cunit_t) ; # something understandable!
94vtime = vtime.astype(dtype='datetime64[D]')
95
96ii=Nt/300
97ib=max(ii-ii%10,1)
98xticks_d=int(30*ib)
99
100rat = 100./float(rDPI)
101params = { 'font.family':'Open Sans',
102           'font.size':       int(15.*rat),
103           'legend.fontsize': int(15.*rat),
104           'xtick.labelsize': int(15.*rat),
105           'ytick.labelsize': int(15.*rat),
106           'axes.labelsize':  int(16.*rat)
107}
108mpl.rcParams.update(params)
109font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':18.*rat }
110font_x   = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':15.*rat }
111
112
113# Now we compare output variables from bulk algorithms between them:
114
115nb_var = len(L_VNEM)
116
117xF  = nmp.zeros((Nt,nb_algos))
118xFa = nmp.zeros((Nt,nb_algos))
119
120
121for jv in range(nb_var):
122    print('\n *** Treating variable: '+L_VARO[jv]+' !')
123
124    for ja in range(nb_algos):
125        #
126        id_in = Dataset(cf_in[ja])
127        xF[:,ja] = id_in.variables[L_VNEM[jv]][jt0:,1,1] # only the center point of the 3x3 spatial domain!
128        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
129        id_in.close()
130
131    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
133    ax1 = plt.axes([0.08, 0.25, 0.9, 0.7])
134    ax1.set_xticks(vtime[::xticks_d])
135    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
136    plt.xticks(rotation='60', **font_x)
137
138    for ja in range(nb_algos):
139        plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
140
141    ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
142    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
143
144    ax1.grid(color='k', linestyle='-', linewidth=0.3)
145    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
146    ax1.annotate(cvar_lnm+', station: '+cstation, xy=(0.3, 1.), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
147    plt.savefig(L_VARO[jv]+'.'+fig_ext, dpi=int(rDPI), transparent=False)
148    plt.close(jv)
149    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
150
151    def symetric_range( pmin, pmax ):
152        # Returns a symetric f-range that makes sense for the anomaly of "f" we're looking at...
153        from math import floor, copysign, log, ceil
154        zmax = max( abs(pmax) , abs(pmin) )
155        romagn = floor(log(zmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
156        rmlt = 10.**(int(romagn)) / 2.
157        frng = copysign( ceil(abs(zmax)/rmlt)*rmlt , zmax)
158        return frng
159
160   
161
162    if L_ANOM[jv]:
163
164        for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
165
166        if nmp.sum(nmp.abs(xFa[:,:])) == 0.0:
167            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
168            print('          ==> skipping anomaly plot...')
169
170        else:
171
172            yrng = symetric_range( nmp.min(xFa) , nmp.max(xFa) )
173
174            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
175            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
176            ax1 = plt.axes([0.08, 0.25, 0.9, 0.7])
177            ax1.set_xticks(vtime[::xticks_d])
178            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
179            plt.xticks(rotation='60', **font_x)
180
181            for ja in range(nb_algos):
182                plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
183
184            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
185            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
186            ax1.grid(color='k', linestyle='-', linewidth=0.3)
187            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
188            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)
189            plt.savefig(L_VARO[jv]+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False)
190            plt.close(10+jv)
191            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
192
193
194
195
Note: See TracBrowser for help on using the repository browser.