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

1import os 

2import numpy as np 

3 

4 

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

6 """ Loading function for Coupled-Earth data. Also adds the data to the dataModel. 

7 

8 :param dataDirectory: directory where the data is located 

9 :type dataDirectory: str (path) 

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

11 :type dataModel: Model 

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

13 :type keepRealisations: bool (default: False) 

14 :return: 0 if everything went well 

15 :rtype: int 

16 """ 

17 print("read MF CED data...") 

18 # Reading MF @ CMB in NS units 

19 gnmE = np.loadtxt(os.path.join(dataDirectory, "Cor_gnm_r_o.dat")) 

20 

21 print("... CED gnm data read !!") 

22 

23 # Detect lmax 

24 lmax = (-2+np.sqrt(4+4*gnmE.shape[1]))/2 

25 print("lmax B=", lmax) 

26 if int(lmax) != lmax: 

27 raise ValueError("Data length %i does not correspond to lmax*(lmax+2)" % gnmE.shape[1]) 

28 else: 

29 lmax = int(lmax) 

30 

31 nt = gnmE.shape[0] 

32 

33 # Load times 

34 times = np.linspace(1, nt, nt) 

35 

36 data = gnmE #* 10**9 # scale to nT ?? 

37 print("number of epochs:", nt) 

38 dataModel.addMeasure("MF", "MF", lmax, "nT", data, times=times) 

39 

40################# deals with gnm below CMB 

41 

42 print("read MF mid-path data below CMB...") 

43 # Reading MF @ CMB in NS units 

44 gnmE = np.loadtxt(os.path.join(dataDirectory,"Cor_gnm_r_o-delta.dat")) 

45 

46 print("... mid-path gnm data below CMB read !!") 

47 

48 # Detect lmax 

49 lmax = (-2+np.sqrt(4+4*gnmE.shape[1]))/2 

50 print("lmax B=",lmax) 

51 if int(lmax) != lmax: 

52 raise ValueError("Data length %i does not correspond to lmax*(lmax+2)" % gnmE.shape[1]) 

53 else: 

54 lmax = int(lmax) 

55 

56 n=lmax*(lmax+2) 

57 nt=gnmE.shape[0] 

58 

59 # Load times 

60 times=np.linspace(1, nt, nt) 

61 

62 data = gnmE #* 10**9 # scale to nT ?? 

63 print("number of epochs:", nt) 

64 dataModel.addMeasure("MFsub", "MF", lmax, "nT", data, times=times) 

65 

66################# deals with dgnm/dt 

67 

68 print("read SV mid-path data...") 

69 # Reading MF @ CMB in NS units 

70 dgnmE = np.loadtxt(os.path.join(dataDirectory,"Cor_dgnm_r_o.dat")) 

71 

72 print("... mid-path dgnm/dt data read !!") 

73 

74 # Detect lmax 

75 lmax = (-2+np.sqrt(4+4*dgnmE.shape[1]))/2 

76 print("lmax B=",lmax) 

77 if int(lmax) != lmax: 

78 raise ValueError("Data length %i does not correspond to lmax*(lmax+2)" % dgnmE.shape[1]) 

79 else: 

80 lmax = int(lmax) 

81 

82 n=lmax*(lmax+2) 

83 nt=dgnmE.shape[0] 

84 

85 # Load times 

86 times=np.linspace(1, nt, nt) 

87 

88 data = dgnmE #* 10**9 # scale to nT ?? 

89 print("number of epochs:", nt) 

90 dataModel.addMeasure("SV", "SV", lmax, "nT/yr", data, times=times) 

91 

92################# now deals with tnm, snm 

93 

94 # Reading flow U @ CMB in NS units 

95 tnmsnm = np.loadtxt(os.path.join(dataDirectory, "Cor_tnmsnm_r_o.dat")) 

96 

97 print("... CED tnm+snm data read !!") 

98 

99 # Detect lmax 

100 lmax = (-2+np.sqrt(4+4*(tnmsnm.shape[1])/2))/2 

101 print("lmax U=", lmax) 

102 if int(lmax) != lmax: 

103 raise ValueError("Data length %i does not correspond to 2*lmax*(lmax+2)" % data.shape[1]) 

104 else: 

105 lmax = int(lmax) 

106 

107 # /!\ not the same sign convention from JA for tnm,snm (need to invert) 

108 data = tnmsnm * (-1) # scale in km/yr 

109 dataModel.addMeasure("U", "U", lmax, "km/yr", data, times=times) 

110 

111 return 0