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

1import os 

2import numpy as np 

3from scipy.interpolate import BSpline 

4import glob 

5 

6 

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

8 """ Loading function for CHAOS-6 model. Also adds the data to the dataModel. 

9 

10 :param dataDirectory: directory where the data is located 

11 :type dataDirectory: str (path) 

12 :param dataModel: Model to which the data should be added 

13 :type dataModel: Model 

14 :param keepRealisations: indicating if realisations should be kept or averaged (not used) 

15 :type keepRealisations: bool (default: False) 

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

17 :rtype: int 

18 """ 

19 # Custom lmax for displaying core field & secular variation (SV) 

20 lmax_custom = {"MF": 13, "SV": 16} 

21 dataDirectory = os.path.normpath(dataDirectory) 

22 chaos_files = glob.glob(os.path.join(dataDirectory, "CHAOS-*_spline-coefficients.dat")) 

23 

24 if len(chaos_files) == 0: 

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

26 

27 # Take the most recent extension 

28 fname = sorted(chaos_files, reverse=True)[0] 

29 

30 # Open the file and print the header 

31 f = open(fname) 

32 print(f.readline()) 

33 

34 # Read parameters 

35 lmax, nspl, jorder = [int(x) for x in f.readline().split()] 

36 print(lmax, nspl, jorder) 

37 

38 # Read times 

39 knotstimes = np.zeros((nspl+jorder,)) 

40 for i in range(nspl+jorder): 

41 knotstimes[i] = float(f.readline()) 

42 

43 # Read g spline coefficient values 

44 gt = np.zeros((lmax*(lmax+2), nspl)) 

45 for i in range(nspl): 

46 for j in range(lmax*(lmax+2)): 

47 gt[j, i] = float(f.readline()) 

48 

49 f.close() 

50 # End of reading of the CHAOS file. 

51 

52 times = np.copy(knotstimes[jorder-1:-(jorder-1)]) 

53 ntimes = times.shape[0] 

54 g = np.zeros((ntimes, lmax*(lmax+2))) 

55 dg = np.zeros((ntimes, lmax*(lmax+2))) 

56 for j in range(lmax*(lmax+2)): 

57 # Build g and its derivative from splines 

58 spl = BSpline(knotstimes, gt[j, :], jorder-1) 

59 dspl = spl.derivative() 

60 g[:, j] = spl(times) 

61 dg[:, j] = dspl(times) 

62 

63 # Add data for MF and SV up to the custom lmax 

64 lmax_b = lmax_custom["MF"] 

65 lmax_sv = lmax_custom["SV"] 

66 dataModel.addMeasure("MF lmax={}".format(lmax_b), "MF", lmax_b, "nT", g, times=times) 

67 dataModel.addMeasure("SV lmax={}".format(lmax_sv), "SV", lmax_sv, "nT/yr", dg, times=times) 

68 # Add data for MF and SV up to the lmax from CHAOS 

69 dataModel.addMeasure("MF", "MF", lmax, "nT", g, times=times) 

70 dataModel.addMeasure("SV", "SV", lmax, "nT/yr", dg, times=times) 

71 

72 return 0