Context Navigation

source: flex_extract.git/source/python/classes/EcFlexpart.py @ 274f9ef

ctbtodev

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

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

Download in other formats: