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

1import os 

2import h5py 

3import numpy as np 

4import glob 

5 

6def load(dataDirectory, dataModel, keepRealisations=False): 

7 

8 measures_to_load = ['MF', 'SV'] 

9 

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

11 if len(hdf5_files) == 0: 

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

13 

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

15 hdf_filename = hdf5_files[0] 

16 print('Reading:', hdf_filename) 

17 

18 lmax = 14 

19 

20 with h5py.File(hdf_filename) as hdf_file: 

21 times = np.array(hdf_file['times']) 

22 

23 for measureName in measures_to_load: 

24 

25 if measureName == "SV": 

26 measureData = hdf_file['dgnm'] 

27 units = "nT/yr" 

28 elif measureName == "MF": 

29 measureData = hdf_file['gnm'] 

30 units = "nT" 

31 

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

33 # Remove firstpoints 

34 formattedData = np.moveaxis(measureData, 0, -1) 

35 assert lmax * (lmax + 2) == formattedData.shape[1], 'lmax or shape of data is wrong' 

36 

37 print('keepRealisations', keepRealisations) 

38 if keepRealisations: 

39 # If realisation must be kept, add the data to the model as read 

40 dataModel.addMeasure(measureName, measureName, 

41 lmax, units, formattedData, times=times) 

42 else: 

43 # Else, average the data read on realisations 

44 meanData = formattedData.mean(axis=2) 

45 rmsData = formattedData.std(axis=2) 

46 dataModel.addMeasure(measureName, measureName, lmax, units, 

47 meanData, rmsData, times) 

48 

49 return 0