Context Navigation

source: flex_extract.git/source/python/classes/EcFlexpart.py @ 3f36e42

ctbtodev

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

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

Download in other formats: