Coverage for webgeodyn/inout/pygeodyn_hdf5.py: 0%

1#-*- coding: utf-8 -*- 

2import glob 

3import h5py 

4import os 

5import numpy as np 

6from .default import giveMeasureTypeUnits 

7 

8 

9def load(dataDirectory, dataModel, keepRealisations, state_type='analysed'): 

10 """ Loading function for pygeodyn files of hdf5 format. Also adds the data to the dataModel. 

11 

12 :param dataDirectory: Location of the pygeodyn files 

13 :type dataDirectory: os.path 

14 :param dataModel: Model in which to add the loaded measures 

15 :type dataModel: Model 

16 :param keepRealisations: If True, all realisations are kept in the data. Else, the data is averaged over the realisations 

17 :type keepRealisations: bool 

18 :return: 0 if everything went well, -1 otherwise 

19 :rtype: int 

20 :param state_type: Either forecast, computed or analysed depending on the type of states needed 

21 """ 

22 firstpoint = 3 

23 assert state_type in ('computed', 'analysed', 'analysis', 'forecast') 

24 measures_to_load = ['MF', 'SV', 'ER', 'U', 'S'] 

25 

26 hdf5_files = glob.glob(os.path.join(dataDirectory, '*.hdf5')) 

27 if len(hdf5_files) == 0: 

28 raise IOError('No hdf5 file was found in {} !'.format(dataDirectory)) 

29 # Assuming that the file to read is the first one 

30 hdf_filename = hdf5_files[0] 

31 print('Reading:', hdf_filename, ' state_type:', state_type) 

32 

33 with h5py.File(hdf_filename) as hdf_file: 

34 computed_data = hdf_file[state_type] 

35 

36 times = np.array(computed_data['times'], dtype="float64")[firstpoint:] 

37 

38 for measureName, measureData in computed_data.items(): 

39 measureData = np.float64(measureData) # Convert to float64 

40 if measureName not in measures_to_load: 

41 continue 

42 else: 

43 measureType, units = giveMeasureTypeUnits(measureName) 

44 

45 # Move realisation axis to last place to have data of form [ntimes, ncoef, nreal] (originally [nreal, ntimes, ncoef]) 

46 # Remove firstpoints 

47 formattedData = np.moveaxis(measureData, 0, -1)[firstpoint:] 

48 

49 if measureName == 'MF': 

50 lmax = hdf_file.attrs['Lb'] 

51 elif measureName == 'U': 

52 lmax = hdf_file.attrs['Lu'] 

53 elif measureName == 'S': 

54 lmax = hdf_file.attrs['Lu'] 

55 else: 

56 lmax = hdf_file.attrs['Lsv'] 

57 

58 if keepRealisations: 

59 dataModel.addMeasure(measureName, measureType, lmax, units, 

60 formattedData, times=times) 

61 else: 

62 meanData = formattedData.mean(axis=2) 

63 rmsData = formattedData.std(axis=2) 

64 dataModel.addMeasure(measureName, measureType, lmax, units, 

65 meanData, rmsData, times) 

66 

67 # Returns 0 if everything went well 

68 return 0