source: NEMO/trunk/tests/STATION_ASF/EXPREF/plot_station_asf.py

Last change on this file was 13468, checked in by laurent, 6 weeks ago

Python2.7 to Python3!

  • Property svn:executable set to *
File size: 9.5 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, 2020
5
6import sys
7from os import path as path
8import math
9import numpy as nmp
10from netCDF4 import Dataset
11import matplotlib as mpl
12mpl.use('Agg')
13import matplotlib.pyplot as plt
14import matplotlib.dates as mdates
15
16cy1     = '2018' ; # First year
17cy2     = '2018' ; # Last year
18
19dir_figs='.'
20size_fig=(13,7)
21fig_ext='png'
22
23clr_red = '#AD0000'
24clr_sat = '#ffed00'
25clr_mod = '#008ab8'
26
27rDPI=200.
28
29L_ALGOS = [ 'COARE3p6' , 'ECMWF'   , 'NCAR' ]
30l_xtrns = [ '-noskin'  , '-noskin' ,  ''    ] ; # string to add to algo name (L_ALGOS) to get version without skin params turned on
31l_color = [  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot
32l_width = [     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot
33l_style = [    '-'     ,   '-'     , '--'   ] ; # line-style
34
35L_VNEM  = [   'qla'     ,     'qsb'     ,     'qt'     ,   'qlw'     ,  'taum'     ,    'dt_skin'         ]
36L_VARO  = [   'Qlat'    ,    'Qsen'     ,     'Qnet'   ,   'Qlw'     ,  'Tau'      ,    'dT_skin'         ] ; # name of variable on figure
37L_VARL  = [ r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$|\tau|$' , r'$\Delta T_{skin}$' ] ; # name of variable in latex mode
38L_VUNT  = [ r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$N/m^2$'  ,      'K'             ]
39L_VMAX  = [     75.     ,     75.       ,    800.      ,     25.     ,    1.2      ,       0.7            ]
40L_VMIN  = [   -250.     ,   -125.       ,   -400.      ,   -150.     ,    0.       ,      -0.7            ]
41L_ANOM  = [   True      ,    True       ,    True      ,    True     ,   True      ,      False           ]
42
43#L_VNEM  = [   'qlw'      ]
44#L_VARO  = [   'Qlw'      ] ; # name of variable on figure
45#L_VARL  = [ r'$Q_{lw}$'  ] ; # name of variable in latex mode
46#L_VUNT  = [ r'$W/m^2$'   ]
47#L_VMAX  = [     25.      ]
48#L_VMIN  = [   -150.      ]
49#L_ANOM  = [    True      ]
50
51
52
53nb_algos = len(L_ALGOS) ; print(nb_algos)
54
55# Getting arguments:
56narg = len(sys.argv)
57if narg != 2:
58    print('Usage: '+sys.argv[0]+' <DIR_OUT_SASF>'); sys.exit(0)
59cdir_data = sys.argv[1]
60
61
62
63# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
64# Populating and checking existence of files to be read
65# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
66def chck4f(cf):
67    cmesg = 'ERROR: File '+cf+' does not exist !!!'
68    if not path.exists(cf): print(cmesg); sys.exit(0)
69
70###cf_in = nmp.empty((), dtype="S10")
71cf_in = [] ; cf_in_ns = []
72for ja in range(nb_algos):
73    cfi = cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc'
74    chck4f(cfi)
75    cf_in.append(cfi)
76# Same but without skin params:
77for ja in range(nb_algos):
78    cfi = cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+l_xtrns[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc'
79    chck4f(cfi)
80    cf_in_ns.append(cfi)
81print('Files we are goin to use:')
82for ja in range(nb_algos): print(cf_in[ja])
83print('   --- same without cool-skin/warm-layer:')
84for ja in range(nb_algos): print(cf_in_ns[ja])
85#-----------------------------------------------------------------
86
87
88# Getting time array from the first file:
89id_in = Dataset(cf_in[0])
90vt = id_in.variables['time_counter'][:]
91cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
92id_in.close()
93nbr = len(vt)
94
95
96
97
98vtime = nmp.zeros(nbr)
99
100vt = vt + 1036800. + 30.*60. # BUG!??? don't get why false in epoch to date conversion, and yet ncview gets it right!
101for jt in range(nbr): vtime[jt] = mdates.epoch2num(vt[jt])
102
103ii=nbr/300
104ib=max(ii-ii%10,1)
105xticks_d=int(30*ib)
106
107font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':14 }
108
109nb_var = len(L_VNEM)
110
111xF  = nmp.zeros((nbr,nb_algos))
112xFa = nmp.zeros((nbr,nb_algos))
113
114
115for ctest in ['skin','noskin']:
116
117    for jv in range(nb_var):
118        print('\n *** Treating variable: '+L_VARO[jv]+' ('+ctest+') !')
119
120        for ja in range(nb_algos):
121            #
122            if ctest == 'skin':   id_in = Dataset(cf_in[ja])
123            if ctest == 'noskin': id_in = Dataset(cf_in_ns[ja])
124            xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
125            if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
126            id_in.close()
127
128        fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
129
130        ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
131
132        ax1.set_xticks(vtime[::xticks_d])
133        ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
134        plt.xticks(rotation='60')
135
136        for ja in range(nb_algos):
137            plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
138
139        ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
140        plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
141
142        ax1.grid(color='k', linestyle='-', linewidth=0.3)
143        plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
144        ax1.annotate(cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
145        plt.savefig(L_VARO[jv]+'_'+ctest+'.'+fig_ext, dpi=int(rDPI), transparent=False)
146        plt.close(jv)
147
148
149
150        if L_ANOM[jv]:
151
152            for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
153
154            if nmp.sum(xFa[:,:]) == 0.0:
155                print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
156                print('          ==> skipping anomaly plot...')
157
158            else:
159
160                # Want a symetric y-range that makes sense for the anomaly we're looking at:
161                rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
162                rmax = max( abs(rmax) , abs(rmin) )
163                romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
164                rmlt = 10.**(int(romagn)) / 2.
165                yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax)
166
167                fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
168                ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
169
170                ax1.set_xticks(vtime[::xticks_d])
171                ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
172                plt.xticks(rotation='60')
173
174                for ja in range(nb_algos):
175                    plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
176
177                ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
178                plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
179                ax1.grid(color='k', linestyle='-', linewidth=0.3)
180                plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
181                ax1.annotate('Anomaly of '+cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
182                plt.savefig(L_VARO[jv]+'_'+ctest+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False)
183                plt.close(10+jv)
184
185
186
187
188# Difference skin vs noskin:
189xFns = nmp.zeros((nbr,nb_algos))
190
191for jv in range(nb_var-1):
192    print('\n *** Treating variable: '+L_VARO[jv]+' ('+ctest+') !')
193
194    for ja in range(nb_algos-1):
195        id_in = Dataset(cf_in[ja])
196        xF[:,ja]   = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
197        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
198        id_in.close()
199        #
200        id_in = Dataset(cf_in_ns[ja])
201        xFns[:,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        xFa[:,ja] = xF[:,ja] - xFns[:,ja]  ; # difference!
206
207
208    # Want a symetric y-range that makes sense for the anomaly we're looking at:
209    rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
210    rmax = max( abs(rmax) , abs(rmin) )
211    romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
212    rmlt = 10.**(int(romagn)) / 2.
213    yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax)
214
215   
216    for ja in range(nb_algos-1):
217
218        calgo = L_ALGOS[ja]
219
220        if nmp.sum(xFa[:,ja]) == 0.0:
221            print('     Well! Seems that for variable '+L_VARO[jv]+', and algo '+calgo+', skin param has no impact')
222            print('          ==> skipping difference plot...')
223
224        else:
225
226            fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
227            ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
228
229            ax1.set_xticks(vtime[::xticks_d])
230            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
231            plt.xticks(rotation='60')
232
233            plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=None, zorder=10+ja)
234
235            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
236            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
237
238            ax1.grid(color='k', linestyle='-', linewidth=0.3)
239            #plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
240            ax1.annotate(cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
241            plt.savefig('diff_skin-noskin_'+L_VARO[jv]+'_'+calgo+'_'+ctest+'.'+fig_ext, dpi=int(rDPI), transparent=False)
242            plt.close(jv)
Note: See TracBrowser for help on using the repository browser.