Coverage for webgeodyn/inout/pygeodyn_hdf5_forecasts.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): 

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 realisationsr 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 """ 

21 firstpoint = 0 # Do not cut any point 

22 

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

24 if len(hdf5_files) == 0: 

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

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

27 hdf_filename = hdf5_files[0] 

28 print('Reading:', hdf_filename) 

29 

30 with h5py.File(hdf_filename) as hdf_file: 

31 computed_data = hdf_file['computed'] 

32 

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

34 try: 

35 times_analysis = np.array(hdf_file['analysed']['times'], dtype="float64")[firstpoint:] 

36 except KeyError: 

37 times_analysis = np.array(hdf_file['analysis']['times'], dtype="float64")[firstpoint:] 

38 indexes_analysis = np.intersect1d(times, times_analysis, return_indices=True)[1] 

39 times = np.insert(times, indexes_analysis, times_analysis) 

40 

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

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

43 if measureName == 'times': 

44 continue 

45 else: 

46 measureType, units = giveMeasureTypeUnits(measureName) 

47 

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

49 # Remove firstpoints 

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

51 formattedForecastData = np.moveaxis(hdf_file['forecast'][measureName], 0, -1)[firstpoint:] 

52 

53 if measureName == 'MF': 

54 lmax = hdf_file.attrs['Lb'] 

55 elif measureName == 'U': 

56 lmax = hdf_file.attrs['Lu'] 

57 else: 

58 lmax = hdf_file.attrs['Lsv'] 

59 

60 # Insert the forecast at analysis times 

61 formattedData = np.insert(formattedData, indexes_analysis, formattedForecastData[indexes_analysis], axis=0) 

62 

63 if keepRealisations: 

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

65 formattedData, times=times) 

66 else: 

67 meanData = formattedData.mean(axis=2) 

68 rmsData = formattedData.std(axis=2) 

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

70 meanData, rmsData, times) 

71 

72 # Returns 0 if everything went well 

73 return 0