Context Navigation

source: flex_extract.git/source/python/classes/EcFlexpart.py @ 70fee58

ctbtodev

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

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

Download in other formats: