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

Last change on this file since 13690 was 13690, checked in by laurent, 6 months ago

Improvements for "oce+ice capable" STATION_ASF!

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