Context Navigation

source: flex_extract.git/source/python/classes/EcFlexpart.py @ 27fe969

ctbtodev

Last change on this file since 27fe969 was 27fe969, checked in by Anne Philipp <anne.philipp@…>, 6 years ago
put grib2flexpart part of code into its own function; changed function documentation of input controlfile parameter
Property mode set to `100644`
File size: 52.7 KB

Line
1	#!/usr/bin/env python
2	# -- coding: utf-8 --
3	#*******************************************************************************
4	# @Author: Anne Fouilloux (University of Oslo)
5	#
6	# @Date: October 2014
7	#
8	# @Change History:
9	#
10	# November 2015 - Leopold Haimberger (University of Vienna):
11	# - extended with class Control
12	# - removed functions mkdir_p, daterange, years_between, months_between
13	# - added functions darain, dapoly, to_param_id, init128, normal_exit,
14	# my_error, clean_up, install_args_and_control,
15	# interpret_args_and_control,
16	# - removed function __del__ in class EIFLexpart
17	# - added the following functions in EIFlexpart:
18	# - create_namelist
19	# - process_output
20	# - deacc_fluxes
21	# - modified existing EIFlexpart - functions for the use in
22	# flex_extract
23	# - retrieve also longer term forecasts, not only analyses and
24	# short term forecast data
25	# - added conversion into GRIB2
26	# - added conversion into .fp format for faster execution of FLEXPART
27	# (see https://www.flexpart.eu/wiki/FpCtbtoWo4FpFormat)
28	#
29	# February 2018 - Anne Philipp (University of Vienna):
30	# - applied PEP8 style guide
31	# - added documentation
32	# - removed function getFlexpartTime in class EcFlexpart
33	# - outsourced class ControlFile
34	# - outsourced class MarsRetrieval
35	# - changed class name from EIFlexpart to EcFlexpart
36	# - applied minor code changes (style)
37	# - removed "dead code" , e.g. retrieval of Q since it is not needed
38	# - removed "times" parameter from retrieve-method since it is not used
39	# - seperated function "retrieve" into smaller functions (less code
40	# duplication, easier testing)
41	#
42	# @License:
43	# (C) Copyright 2014-2018.
44	#
45	# This software is licensed under the terms of the Apache Licence Version 2.0
46	# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
47	#
48	# @Class Description:
49	# FLEXPART needs grib files in a specifc format. All necessary data fields
50	# for one time step are stored in a single file. The class represents an
51	# instance with all the parameter and settings necessary for retrieving
52	# MARS data and modifing them so they are fitting FLEXPART need. The class
53	# is able to disaggregate the fluxes and convert grid types to the one needed
54	# by FLEXPART, therefore using the FORTRAN program.
55	#
56	# @Class Content:
57	# - __init__
58	# - write_namelist
59	# - retrieve
60	# - process_output
61	# - create
62	# - deacc_fluxes
63	#
64	# @Class Attributes:
65	#
66	# TODO
67	#
68	#*******************************************************************************
69	#pylint: disable=unsupported-assignment-operation
70	# this is disabled because for this specific case its an error in pylint
71	#pylint: disable=consider-using-enumerate
72	# this is not useful in this case
73	# ------------------------------------------------------------------------------
74	# MODULES
75	# ------------------------------------------------------------------------------
76	import subprocess
77	import shutil
78	import os
79	import glob
80	from datetime import datetime, timedelta
81	import numpy as np
82	from gribapi import grib_set, grib_index_select, grib_new_from_index, grib_get,\
83	grib_write, grib_get_values, grib_set_values, grib_release,\
84	grib_index_release, grib_index_get
85
86	# software specific classes and modules from flex_extract
87	import _config
88	from GribTools import GribTools
89	from mods.tools import init128, to_param_id, silent_remove, product, my_error
90	from MarsRetrieval import MarsRetrieval
91	import mods.disaggregation
92
93	# ------------------------------------------------------------------------------
94	# CLASS
95	# ------------------------------------------------------------------------------
96	class EcFlexpart(object):
97	'''
98	Class to retrieve FLEXPART specific ECMWF data.
99	'''
100	# --------------------------------------------------------------------------
101	# CLASS FUNCTIONS
102	# --------------------------------------------------------------------------
103	def __init__(self, c, fluxes=False):
104	'''
105	@Description:
106	Creates an object/instance of EcFlexpart with the
107	associated settings of its attributes for the retrieval.
108
109	@Input:
110	self: instance of EcFlexpart
111	The current object of the class.
112
113	c: instance of class ControlFile
114	Contains all the parameters of CONTROL file and
115	command line.
116	For more information about format and content of the parameter
117	see documentation.
118
119	fluxes: boolean, optional
120	Decides if the flux parameter settings are stored or
121	the rest of the parameter list.
122	Default value is False.
123
124	@Return:
125	<nothing>
126	'''
127	# set a counter for the number of mars requests generated
128	self.mreq_count = 0
129
130	# different mars types for retrieving data for flexpart
131	self.types = dict()
132
133	if c.maxstep > len(c.type): # Pure forecast mode
134	c.type = [c.type[1]]
135	c.step = ['{:0>3}'.format(int(c.step[0]))]
136	c.time = [c.time[0]]
137	for i in range(1, c.maxstep + 1):
138	c.type.append(c.type[0])
139	c.step.append('{:0>3}'.format(i))
140	c.time.append(c.time[0])
141
142	self.inputdir = c.inputdir
143	self.basetime = c.basetime
144	self.dtime = c.dtime
145	i = 0
146	if fluxes and c.maxstep <= 24:
147	# no forecast beyond one day is needed!
148	# Thus, prepare flux data manually as usual
149	# with only forecast fields with start times at 00/12
150	# (but without 00/12 fields since these are
151	# the initialisation times of the flux fields
152	# and therefore are zero all the time)
153	self.types[c.type[1]] = {'times': '00/12', 'steps':
154	'{}/to/12/by/{}'.format(c.dtime, c.dtime)}
155	else:
156	for ty, st, ti in zip(c.type, c.step, c.time):
157	btlist = range(24)
158	if c.basetime == '12':
159	btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
160	if c.basetime == '00':
161	btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
162
163	if i % int(c.dtime) == 0 and (i in btlist or c.maxstep > 24):
164
165	if ty not in self.types.keys():
166	self.types[ty] = {'times': '', 'steps': ''}
167
168	if ti not in self.types[ty]['times']:
169	if self.types[ty]['times']:
170	self.types[ty]['times'] += '/'
171	self.types[ty]['times'] += ti
172
173	if st not in self.types[ty]['steps']:
174	if self.types[ty]['steps']:
175	self.types[ty]['steps'] += '/'
176	self.types[ty]['steps'] += st
177	i += 1
178
179	self.marsclass = c.marsclass
180	self.stream = c.stream
181	self.number = c.number
182	self.resol = c.resol
183	self.accuracy = c.accuracy
184	self.level = c.level
185	self.expver = c.expver
186	self.levelist = c.levelist
187	# for gaussian grid retrieval
188	self.glevelist = '1/to/' + c.level
189
190	if hasattr(c, 'gaussian') and c.gaussian:
191	self.gaussian = c.gaussian
192	else:
193	self.gaussian = ''
194
195	if 'N' in c.grid: # Gaussian output grid
196	self.grid = c.grid
197	self.area = 'G'
198	else:
199	self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.)
200	self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000.,
201	int(c.left) / 1000.,
202	int(c.lower) / 1000.,
203	int(c.right) / 1000.)
204
205	self.outputfilelist = []
206
207
208	# Now comes the nasty part that deals with the different
209	# scenarios we have:
210	# 1) Calculation of etadot on
211	# a) Gaussian grid
212	# b) Output grid
213	# c) Output grid using parameter 77 retrieved from MARS
214	# 3) Calculation/Retrieval of omega
215	# 4) Download also data for WRF
216
217	# Different grids need different retrievals
218	# SH = Spherical Harmonics, GG = Gaussian Grid,
219	# OG = Output Grid, ML = MultiLevel, SL = SingleLevel
220	self.params = {'SH__ML': '', 'SH__SL': '',
221	'GG__ML': '', 'GG__SL': '',
222	'OG__ML': '', 'OG__SL': '',
223	'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
224	# the self.params dictionary stores a list of
225	# [param, levtype, levelist, grid] per key
226
227	if fluxes is False:
228	self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
229	# "SD/MSL/TCC/10U/10V/2T/2D/129/172"
230	self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
231	'SFC', '1', self.grid]
232	if c.addpar:
233	if c.addpar[0] == '/':
234	c.addpar = c.addpar[1:]
235	self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
236
237	self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
238	'SFC', '1', self.grid]
239
240	self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
241
242	#if c.gauss == '0' and c.eta == '1':
243	if not c.gauss and c.eta:
244	# the simplest case
245	self.params['OG__ML'][0] += '/U/V/77'
246	#elif c.gauss == '0' and c.eta == '0':
247	elif not c.gauss and not c.eta:
248	# this is not recommended (inaccurate)
249	self.params['OG__ML'][0] += '/U/V'
250	#elif c.gauss == '1' and c.eta == '0':
251	elif c.gauss and not c.eta:
252	# this is needed for data before 2008, or for reanalysis data
253	self.params['GG__SL'] = ['Q', 'ML', '1', \
254	'{}'.format((int(self.resol) + 1) / 2)]
255	self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
256	else:
257	print('Warning: This is a very costly parameter combination, \
258	use only for debugging!')
259	self.params['GG__SL'] = ['Q', 'ML', '1', \
260	'{}'.format((int(self.resol) + 1) / 2)]
261	self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
262	'{}'.format((int(self.resol) + 1) / 2)]
263
264	if c.omega:
265	self.params['OG__ML'][0] += '/W'
266
267	# add cloud water content if necessary
268	if c.cwc:
269	self.params['OG__ML'][0] += '/CLWC/CIWC'
270
271	# add vorticity and geopotential height for WRF if necessary
272	if c.wrf:
273	self.params['OG__ML'][0] += '/Z/VO'
274	if '/D' not in self.params['OG__ML'][0]:
275	self.params['OG__ML'][0] += '/D'
276	#wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
277	# stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
278	wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
279	139/170/183/236/39/40/41/42'
280	lwrt_sfc = wrf_sfc.split('/')
281	for par in lwrt_sfc:
282	if par not in self.params['OG__SL'][0]:
283	self.params['OG__SL'][0] += '/' + par
284
285	else:
286	self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
287	'SFC', '1', self.grid]
288
289	# if needed, add additional WRF specific parameters here
290
291	return
292
293
294	def _mk_targetname(self, ftype, param, date):
295	'''
296	@Description:
297	Creates the filename for the requested grib data to be stored in.
298	This name is passed as the "target" parameter in the request.
299
300	@Input:
301	ftype: string
302	Shortcut name of the type of the field. E.g. AN, FC, PF, ...
303
304	param: string
305	Shortcut of the grid type. E.g. SH__ML, SH__SL, GG__ML,
306	GG__SL, OG__ML, OG__SL, OG_OROLSM_SL, OG_acc_SL
307
308	date: string
309	The date period of the grib data to be stored in this file.
310
311	@Return:
312	targetname: string
313	The target filename for the grib data.
314	'''
315	targetname = (self.inputdir + '/' + ftype + param + '.' + date + '.' +
316	str(os.getppid()) + '.' + str(os.getpid()) + '.grb')
317
318	return targetname
319
320
321	def _start_retrievement(self, request, par_dict):
322	'''
323	@Description:
324	Creates the Mars Retrieval and prints or submits the request
325	depending on the status of the request variable.
326
327	@Input:
328	self: instance of EcFlexpart
329	The current object of the class.
330
331	request: integer
332	Selects the mode of retrieval.
333	0: Retrieves the data from ECMWF.
334	1: Prints the mars requests to an output file.
335	2: Retrieves the data and prints the mars request.
336
337	par_dict: dictionary
338	Contains all parameter which have to be set for creating the
339	Mars Retrievals. The parameter are:
340	marsclass, stream, type, levtype, levelist, resol, gaussian,
341	accuracy, grid, target, area, date, time, number, step, expver,
342	param
343
344	@Return:
345	<nothing>
346	'''
347	# increase number of mars requests
348	self.mreq_count += 1
349
350	MR = MarsRetrieval(self.server,
351	marsclass=par_dict['marsclass'],
352	stream=par_dict['stream'],
353	type=par_dict['type'],
354	levtype=par_dict['levtype'],
355	levelist=par_dict['levelist'],
356	resol=par_dict['resol'],
357	gaussian=par_dict['gaussian'],
358	accuracy=par_dict['accuracy'],
359	grid=par_dict['grid'],
360	target=par_dict['target'],
361	area=par_dict['area'],
362	date=par_dict['date'],
363	time=par_dict['time'],
364	number=par_dict['number'],
365	step=par_dict['step'],
366	expver=par_dict['expver'],
367	param=par_dict['param'])
368
369	if request == 0:
370	MR.display_info()
371	MR.data_retrieve()
372	elif request == 1:
373	MR.print_infodata_csv(self.inputdir, self.mreq_count)
374	elif request == 2:
375	MR.print_infodata_csv(self.inputdir, self.mreq_count)
376	MR.display_info()
377	MR.data_retrieve()
378	else:
379	print('Failure')
380
381	return
382
383
384	def retrieve(self, server, dates, request, inputdir='.'):
385	'''
386	@Description:
387	Finalizing the retrieval information by setting final details
388	depending on grid type.
389	Prepares MARS retrievals per grid type and submits them.
390
391	@Input:
392	self: instance of EcFlexpart
393	The current object of the class.
394
395	server: instance of ECMWFService or ECMWFDataServer
396	The connection to the ECMWF server. This is different
397	for member state users which have full access and non
398	member state users which have only access to the public
399	data sets. The decision is made from command line argument
400	"public"; for public access its True (ECMWFDataServer)
401	for member state users its False (ECMWFService)
402
403	dates: string
404	Contains start and end date of the retrieval in the format
405	"YYYYMMDD/to/YYYYMMDD"
406
407	request: integer
408	Selects the mode of retrieval.
409	0: Retrieves the data from ECMWF.
410	1: Prints the mars requests to an output file.
411	2: Retrieves the data and prints the mars request.
412
413	inputdir: string, optional
414	Path to the directory where the retrieved data is about
415	to be stored. The default is the current directory ('.').
416
417	@Return:
418	<nothing>
419	'''
420	self.dates = dates
421	self.server = server
422	self.inputdir = inputdir
423	oro = False
424
425	# define times with datetime module
426	t12h = timedelta(hours=12)
427	t24h = timedelta(hours=24)
428
429	# dictionary which contains all parameter for the mars request
430	# entries with a "None" will change in different requests and will
431	# therefore be set in each request seperately
432	retr_param_dict = {'marsclass':self.marsclass,
433	'stream':None,
434	'type':None,
435	'levtype':None,
436	'levelist':None,
437	'resol':self.resol,
438	'gaussian':None,
439	'accuracy':self.accuracy,
440	'grid':None,
441	'target':None,
442	'area':None,
443	'date':None,
444	'time':None,
445	'number':self.number,
446	'step':None,
447	'expver':self.expver,
448	'param':None}
449
450	for ftype in self.types:
451	# fk contains field types such as
452	# [AN, FC, PF, CV]
453	# fv contains all of the items of the belonging key
454	# [times, steps]
455	for pk, pv in self.params.iteritems():
456	# pk contains one of these keys of params
457	# [SH__ML, SH__SL, GG__ML, GG__SL, OG__ML, OG__SL,
458	# OG_OROLSM_SL, OG_acc_SL]
459	# pv contains all of the items of the belonging key
460	# [param, levtype, levelist, grid]
461	if isinstance(pv, str):
462	continue
463	retr_param_dict['type'] = '' + ftype
464	retr_param_dict['time'] = self.types[ftype]['times']
465	retr_param_dict['step'] = self.types[ftype]['steps']
466	retr_param_dict['date'] = self.dates
467	retr_param_dict['stream'] = self.stream
468	retr_param_dict['target'] = \
469	self._mk_targetname(ftype, pk,
470	retr_param_dict['date'].split('/')[0])
471	retr_param_dict['param'] = pv[0]
472	retr_param_dict['levtype'] = pv[1]
473	retr_param_dict['levelist'] = pv[2]
474	retr_param_dict['grid'] = pv[3]
475	retr_param_dict['area'] = self.area
476	retr_param_dict['gaussian'] = self.gaussian
477
478	if pk == 'OG__SL':
479	pass
480	if pk == 'OG_OROLSM__SL' and not oro:
481	oro = True
482	retr_param_dict['stream'] = 'OPER'
483	retr_param_dict['type'] = 'AN'
484	retr_param_dict['time'] = '00'
485	retr_param_dict['step'] = '000'
486	retr_param_dict['date'] = self.dates.split('/')[0]
487	retr_param_dict['target'] = self._mk_targetname('',
488	pk, retr_param_dict['date'])
489	elif pk == 'OG_OROLSM__SL' and oro:
490	continue
491	if pk == 'GG__SL' and pv[0] == 'Q':
492	retr_param_dict['area'] = ""
493	retr_param_dict['gaussian'] = 'reduced'
494
495	# ------ on demand path --------------------------------------------------
496	if not self.basetime:
497	# ******* start retrievement
498	self._start_retrievement(request, retr_param_dict)
499	# ------ operational path ------------------------------------------------
500	else:
501	# check if mars job requests fields beyond basetime.
502	# if yes eliminate those fields since they may not
503	# be accessible with user's credentials
504
505	enddate = retr_param_dict['date'].split('/')[-1]
506	elimit = datetime.strptime(enddate + self.basetime,
507	'%Y%m%d%H')
508
509	if self.basetime == '12':
510	# -------------- flux data ----------------------------
511	if 'acc' in pk:
512
513	# Strategy:
514	# if maxtime-elimit >= 24h reduce date by 1,
515	# if 12h <= maxtime-elimit<12h reduce time for last date
516	# if maxtime-elimit<12h reduce step for last time
517	# A split of the MARS job into 2 is likely necessary.
518
519
520	startdate = retr_param_dict['date'].split('/')[0]
521	enddate = datetime.strftime(elimit - t24h,'%Y%m%d')
522	retr_param_dict['date'] = '/'.join([startdate,
523	'to',
524	enddate])
525
526	# ******* start retrievement
527	self._start_retrievement(request, retr_param_dict)
528
529	retr_param_dict['date'] = \
530	datetime.strftime(elimit - t12h, '%Y%m%d')
531	retr_param_dict['time'] = '00'
532	retr_param_dict['target'] = \
533	self._mk_targetname(ftype, pk,
534	retr_param_dict['date'])
535
536	# ******* start retrievement
537	self._start_retrievement(request, retr_param_dict)
538
539	# -------------- non flux data ------------------------
540	else:
541	# ******* start retrievement
542	self._start_retrievement(request, retr_param_dict)
543
544	else: # basetime = 0
545	retr_param_dict['date'] = \
546	datetime.strftime(elimit - t24h, '%Y%m%d')
547
548	timesave = ''.join(retr_param_dict['time'])
549
550	if '/' in retr_param_dict['time']:
551	times = retr_param_dict['time'].split('/')
552	steps = retr_param_dict['step'].split('/')
553	while (pk != 'OG_OROLSM__SL' and
554	'acc' not in pk and
555	(int(times[0]) + int(steps[0])) <= 12):
556	times = times[1:]
557
558	if len(times) > 1:
559	retr_param_dict['time'] = '/'.join(times)
560	else:
561	retr_param_dict['time'] = times[0]
562
563	# ******* start retrievement
564	self._start_retrievement(request, retr_param_dict)
565
566	if (pk != 'OG_OROLSM__SL' and
567	int(retr_param_dict['step'].split('/')[0]) == 0 and
568	int(timesave.split('/')[0]) == 0):
569
570	retr_param_dict['date'] = \
571	datetime.strftime(elimit, '%Y%m%d')
572	retr_param_dict['time'] = '00'
573	retr_param_dict['step'] = '000'
574	retr_param_dict['target'] = \
575	self._mk_targetname(ftype, pk,
576	retr_param_dict['date'])
577
578	# ******* start retrievement
579	self._start_retrievement(request, retr_param_dict)
580
581	if request == 0 or request == 2:
582	print('MARS retrieve done ... ')
583	elif request == 1:
584	print('MARS request printed ...')
585
586	return
587
588
589	def write_namelist(self, c, filename):
590	'''
591	@Description:
592	Creates a namelist file in the temporary directory and writes
593	the following values to it: maxl, maxb, mlevel,
594	mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
595	momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
596
597	@Input:
598	self: instance of EcFlexpart
599	The current object of the class.
600
601	c: instance of class ControlFile
602	Contains all the parameters of CONTROL file and
603	command line.
604	For more information about format and content of the parameter
605	see documentation.
606
607	filename: string
608	Name of the namelist file.
609
610	@Return:
611	<nothing>
612	'''
613
614	self.inputdir = c.inputdir
615	area = np.asarray(self.area.split('/')).astype(float)
616	grid = np.asarray(self.grid.split('/')).astype(float)
617
618	if area[1] > area[3]:
619	area[1] -= 360
620	maxl = int((area[3] - area[1]) / grid[1]) + 1
621	maxb = int((area[0] - area[2]) / grid[0]) + 1
622
623	with open(self.inputdir + '/' + filename, 'w') as f:
624	f.write('&NAMGEN\n')
625	f.write(',\n '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
626	'mlevel = ' + str(self.level),
627	'mlevelist = ' + '"' + str(self.levelist)
628	+ '"',
629	'mnauf = ' + str(self.resol),
630	'metapar = ' + '77',
631	'rlo0 = ' + str(area[1]),
632	'rlo1 = ' + str(area[3]),
633	'rla0 = ' + str(area[2]),
634	'rla1 = ' + str(area[0]),
635	'momega = ' + str(c.omega),
636	'momegadiff = ' + str(c.omegadiff),
637	'mgauss = ' + str(c.gauss),
638	'msmooth = ' + str(c.smooth),
639	'meta = ' + str(c.eta),
640	'metadiff = ' + str(c.etadiff),
641	'mdpdeta = ' + str(c.dpdeta)]))
642
643	f.write('\n/\n')
644
645	return
646
647
648	def deacc_fluxes(self, inputfiles, c):
649	'''
650	@Description:
651	Goes through all flux fields in ordered time and de-accumulate
652	the fields. Afterwards the fields are disaggregated in time.
653	Different versions of disaggregation is provided for rainfall
654	data (darain, modified linear) and the surface fluxes and
655	stress data (dapoly, cubic polynomial).
656
657	@Input:
658	self: instance of EcFlexpart
659	The current object of the class.
660
661	inputfiles: instance of UioFiles
662	Contains a list of files.
663
664	c: instance of class ControlFile
665	Contains all the parameters of CONTROL file and
666	command line.
667	For more information about format and content of the parameter
668	see documentation.
669
670	@Return:
671	<nothing>
672	'''
673
674	table128 = init128(_config.PATH_GRIBTABLE)
675	pars = to_param_id(self.params['OG_acc_SL'][0], table128)
676	index_keys = ["date", "time", "step"]
677	indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
678	silent_remove(indexfile)
679	grib = GribTools(inputfiles.files)
680	# creates new index file
681	iid = grib.index(index_keys=index_keys, index_file=indexfile)
682
683	# read values of index keys
684	index_vals = []
685	for key in index_keys:
686	key_vals = grib_index_get(iid, key)
687	print(key_vals)
688	# have to sort the steps for disaggregation,
689	# therefore convert to int first
690	if key == 'step':
691	key_vals = [int(k) for k in key_vals]
692	key_vals.sort()
693	key_vals = [str(k) for k in key_vals]
694	index_vals.append(key_vals)
695	# index_vals looks for example like:
696	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
697	# index_vals[1]: ('0', '1200') ; time
698	# index_vals[2]: (3', '6', '9', '12') ; stepRange
699
700	valsdict = {}
701	svalsdict = {}
702	stepsdict = {}
703	for p in pars:
704	valsdict[str(p)] = []
705	svalsdict[str(p)] = []
706	stepsdict[str(p)] = []
707
708	print('maxstep: ', c.maxstep)
709
710	for prod in product(*index_vals):
711	# e.g. prod = ('20170505', '0', '12')
712	# ( date ,time, step)
713	# per date e.g. time = 0, 1200
714	# per time e.g. step = 3, 6, 9, 12
715	print('current prod: ', prod)
716	for i in range(len(index_keys)):
717	grib_index_select(iid, index_keys[i], prod[i])
718
719	# get first id from current product
720	gid = grib_new_from_index(iid)
721
722	# if there is data for this product combination
723	# prepare some date and time parameter before reading the data
724	if gid is not None:
725	cdate = grib_get(gid, 'date')
726	time = grib_get(gid, 'time')
727	step = grib_get(gid, 'step')
728	# date+time+step-2*dtime
729	# (since interpolated value valid for step-2*dtime)
730	sdate = datetime(year=cdate/10000,
731	month=(cdate % 10000)/100,
732	day=(cdate % 100),
733	hour=time/100)
734	fdate = sdate + timedelta(hours=step-2*int(c.dtime))
735	sdates = sdate + timedelta(hours=step)
736	elimit = None
737	else:
738	break
739
740	if c.maxstep > 12:
741	fnout = os.path.join(c.inputdir, 'flux' +
742	sdate.strftime('%Y%m%d') +
743	'.{:0>2}'.format(time/100) +
744	'.{:0>3}'.format(step-2*int(c.dtime)))
745	gnout = os.path.join(c.inputdir, 'flux' +
746	sdate.strftime('%Y%m%d') +
747	'.{:0>2}'.format(time/100) +
748	'.{:0>3}'.format(step-int(c.dtime)))
749	hnout = os.path.join(c.inputdir, 'flux' +
750	sdate.strftime('%Y%m%d') +
751	'.{:0>2}'.format(time/100) +
752	'.{:0>3}'.format(step))
753	else:
754	fnout = os.path.join(c.inputdir, 'flux' +
755	fdate.strftime('%Y%m%d%H'))
756	gnout = os.path.join(c.inputdir, 'flux' +
757	(fdate + timedelta(hours=int(c.dtime))).
758	strftime('%Y%m%d%H'))
759	hnout = os.path.join(c.inputdir, 'flux' +
760	sdates.strftime('%Y%m%d%H'))
761
762	print("outputfile = " + fnout)
763	f_handle = open(fnout, 'w')
764	g_handle = open(gnout, 'w')
765	h_handle = open(hnout, 'w')
766
767	# read message for message and store relevant data fields
768	# data keywords are stored in pars
769	while 1:
770	if gid is None:
771	break
772	cparamId = str(grib_get(gid, 'paramId'))
773	step = grib_get(gid, 'step')
774	atime = grib_get(gid, 'time')
775	ni = grib_get(gid, 'Ni')
776	nj = grib_get(gid, 'Nj')
777	if cparamId in valsdict.keys():
778	values = grib_get_values(gid)
779	vdp = valsdict[cparamId]
780	svdp = svalsdict[cparamId]
781	sd = stepsdict[cparamId]
782
783	if cparamId == '142' or cparamId == '143':
784	fak = 1. / 1000.
785	else:
786	fak = 3600.
787
788	values = (np.reshape(values, (nj, ni))).flatten() / fak
789	vdp.append(values[:]) # save the accumulated values
790	if step <= int(c.dtime):
791	svdp.append(values[:] / int(c.dtime))
792	else: # deaccumulate values
793	svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
794
795	print(cparamId, atime, step, len(values),
796	values[0], np.std(values))
797	# save the 1/3-hourly or specific values
798	# svdp.append(values[:])
799	sd.append(step)
800	# len(svdp) correspond to the time
801	if len(svdp) >= 3:
802	if len(svdp) > 3:
803	if cparamId == '142' or cparamId == '143':
804	values = disaggregation.darain(svdp)
805	else:
806	values = disaggregation.dapoly(svdp)
807
808	if not (step == c.maxstep and c.maxstep > 12 \
809	or sdates == elimit):
810	vdp.pop(0)
811	svdp.pop(0)
812	else:
813	if c.maxstep > 12:
814	values = svdp[1]
815	else:
816	values = svdp[0]
817
818	grib_set_values(gid, values)
819	if c.maxstep > 12:
820	grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
821	else:
822	grib_set(gid, 'step', 0)
823	grib_set(gid, 'time', fdate.hour*100)
824	grib_set(gid, 'date', fdate.year*10000 +
825	fdate.month*100+fdate.day)
826	grib_write(gid, f_handle)
827
828	if c.basetime:
829	elimit = datetime.strptime(c.end_date +
830	c.basetime, '%Y%m%d%H')
831	else:
832	elimit = sdate + timedelta(2*int(c.dtime))
833
834	# squeeze out information of last two steps contained
835	# in svdp
836	# if step+int(c.dtime) == c.maxstep and c.maxstep>12
837	# or sdates+timedelta(hours = int(c.dtime))
838	# >= elimit:
839	# Note that svdp[0] has not been popped in this case
840
841	if step == c.maxstep and c.maxstep > 12 or \
842	sdates == elimit:
843
844	values = svdp[3]
845	grib_set_values(gid, values)
846	grib_set(gid, 'step', 0)
847	truedatetime = fdate + timedelta(hours=
848	2*int(c.dtime))
849	grib_set(gid, 'time', truedatetime.hour * 100)
850	grib_set(gid, 'date', truedatetime.year * 10000 +
851	truedatetime.month * 100 +
852	truedatetime.day)
853	grib_write(gid, h_handle)
854
855	#values = (svdp[1]+svdp[2])/2.
856	if cparamId == '142' or cparamId == '143':
857	values = disaggregation.darain(list(reversed(svdp)))
858	else:
859	values = disaggregation.dapoly(list(reversed(svdp)))
860
861	grib_set(gid, 'step', 0)
862	truedatetime = fdate + timedelta(hours=int(c.dtime))
863	grib_set(gid, 'time', truedatetime.hour * 100)
864	grib_set(gid, 'date', truedatetime.year * 10000 +
865	truedatetime.month * 100 +
866	truedatetime.day)
867	grib_set_values(gid, values)
868	grib_write(gid, g_handle)
869
870	grib_release(gid)
871
872	gid = grib_new_from_index(iid)
873
874	f_handle.close()
875	g_handle.close()
876	h_handle.close()
877
878	grib_index_release(iid)
879
880	return
881
882
883	def create(self, inputfiles, c):
884	'''
885	@Description:
886	This method is based on the ECMWF example index.py
887	https://software.ecmwf.int/wiki/display/GRIB/index.py
888
889	An index file will be created which depends on the combination
890	of "date", "time" and "stepRange" values. This is used to iterate
891	over all messages in each grib file which were passed through the
892	parameter "inputfiles" to seperate specific parameters into fort.*
893	files. Afterwards the FORTRAN program is called to convert
894	the data fields all to the same grid and put them in one file
895	per unique time step (combination of "date", "time" and
896	"stepRange").
897
898	@Input:
899	self: instance of EcFlexpart
900	The current object of the class.
901
902	inputfiles: instance of UioFiles
903	Contains a list of files.
904
905	c: instance of class ControlFile
906	Contains all the parameters of CONTROL file and
907	command line.
908	For more information about format and content of the parameter
909	see documentation.
910
911	@Return:
912	<nothing>
913	'''
914
915	table128 = init128(_config.PATH_GRIBTABLE)
916	wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
917	stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
918	table128)
919
920	index_keys = ["date", "time", "step"]
921	indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
922	silent_remove(indexfile)
923	grib = GribTools(inputfiles.files)
924	# creates new index file
925	iid = grib.index(index_keys=index_keys, index_file=indexfile)
926
927	# read values of index keys
928	index_vals = []
929	for key in index_keys:
930	index_vals.append(grib_index_get(iid, key))
931	print(index_vals[-1])
932	# index_vals looks for example like:
933	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
934	# index_vals[1]: ('0', '1200', '1800', '600') ; time
935	# index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
936
937	fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
938	'17':None, '19':None, '21':None, '22':None, '20':None}
939
940	for prod in product(*index_vals):
941	# e.g. prod = ('20170505', '0', '12')
942	# ( date ,time, step)
943	# per date e.g. time = 0, 1200
944	# per time e.g. step = 3, 6, 9, 12
945	# flag for Fortran program and file merging
946	convertFlag = False
947	print('current prod: ', prod)
948	# e.g. prod = ('20170505', '0', '12')
949	# ( date ,time, step)
950	# per date e.g. time = 0, 600, 1200, 1800
951	# per time e.g. step = 0, 3, 6, 9, 12
952	for i in range(len(index_keys)):
953	grib_index_select(iid, index_keys[i], prod[i])
954
955	# get first id from current product
956	gid = grib_new_from_index(iid)
957
958	# if there is data for this product combination
959	# prepare some date and time parameter before reading the data
960	if gid is not None:
961	# Combine all temporary data files into final grib file if
962	# gid is at least one time not None. Therefore set convertFlag
963	# to save information. The Fortran program is also
964	# only executed if convertFlag is True
965	convertFlag = True
966	# remove old fort.* files and open new ones
967	# they are just valid for a single product
968	for k, f in fdict.iteritems():
969	fortfile = os.path.join(c.inputdir, 'fort.' + k)
970	silent_remove(fortfile)
971	fdict[k] = open(fortfile, 'w')
972
973	cdate = str(grib_get(gid, 'date'))
974	time = grib_get(gid, 'time')
975	step = grib_get(gid, 'step')
976	# create correct timestamp from the three time informations
977	# date, time, step
978	timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100),
979	'%Y%m%d%H')
980	timestamp += timedelta(hours=int(step))
981	cdateH = datetime.strftime(timestamp, '%Y%m%d%H')
982
983	if c.basetime:
984	slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
985	bt = '23'
986	if c.basetime == '00':
987	bt = '00'
988	slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
989	- timedelta(hours=12-int(c.dtime))
990	if c.basetime == '12':
991	bt = '12'
992	slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
993	- timedelta(hours=12-int(c.dtime))
994
995	elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')
996
997	if timestamp < slimit or timestamp > elimit:
998	continue
999	else:
1000	pass
1001
1002	try:
1003	if c.wrf:
1004	if 'olddate' not in locals() or cdate != olddate:
1005	fwrf = open(os.path.join(c.outputdir,
1006	'WRF' + cdate + '.{:0>2}'.format(time) +
1007	'.000.grb2'), 'w')
1008	olddate = cdate[:]
1009	except AttributeError:
1010	pass
1011
1012	# helper variable to remember which fields were already used.
1013	savedfields = []
1014	while 1:
1015	if gid is None:
1016	break
1017	paramId = grib_get(gid, 'paramId')
1018	gridtype = grib_get(gid, 'gridType')
1019	levtype = grib_get(gid, 'typeOfLevel')
1020	if paramId == 133 and gridtype == 'reduced_gg':
1021	# Specific humidity (Q.grb) is used as a template only
1022	# so we need the first we "meet"
1023	with open(os.path.join(c.inputdir, 'fort.18'), 'w') as fout:
1024	grib_write(gid, fout)
1025	elif paramId == 131 or paramId == 132:
1026	grib_write(gid, fdict['10'])
1027	elif paramId == 130:
1028	grib_write(gid, fdict['11'])
1029	elif paramId == 133 and gridtype != 'reduced_gg':
1030	grib_write(gid, fdict['17'])
1031	elif paramId == 152:
1032	grib_write(gid, fdict['12'])
1033	elif paramId == 155 and gridtype == 'sh':
1034	grib_write(gid, fdict['13'])
1035	elif paramId in [129, 138, 155] and levtype == 'hybrid' \
1036	and c.wrf:
1037	pass
1038	elif paramId == 246 or paramId == 247:
1039	# cloud liquid water and ice
1040	if paramId == 246:
1041	clwc = grib_get_values(gid)
1042	else:
1043	clwc += grib_get_values(gid)
1044	grib_set_values(gid, clwc)
1045	grib_set(gid, 'paramId', 201031)
1046	grib_write(gid, fdict['22'])
1047	elif paramId == 135:
1048	grib_write(gid, fdict['19'])
1049	elif paramId == 77:
1050	grib_write(gid, fdict['21'])
1051	else:
1052	if paramId not in savedfields:
1053	grib_write(gid, fdict['16'])
1054	savedfields.append(paramId)
1055	else:
1056	print('duplicate ' + str(paramId) + ' not written')
1057
1058	try:
1059	if c.wrf:
1060	# model layer
1061	if levtype == 'hybrid' and \
1062	paramId in [129, 130, 131, 132, 133, 138, 155]:
1063	grib_write(gid, fwrf)
1064	# sfc layer
1065	elif paramId in wrfpars:
1066	grib_write(gid, fwrf)
1067	except AttributeError:
1068	pass
1069
1070	grib_release(gid)
1071	gid = grib_new_from_index(iid)
1072
1073	for f in fdict.values():
1074	f.close()
1075
1076	# call for Fortran program if flag is True
1077	if convertFlag:
1078	pwd = os.getcwd()
1079	os.chdir(c.inputdir)
1080	if os.stat('fort.21').st_size == 0 and c.eta:
1081	print('Parameter 77 (etadot) is missing, most likely it is \
1082	not available for this type or date/time\n')
1083	print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1084	my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1085	is set to 1 in CONTROL file')
1086
1087	# create the corresponding output file fort.15
1088	# (generated by Fortran program) + fort.16 (paramId 167 and 168)
1089	p = subprocess.check_call([os.path.join(
1090	c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
1091	os.chdir(pwd)
1092
1093	# create final output filename, e.g. EN13040500 (ENYYMMDDHH)
1094	fnout = os.path.join(c.inputdir, c.prefix)
1095	if c.maxstep > 12:
1096	suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
1097	'.{:0>3}'.format(step)
1098	else:
1099	suffix = cdateH[2:10]
1100	fnout += suffix
1101	print("outputfile = " + fnout)
1102	self.outputfilelist.append(fnout) # needed for final processing
1103
1104	# create outputfile and copy all data from intermediate files
1105	# to the outputfile (final GRIB files)
1106	orolsm = os.path.basename(glob.glob(
1107	c.inputdir + '/OG_OROLSM__SL..' + c.ppid + '')[0])
1108	fluxfile = 'flux' + cdate[0:2] + suffix
1109	if not c.cwc:
1110	flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1111	else:
1112	flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1113
1114	with open(fnout, 'wb') as fout:
1115	for f in flist:
1116	shutil.copyfileobj(open(os.path.join(c.inputdir, f),
1117	'rb'), fout)
1118
1119	if c.omega:
1120	with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1121	shutil.copyfileobj(
1122	open(os.path.join(c.inputdir, 'fort.25'),
1123	'rb'), fout)
1124	else:
1125	pass
1126
1127	if c.wrf:
1128	fwrf.close()
1129
1130	grib_index_release(iid)
1131
1132	return
1133
1134
1135	def process_output(self, c):
1136	'''
1137	@Description:
1138	The grib files are postprocessed depending on the selection in
1139	CONTROL file. The resulting files are moved to the output
1140	directory if its not equal to the input directory.
1141	The following modifications might be done if
1142	properly switched in CONTROL file:
1143	GRIB2 - Conversion to GRIB2
1144	ECTRANS - Transfer of files to gateway server
1145	ECSTORAGE - Storage at ECMWF server
1146
1147	@Input:
1148	self: instance of EcFlexpart
1149	The current object of the class.
1150
1151	c: instance of class ControlFile
1152	Contains all the parameters of CONTROL file and
1153	command line.
1154	For more information about format and content of the parameter
1155	see documentation.
1156
1157	@Return:
1158	<nothing>
1159
1160	'''
1161
1162	print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
1163
1164	if not c.ecapi:
1165	print('ecstorage: {}\n ecfsdir: {}\n'.
1166	format(c.ecstorage, c.ecfsdir))
1167	print('ectrans: {}\n gateway: {}\n destination: {}\n '
1168	.format(c.ectrans, c.gateway, c.destination))
1169
1170	print('Output filelist: \n')
1171	print(self.outputfilelist)
1172
1173	if c.format.lower() == 'grib2':
1174	for ofile in self.outputfilelist:
1175	p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
1176	productDefinitionTemplateNumber=8',
1177	ofile, ofile + '_2'])
1178	p = subprocess.check_call(['mv', ofile + '_2', ofile])
1179
1180	if c.ectrans and not c.ecapi:
1181	for ofile in self.outputfilelist:
1182	p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
1183	c.gateway, '-remote', c.destination,
1184	'-source', ofile])
1185
1186	if c.ecstorage and not c.ecapi:
1187	for ofile in self.outputfilelist:
1188	p = subprocess.check_call(['ecp', '-o', ofile,
1189	os.path.expandvars(c.ecfsdir)])
1190
1191	if c.outputdir != c.inputdir:
1192	for ofile in self.outputfilelist:
1193	p = subprocess.check_call(['mv', ofile, c.outputdir])
1194
1195	return
1196
1197
1198	def prepare_fp_files(self, c):
1199	'''
1200	@Description:
1201	Conversion of GRIB files to FLEXPART binary format.
1202
1203	@Input:
1204	c: instance of class ControlFile
1205	Contains all the parameters of CONTROL file and
1206	command line.
1207	For more information about format and content of the parameter
1208	see documentation.
1209
1210
1211	@Return:
1212	<nothing>
1213	'''
1214	# generate AVAILABLE file
1215	# Example of AVAILABLE file data:
1216	# 20131107 000000 EN13110700 ON DISC
1217	clist = []
1218	for ofile in self.outputfilelist:
1219	fname = ofile.split('/')
1220	if '.' in fname[-1]:
1221	l = fname[-1].split('.')
1222	timestamp = datetime.strptime(l[0][-6:] + l[1],
1223	'%y%m%d%H')
1224	timestamp += timedelta(hours=int(l[2]))
1225	cdate = datetime.strftime(timestamp, '%Y%m%d')
1226	chms = datetime.strftime(timestamp, '%H%M%S')
1227	else:
1228	cdate = '20' + fname[-1][-8:-2]
1229	chms = fname[-1][-2:] + '0000'
1230	clist.append(cdate + ' ' + chms + ' '*6 +
1231	fname[-1] + ' '*14 + 'ON DISC')
1232	clist.sort()
1233	with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
1234	f.write('\n'.join(clist) + '\n')
1235
1236	# generate pathnames file
1237	pwd = os.path.abspath(c.outputdir)
1238	with open(pwd + '/pathnames', 'w') as f:
1239	f.write(pwd + '/Options/\n')
1240	f.write(pwd + '/\n')
1241	f.write(pwd + '/\n')
1242	f.write(pwd + '/AVAILABLE\n')
1243	f.write(' = == = == = == = == = == == = \n')
1244
1245	# create Options dir if necessary
1246	if not os.path.exists(pwd + '/Options'):
1247	os.makedirs(pwd+'/Options')
1248
1249	# read template COMMAND file
1250	with open(os.path.expandvars(os.path.expanduser(
1251	c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
1252	lflist = f.read().split('\n')
1253
1254	# find index of list where to put in the
1255	# date and time information
1256	# usually after the LDIRECT parameter
1257	i = 0
1258	for l in lflist:
1259	if 'LDIRECT' in l.upper():
1260	break
1261	i += 1
1262
1263	# insert the date and time information of run start and end
1264	# into the list of lines of COMMAND file
1265	lflist = lflist[:i+1] + \
1266	[clist[0][:16], clist[-1][:16]] + \
1267	lflist[i+3:]
1268
1269	# write the new COMMAND file
1270	with open(pwd + '/Options/COMMAND', 'w') as g:
1271	g.write('\n'.join(lflist) + '\n')
1272
1273	# change to outputdir and start the grib2flexpart run
1274	# afterwards switch back to the working dir
1275	os.chdir(c.outputdir)
1276	p = subprocess.check_call([
1277	os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
1278	+ '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
1279	os.chdir(pwd)
1280
1281	return

Note: See TracBrowser for help on using the repository browser.

Download in other formats: