source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/plot_station_asf.py @ 11930

Last change on this file since 11930 was 11930, checked in by laurent, 10 months ago

New functional "stp_ctl" for STATION_ASF, test-case launch script and python diagnostics…

  • Property svn:executable set to *
File size: 5.6 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
7#from 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     = '2016' ; # 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_color = [  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot
43l_width = [     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot
44l_style = [    '-'     ,   '-'     , '--'   ] ; # line-style
45
46L_VNEM  = [ 'qla'  , 'qsb'   , 'qt'   , 'qlw' , 'taum'  , 'dt_skin'  ]
47L_VARO  = [  'Qlat'    ,   'Qsen'    ,  'Qnet'    ,  'Qlw'    , 'Tau'   , 'dT_skin'  ] ; # name of variable on figure
48L_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
49L_VUNT  = [  r'$W/m^2$'   ,  r'$W/m^2$'    ,  r'$W/m^2$'   ,  r'$W/m^2$'  , r'$N/m^2$' ,   'K'     ]
50L_VMAX  = [    75.     ,    75.      ,   800.     ,    25.    ,   1.2   ,   -0.7     ]
51L_VMIN  = [  -250.     ,  -125.      ,  -400.     ,  -150.    ,   0.    ,    0.7     ]
52L_ANOM  = [  True      ,   True      ,   True     ,   True    ,  True   ,  False     ]
53
54nb_algos = len(L_ALGOS) ; print(nb_algos)
55
56# Getting arguments:
57narg = len(sys.argv)
58if narg != 3:
59    print 'Usage: '+sys.argv[0]+' <DIR_OUT_SASF> <year>'; sys.exit(0)
60cdir_in = sys.argv[1]
61cyyyy   = sys.argv[2]
62
63
64
65# Files to be read:
66cf_in = []
67for ja in range(nb_algos): cf_in.append(cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc')
68print('Files we are goin to use:')
69for ja in range(nb_algos): print(cf_in[ja])
70
71
72# Getting time array from the first file:
73id_in = Dataset(cf_in[0])
74vt = id_in.variables['time_counter'][:]
75cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
76id_in.close()
77nbr = len(vt)
78
79
80
81
82vtime = nmp.zeros(nbr)
83
84vt = vt + 1036800. # BUG!??? don't get why false in epoch to date conversion, and yet ncview gets it right!
85for jt in range(nbr): vtime[jt] = mdates.epoch2num(vt[jt])
86
87ii=nbr/300
88ib=max(ii-ii%10,1)
89xticks_d=int(30*ib)
90
91font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':14 }
92
93nb_var = len(L_VNEM)
94
95xF  = nmp.zeros((nbr,nb_algos))
96xFa = nmp.zeros((nbr,nb_algos))
97
98
99
100
101
102for jv in range(nb_var):
103    print('\n *** Treating variable: '+L_VARO[jv]+' !')
104
105    for ja in range(nb_algos):
106        id_in = Dataset(cf_in[ja])
107        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain!
108        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
109        id_in.close()   
110        #print('    => '+L_ALGOS[ja]+'  => ok!')
111
112
113    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
114
115    ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
116
117    ax1.set_xticks(vtime[::xticks_d])
118    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
119    plt.xticks(rotation='60')
120
121    for ja in range(nb_algos):
122        plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
123
124    ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
125    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
126   
127    ax1.grid(color='k', linestyle='-', linewidth=0.3)
128    plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
129    ax1.annotate(cvar_lnm, xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf)
130    plt.savefig(L_VARO[jv]+'.'+fig_ext, dpi=int(rDPI), transparent=False)
131    plt.close(jv)
132
133
134
135    if L_ANOM[jv]:
136
137        for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
138
139        # Want a symetric y-range that makes sense for the anomaly we're looking at:
140        rmax = nmp.max(xFa) ; rmin = nmp.min(xFa)
141        rmax = max( abs(rmax) , abs(rmin) )
142        romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
143        rmlt = 10.**(int(romagn)) / 2.
144        yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax)
145        #print 'yrng = ', yrng ;  #sys.exit(0)
146       
147        fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
148        ax1 = plt.axes([0.07, 0.22, 0.9, 0.75])
149       
150        ax1.set_xticks(vtime[::xticks_d])
151        ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
152        plt.xticks(rotation='60')
153       
154        for ja in range(nb_algos):
155            plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja)
156
157        ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1])
158        plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
159        ax1.grid(color='k', linestyle='-', linewidth=0.3)
160        plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
161        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)
162        plt.savefig(L_VARO[jv]+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False)
163        plt.close(10+jv)
164
Note: See TracBrowser for help on using the repository browser.