Context Navigation

source: flex_extract.git/source/python/classes/EcFlexpart.py @ 96e1533

ctbtodev

Last change on this file since 96e1533 was 96e1533, checked in by Anne Philipp <anne.philipp@…>, 5 years ago
redefined test data dir and completed unittests for tools module
Property mode set to `100644`
File size: 52.3 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	# initialize dictionaries
758	valsdict = {}
759	svalsdict = {}
760	for p in pars:
761	valsdict[str(p)] = []
762	svalsdict[str(p)] = []
763
764	# "product" genereates each possible combination between the
765	# values of the index keys
766	for prod in product(*index_vals):
767	# e.g. prod = ('20170505', '0', '12')
768	# ( date ,time, step)
769
770	print('current product: ', prod)
771
772	for i in range(len(index_keys)):
773	grib_index_select(iid, index_keys[i], prod[i])
774
775	# get first id from current product
776	gid = grib_new_from_index(iid)
777
778	# if there is no data for this specific time combination / product
779	# skip the rest of the for loop and start with next timestep/product
780	if not gid:
781	continue
782
783	# create correct timestamp from the three time informations
784	cdate = str(grib_get(gid, 'date'))
785	ctime = '{:0>2}'.format(grib_get(gid, 'time')/100)
786	cstep = '{:0>3}'.format(grib_get(gid, 'step'))
787	t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
788	t_dt = t_date + timedelta(hours=int(cstep))
789	t_m1dt = t_date + timedelta(hours=int(cstep)-int(c.dtime))
790	t_m2dt = t_date + timedelta(hours=int(cstep)-2*int(c.dtime))
791	t_enddate = None
792
793	if c.maxstep > 12:
794	fnout = os.path.join(c.inputdir, 'flux' +
795	t_date.strftime('%Y%m%d.%H') +
796	'.{:0>3}'.format(step-2*int(c.dtime)))
797	gnout = os.path.join(c.inputdir, 'flux' +
798	t_date.strftime('%Y%m%d.%H') +
799	'.{:0>3}'.format(step-int(c.dtime)))
800	hnout = os.path.join(c.inputdir, 'flux' +
801	t_date.strftime('%Y%m%d.%H') +
802	'.{:0>3}'.format(step))
803	else:
804	fnout = os.path.join(c.inputdir, 'flux' +
805	t_m2dt.strftime('%Y%m%d%H'))
806	gnout = os.path.join(c.inputdir, 'flux' +
807	t_m1dt.strftime('%Y%m%d%H'))
808	hnout = os.path.join(c.inputdir, 'flux' +
809	t_dt.strftime('%Y%m%d%H'))
810
811	print("outputfile = " + fnout)
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(grib_get(gid, 'paramId'))
819	step = grib_get(gid, 'step')
820	time = grib_get(gid, 'time')
821	ni = grib_get(gid, 'Ni')
822	nj = grib_get(gid, 'Nj')
823	if cparamId in valsdict.keys():
824	values = grib_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) correspond 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	grib_set_values(gid, values)
863
864	if c.maxstep > 12:
865	grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
866	else:
867	grib_set(gid, 'step', 0)
868	grib_set(gid, 'time', t_m2dt.hour*100)
869	grib_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d')))
870
871	with open(fnout, 'w') as f_handle:
872	grib_write(gid, f_handle)
873
874	if c.basetime:
875	t_enddate = datetime.strptime(c.end_date +
876	c.basetime,
877	'%Y%m%d%H')
878	else:
879	t_enddate = t_date + timedelta(2*int(c.dtime))
880
881	# squeeze out information of last two steps contained
882	# in svdp
883	# if step+int(c.dtime) == c.maxstep and c.maxstep>12
884	# or t_dt+timedelta(hours = int(c.dtime))
885	# >= t_enddate:
886	# Note that svdp[0] has not been popped in this case
887
888	if step == c.maxstep and c.maxstep > 12 or \
889	t_dt == t_enddate:
890
891	values = svdp[3]
892	grib_set_values(gid, values)
893	grib_set(gid, 'step', 0)
894	truedatetime = t_m2dt + timedelta(hours=
895	2*int(c.dtime))
896	grib_set(gid, 'time', truedatetime.hour * 100)
897	grib_set(gid, 'date', truedatetime.year * 10000 +
898	truedatetime.month * 100 +
899	truedatetime.day)
900	with open(hnout, 'w') as h_handle:
901	grib_write(gid, h_handle)
902
903	#values = (svdp[1]+svdp[2])/2.
904	if cparamId == '142' or cparamId == '143':
905	values = disaggregation.darain(list(reversed(svdp)))
906	else:
907	values = disaggregation.dapoly(list(reversed(svdp)))
908
909	grib_set(gid, 'step', 0)
910	truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
911	grib_set(gid, 'time', truedatetime.hour * 100)
912	grib_set(gid, 'date', truedatetime.year * 10000 +
913	truedatetime.month * 100 +
914	truedatetime.day)
915	grib_set_values(gid, values)
916	with open(gnout, 'w') as g_handle:
917	grib_write(gid, g_handle)
918
919	grib_release(gid)
920
921	gid = grib_new_from_index(iid)
922
923	grib_index_release(iid)
924
925	return
926
927
928	def create(self, inputfiles, c):
929	'''An index file will be created which depends on the combination
930	of "date", "time" and "stepRange" values. This is used to iterate
931	over all messages in each grib file which were passed through the
932	parameter "inputfiles" to seperate specific parameters into fort.*
933	files. Afterwards the FORTRAN program is called to convert
934	the data fields all to the same grid and put them in one file
935	per unique time step (combination of "date", "time" and
936	"stepRange").
937
938	Note
939	----
940	This method is based on the ECMWF example index.py
941	https://software.ecmwf.int/wiki/display/GRIB/index.py
942
943	Parameters
944	----------
945	inputfiles : :obj:`UioFiles`
946	Contains a list of files.
947
948	c : :obj:`ControlFile`
949	Contains all the parameters of CONTROL file and
950	command line.
951
952	Return
953	------
954
955	'''
956
957	if c.wrf:
958	table128 = init128(_config.PATH_GRIBTABLE)
959	wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
960	stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
961	table128)
962
963	# these numbers are indices for the temporary files "fort.xx"
964	# which are used to seperate the grib fields to,
965	# for the Fortran program input
966	# 10: U,V \| 11: T \| 12: lnsp \| 13: D \| 16: sfc fields
967	# 17: Q \| 18: Q , gaussian\| 19: w \| 21: etadot \| 22: clwc+ciwc
968	fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
969	'17':None, '18':None, '19':None, '21':None, '22':None}
970
971	iid = None
972	index_vals = None
973
974	# get the values of the keys which are used for distinct access
975	# of grib messages via product
976	index_keys = ["date", "time", "step"]
977	iid, index_vals = self._mk_index_values(c.inputdir,
978	inputfiles,
979	index_keys)
980	# index_vals looks like e.g.:
981	# index_vals[0]: ('20171106', '20171107', '20171108') ; date
982	# index_vals[1]: ('0', '1200', '1800', '600') ; time
983	# index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
984
985	# "product" genereates each possible combination between the
986	# values of the index keys
987	for prod in product(*index_vals):
988	# e.g. prod = ('20170505', '0', '12')
989	# ( date ,time, step)
990
991	print('current product: ', prod)
992
993	for i in range(len(index_keys)):
994	grib_index_select(iid, index_keys[i], prod[i])
995
996	# get first id from current product
997	gid = grib_new_from_index(iid)
998
999	# if there is no data for this specific time combination / product
1000	# skip the rest of the for loop and start with next timestep/product
1001	if not gid:
1002	continue
1003
1004	# remove old fort.* files and open new ones
1005	# they are just valid for a single product
1006	for k, f in fdict.iteritems():
1007	fortfile = os.path.join(c.inputdir, 'fort.' + k)
1008	silent_remove(fortfile)
1009	fdict[k] = open(fortfile, 'w')
1010
1011	# create correct timestamp from the three time informations
1012	cdate = str(grib_get(gid, 'date'))
1013	ctime = '{:0>2}'.format(grib_get(gid, 'time')/100)
1014	cstep = '{:0>3}'.format(grib_get(gid, 'step'))
1015	timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H')
1016	timestamp += timedelta(hours=int(cstep))
1017	cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H')
1018
1019	# if the timestamp is out of basetime start/end date period,
1020	# skip this specific product
1021	if c.basetime:
1022	start_time = datetime.strptime(c.end_date + c.basetime,
1023	'%Y%m%d%H') - time_delta
1024	end_time = datetime.strptime(c.end_date + c.basetime,
1025	'%Y%m%d%H')
1026	if timestamp < start_time or timestamp > end_time:
1027	continue
1028
1029	if c.wrf:
1030	if 'olddate' not in locals() or cdate != olddate:
1031	fwrf = open(os.path.join(c.outputdir,
1032	'WRF' + cdate + '.' + ctime + '.000.grb2'), 'w')
1033	olddate = cdate[:]
1034
1035	# savedfields remembers which fields were already used.
1036	savedfields = []
1037	# sum of cloud liquid and ice water content
1038	scwc = None
1039	while 1:
1040	if not gid:
1041	break
1042	paramId = grib_get(gid, 'paramId')
1043	gridtype = grib_get(gid, 'gridType')
1044	levtype = grib_get(gid, 'typeOfLevel')
1045	if paramId == 77: # ETADOT
1046	grib_write(gid, fdict['21'])
1047	elif paramId == 130: # T
1048	grib_write(gid, fdict['11'])
1049	elif paramId == 131 or paramId == 132: # U, V wind component
1050	grib_write(gid, fdict['10'])
1051	elif paramId == 133 and gridtype != 'reduced_gg': # Q
1052	grib_write(gid, fdict['17'])
1053	elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian
1054	grib_write(gid, fdict['18'])
1055	elif paramId == 135: # W
1056	grib_write(gid, fdict['19'])
1057	elif paramId == 152: # LNSP
1058	grib_write(gid, fdict['12'])
1059	elif paramId == 155 and gridtype == 'sh': # D
1060	grib_write(gid, fdict['13'])
1061	elif paramId == 246 or paramId == 247: # CLWC, CIWC
1062	# sum cloud liquid water and ice
1063	if scwc is None:
1064	scwc = grib_get_values(gid)
1065	else:
1066	scwc += grib_get_values(gid)
1067	grib_set_values(gid, scwc)
1068	grib_set(gid, 'paramId', 201031)
1069	grib_write(gid, fdict['22'])
1070	elif c.wrf and paramId in [129, 138, 155] and \
1071	levtype == 'hybrid': # Z, VO, D
1072	# do not do anything right now
1073	# these are specific parameter for WRF
1074	pass
1075	else:
1076	if paramId not in savedfields:
1077	# SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR
1078	# and all ADDPAR parameter
1079	grib_write(gid, fdict['16'])
1080	savedfields.append(paramId)
1081	else:
1082	print('duplicate ' + str(paramId) + ' not written')
1083
1084	try:
1085	if c.wrf:
1086	# model layer
1087	if levtype == 'hybrid' and \
1088	paramId in [129, 130, 131, 132, 133, 138, 155]:
1089	grib_write(gid, fwrf)
1090	# sfc layer
1091	elif paramId in wrfpars:
1092	grib_write(gid, fwrf)
1093	except AttributeError:
1094	pass
1095
1096	grib_release(gid)
1097	gid = grib_new_from_index(iid)
1098
1099	for f in fdict.values():
1100	f.close()
1101
1102	# call for Fortran program to convert e.g. reduced_gg grids to
1103	# regular_ll and calculate detadot/dp
1104	pwd = os.getcwd()
1105	os.chdir(c.inputdir)
1106	if os.stat('fort.21').st_size == 0 and c.eta:
1107	print('Parameter 77 (etadot) is missing, most likely it is \
1108	not available for this type or date/time\n')
1109	print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
1110	my_error(c.mailfail, 'fort.21 is empty while parameter eta \
1111	is set to 1 in CONTROL file')
1112
1113	# Fortran program creates file fort.15 (with u,v,etadot,t,sp,q)
1114	p = subprocess.check_call([os.path.join(
1115	c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
1116	os.chdir(pwd)
1117
1118	# create name of final output file, e.g. EN13040500 (ENYYMMDDHH)
1119	if c.maxstep > 12:
1120	suffix = cdate[2:8] + '.' + ctime + '.' + cstep
1121	else:
1122	suffix = cdate_hour[2:10]
1123	fnout = os.path.join(c.inputdir, c.prefix + suffix)
1124	print("outputfile = " + fnout)
1125	# collect for final processing
1126	self.outputfilelist.append(os.path.basename(fnout))
1127
1128	# create outputfile and copy all data from intermediate files
1129	# to the outputfile (final GRIB input files for FLEXPART)
1130	orolsm = os.path.basename(glob.glob(c.inputdir +
1131	'/OG_OROLSM__SL..' + c.ppid + '')[0])
1132	fluxfile = 'flux' + cdate[0:2] + suffix
1133	if not c.cwc:
1134	flist = ['fort.15', fluxfile, 'fort.16', orolsm]
1135	else:
1136	flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
1137
1138	with open(fnout, 'wb') as fout:
1139	for f in flist:
1140	shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'),
1141	fout)
1142
1143	if c.omega:
1144	with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
1145	shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'),
1146	'rb'), fout)
1147
1148	if c.wrf:
1149	fwrf.close()
1150
1151	grib_index_release(iid)
1152
1153	return
1154
1155
1156	def process_output(self, c):
1157	'''The grib files are postprocessed depending on the selection in
1158	CONTROL file. The resulting files are moved to the output
1159	directory if its not equal to the input directory.
1160	The following modifications might be done if
1161	properly switched in CONTROL file:
1162	GRIB2 - Conversion to GRIB2
1163	ECTRANS - Transfer of files to gateway server
1164	ECSTORAGE - Storage at ECMWF server
1165
1166	Parameters
1167	----------
1168	c : :obj:`ControlFile`
1169	Contains all the parameters of CONTROL file and
1170	command line.
1171
1172	Return
1173	------
1174
1175	'''
1176
1177	print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
1178
1179	if not c.ecapi:
1180	print('ecstorage: {}\n ecfsdir: {}\n'.
1181	format(c.ecstorage, c.ecfsdir))
1182	print('ectrans: {}\n gateway: {}\n destination: {}\n '
1183	.format(c.ectrans, c.gateway, c.destination))
1184
1185	print('Output filelist: ')
1186	print(self.outputfilelist)
1187
1188	for ofile in self.outputfilelist:
1189	ofile = os.path.join(self.inputdir, ofile)
1190
1191	if c.format.lower() == 'grib2':
1192	p = subprocess.check_call(['grib_set', '-s', 'edition=2,',
1193	'productDefinitionTemplateNumber=8',
1194	ofile, ofile + '_2'])
1195	p = subprocess.check_call(['mv', ofile + '_2', ofile])
1196
1197	if c.ectrans and not c.ecapi:
1198	p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
1199	c.gateway, '-remote', c.destination,
1200	'-source', ofile])
1201
1202	if c.ecstorage and not c.ecapi:
1203	p = subprocess.check_call(['ecp', '-o', ofile,
1204	os.path.expandvars(c.ecfsdir)])
1205
1206	if c.outputdir != c.inputdir:
1207	p = subprocess.check_call(['mv',
1208	os.path.join(c.inputdir, ofile),
1209	c.outputdir])
1210
1211	return
1212
1213
1214	def prepare_fp_files(self, c):
1215	'''Conversion of GRIB files to FLEXPART binary format.
1216
1217	Parameters
1218	----------
1219	c : :obj:`ControlFile`
1220	Contains all the parameters of CONTROL file and
1221	command line.
1222
1223	Return
1224	------
1225
1226	'''
1227	# generate AVAILABLE file
1228	# Example of AVAILABLE file data:
1229	# 20131107 000000 EN13110700 ON DISC
1230	clist = []
1231	for ofile in self.outputfilelist:
1232	fname = ofile.split('/')
1233	if '.' in fname[-1]:
1234	l = fname[-1].split('.')
1235	timestamp = datetime.strptime(l[0][-6:] + l[1],
1236	'%y%m%d%H')
1237	timestamp += timedelta(hours=int(l[2]))
1238	cdate = datetime.strftime(timestamp, '%Y%m%d')
1239	chms = datetime.strftime(timestamp, '%H%M%S')
1240	else:
1241	cdate = '20' + fname[-1][-8:-2]
1242	chms = fname[-1][-2:] + '0000'
1243	clist.append(cdate + ' ' + chms + ' '*6 +
1244	fname[-1] + ' '*14 + 'ON DISC')
1245	clist.sort()
1246	with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
1247	f.write('\n'.join(clist) + '\n')
1248
1249	# generate pathnames file
1250	pwd = os.path.abspath(c.outputdir)
1251	with open(pwd + '/pathnames', 'w') as f:
1252	f.write(pwd + '/Options/\n')
1253	f.write(pwd + '/\n')
1254	f.write(pwd + '/\n')
1255	f.write(pwd + '/AVAILABLE\n')
1256	f.write(' = == = == = == = == = == == = \n')
1257
1258	# create Options dir if necessary
1259	if not os.path.exists(pwd + '/Options'):
1260	make_dir(pwd+'/Options')
1261
1262	# read template COMMAND file
1263	with open(os.path.expandvars(os.path.expanduser(
1264	c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
1265	lflist = f.read().split('\n')
1266
1267	# find index of list where to put in the
1268	# date and time information
1269	# usually after the LDIRECT parameter
1270	i = 0
1271	for l in lflist:
1272	if 'LDIRECT' in l.upper():
1273	break
1274	i += 1
1275
1276	# insert the date and time information of run start and end
1277	# into the list of lines of COMMAND file
1278	lflist = lflist[:i+1] + \
1279	[clist[0][:16], clist[-1][:16]] + \
1280	lflist[i+3:]
1281
1282	# write the new COMMAND file
1283	with open(pwd + '/Options/COMMAND', 'w') as g:
1284	g.write('\n'.join(lflist) + '\n')
1285
1286	# change to outputdir and start the grib2flexpart run
1287	# afterwards switch back to the working dir
1288	os.chdir(c.outputdir)
1289	p = subprocess.check_call([
1290	os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
1291	+ '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
1292	os.chdir(pwd)
1293
1294	return

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

Download in other formats: