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

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

Last change on this file since 13264 was 13264, checked in by laurent, 23 months ago

Slight improvements in STATION_ASF test-case

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