-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathquarterly_arcmfc_report.py
executable file
·236 lines (179 loc) · 7.46 KB
/
quarterly_arcmfc_report.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
#!/usr/bin/env python
# ./validate 201409
import scipy as sp
import numpy as np
import pylab as pl
import netCDF4
import dataanalysis as da
import os
import datetime as dt
import validationtools as vt
from stationlist import ARCMFClocations as locations
from stationlist import WMsensors, bestWMsensor
import sys
from collectdatatools import validationfile
print("The Python version is %s.%s.%s" % sys.version_info[:3])
timeplist=['201710','201711','201712']
timestr='2017 Oct-Dec'
timestrt='2017_Oct-Dec'
print('time: '+timestr)
# plotpath
ppath = '/lustre/storeA/project/fou/hi/waveverification/Arc-MFC/'+timestrt+'/'
# set color table for models
ct = {'Subjective': 'b', 'WAM10': 'c', 'WAM4':'m', 'ECWAM':'k', 'LAWAM':'0.25', 'AROME': 'b', 'HIRLAM8': 'y', 'MWAM4':'r', 'EXP':'y', 'MWAM4exp':'w', 'MWAM10':'w', 'MWAM8':'g'}
def select_var_from_models(vf,varname):
modeldata={}
for gname in vf.models:
try:
mod = vf.get_modelgroup(gname,create=False)
varraw = mod[varname][:]
try:
var = varraw.data
var[varraw.mask==True]=sp.nan
except AttributeError:
var = varraw
#
# fill data gap in second forecast range
#
var[1] = 0.5*(var[0]+var[2])
#
# check if we are dealing with directions and ensure meteorological convention
#if (G[gname].variables[varname].units[0:6] == 'degree'):
# try:
# if (G[gname].variables[varname].Convention=='oceanographic'):
# var=var+180
# var[var>360.]=var[var>360.]-360.
# except AttributeError:
# var=var
except KeyError:
continue
if sp.isnan(var[0]).all():
continue
modeldata.update({gname: var})
return modeldata
from matplotlib import dates
minorLocator=dates.DayLocator(range(33))
majorLocator=dates.DayLocator(range(5,31,5))
fmt=dates.DateFormatter('%d.%m.%Y')
varname = 'Hs'
obs_all = []
mod_all = {'MWAM8':[]}
time_all = []
os.system('mkdir -p '+ppath)
for station, parameters in locations.iteritems():
print(' ')
print('verification of station '+station+' for '+timestr)
obs_long = []
mod_long = {'MWAM8':[]}
for timep in timeplist: # loop over months
print ' '
print 'read data for station '+station+' for '+timep
# open file
path = '/lustre/storeA/project/fou/hi/waveverification/data'
year,month = int(timep[0:4]),int(timep[4:6])
vf = validationfile(path,station,year,month)
time = vf.time
OBS = vf.get_obs()
if station=='draugen':
time_all = time_all + list(time)
# Specify which WM sensor to use for validation
try:
sensor = bestWMsensor[station]
except KeyError:
sensor = 0
obsraw = OBS[varname][sensor]
try:
obs = obsraw.data
obs[obsraw.mask==True] = sp.nan # make sure all masked values are nan
except AttributeError:
obs = obsraw
units = vf.nc.variables[varname+'_OBS'].units
if all(sp.isnan(obs.data)):
print('no data for '+station+' during '+timestr)
continue
# select variable from each model:
modeldata = select_var_from_models(vf,varname)
#
#nc.close()
#
obs_long.append(obs)
for gname, var in modeldata.iteritems():
if gname in mod_long.keys():
mod_long[gname].append(var)
# make arrays from lists
obs_long = np.concatenate(obs_long)
mod_longa={}
for mod in mod_long.keys():
mod_longa[mod] = np.concatenate(mod_long[mod],axis=1)
print '%s, %s' % (station, str(mod_long['MWAM8'][1].shape))
# append to list for all stations:
obs_all.append(obs_long)
for gname, var in modeldata.iteritems():
if gname in mod_all.keys():
print('append ' +gname+ ' for ' + station)
mod_all[gname].append(mod_longa[gname])
# make arrays from list
obs = np.array(obs_all)
modeldata = {}
for gname in mod_all.keys():
modeldata[gname] = np.dstack(mod_all[gname])
var = modeldata['MWAM8'].transpose(2,1,0)
print('var.shape ', var.shape)
print('obs.shape ', obs.shape)
#
# compute statistics
#
# reshape time axis to join hours
nhours = 6
obsS = obs.reshape(obs.shape[0],nhours,obs.shape[1]/nhours)
varS = var.reshape(var.shape[0],nhours,var.shape[1]/nhours,var.shape[2])
time = time_all[::6]
timeunit = "days since 2001-01-01 12:00:00 UTC"
time_start = time[0].strftime('%Y%m%d')
time_end = time[-1].strftime('%Y%m%d')
print(time_start, time_end)
# produce netcdf file:
nc = netCDF4.Dataset(os.path.join(ppath,'product_quality_stats_ARCTIC_ANALYSIS_FORECAST_WAV_002_006_'+time_start+'-'+time_end+'.nc'),'w')
nc.contact = '[email protected]'
nc.product = 'Arctic wave model WAM'
nc.production_centre = 'Arctic MFC'
nc.production_unit = 'Norwegian Meteorological Institute'
nc.creation_date = str(dt.datetime.now())
nc.thredds_web_site = 'http://thredds.met.no/thredds/myocean/ARC-MFC/mywave-arctic.html'
ncdims = {'string_length':28, 'areas':1, 'metrics':4, 'surface':1, 'forecasts':10} # and time, unlim
metric_names = [name.ljust(28) for name in ["mean of product", "mean of reference", "mean square difference", "number of data values"]]
nc.createDimension('time',size=None)
for name,dim in ncdims.iteritems():
nc.createDimension(name,size=dim)
nc_time = nc.createVariable('time','f8',dimensions=('time'))
nc_time[:] = netCDF4.date2num(time,units=timeunit)
nc_time.units = timeunit
nc_time.long_name = 'validity time'
nc_metricnames = nc.createVariable('metric_names','S1', dimensions=(u'metrics',u'string_length'))
nc_metricnames[:] = netCDF4.stringtochar(np.array(metric_names))
nc_areanames = nc.createVariable('area_names','S1', dimensions=(u'areas',u'string_length'))
nc_areanames[0] = netCDF4.stringtochar(np.array('North Sea and Norwegian Sea'.ljust(28)))
nc_areanames.long_name = 'area names'
nc_areanames.description = 'region over which statistics are aggregated'
nc_leadtime = nc.createVariable('forecasts','f4',dimensions=('forecasts'))
nc_leadtime.long_name = 'forecast lead time'
nc_leadtime.units = 'hours'
nc_leadtime[:] = sp.arange(12,229,24)
varArcMFCname = {'Hs': 'stats_VHM0'}
varstandardname = {'Hs':'sea_surface_wave_significant_height'}
ncvar = nc.createVariable(varArcMFCname[varname],'f4', dimensions=('time', 'forecasts', 'surface', 'metrics', 'areas'))
ncvar.standard_name = varstandardname[varname]
ncvar.parameter = varArcMFCname[varname]
ncvar.units = 'm'
ncvar.reference_source = 'wave data from offshore platforms available from d22 files at the Norwegian Meteorological Institute'
# calculate statistics (bias and root-mean-square difference)
for leadtime in range(10):
# mean of product
ncvar[:,leadtime,0,0,0] = sp.array([sp.mean( vt.returnclean(obsS[:,:,i],varS[:,:,i,leadtime])[1]) for i in range(obsS.shape[2])])
# mean of reference
ncvar[:,leadtime,0,1,0] = sp.array([sp.mean( vt.returnclean(obsS[:,:,i],varS[:,:,i,leadtime])[0]) for i in range(obsS.shape[2])])
# mean square difference
ncvar[:,leadtime,0,2,0] = sp.array([vt.msd( *vt.returnclean(obsS[:,:,i],varS[:,:,i,leadtime])) for i in range(obsS.shape[2])])
# number of data values
ncvar[:,leadtime,0,3,0] = sp.array([ len(vt.returnclean(obsS[:,:,i],varS[:,:,i,leadtime])[1]) for i in range(obsS.shape[2])])
nc.close()