Surface NMR processing and inversion GUI
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182
  1. import os, sys
  2. import numpy as np
  3. from ruamel import yaml
  4. import pyLemma.LemmaCore as lc
  5. import pyLemma.Merlin as mrln
  6. import pyLemma.FDEM1D as em1d
  7. import numpy as np
  8. import matplotlib.pyplot as plt
  9. import seaborn as sns
  10. sns.set(style="ticks")
  11. from ruamel import yaml
  12. #import cmocean
  13. #from SEGPlot import *
  14. #from matplotlib.ticker import FormatStrFormatter
  15. #import matplotlib.ticker as plticker
  16. # Converts Lemma/Merlin/Akvo serialized Eigen arrays into numpy ones for use by Python
  17. class VectorXr(yaml.YAMLObject):
  18. """
  19. Converts Lemma/Merlin/Akvo serialized Eigen arrays into numpy ones for use by Python
  20. """
  21. yaml_tag = u'VectorXr'
  22. def __init__(self, array):
  23. self.size = np.shape(array)[0]
  24. self.data = array.tolist()
  25. def __repr__(self):
  26. # Converts to numpy array on import
  27. return "np.array(%r)" % (self.data)
  28. class AkvoData(yaml.YAMLObject):
  29. """
  30. Reads an Akvo serialized dataset into a standard python dictionary
  31. """
  32. yaml_tag = u'AkvoData'
  33. def __init__(self, array):
  34. pass
  35. #self.size = np.shape(array)[0]
  36. #self.Imp = array.tolist()
  37. def __repr__(self):
  38. # Converts to a dictionary with Eigen vectors represented as Numpy arrays
  39. return self
  40. def loadAkvoData(fnamein):
  41. """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be
  42. corrected in future kernel calculations. The current was reported but not the pulse length.
  43. """
  44. fname = (os.path.splitext(fnamein)[0])
  45. with open(fnamein, 'r') as stream:
  46. try:
  47. AKVO = (yaml.load(stream, Loader=yaml.Loader))
  48. except yaml.YAMLError as exc:
  49. print(exc)
  50. return AKVO
  51. def main():
  52. if len(sys.argv) < 2:
  53. print ("usage python calcAkvoKernel.py AkvoDataset.yaml Coil1.yaml kparams.yaml SaveString.yaml " )
  54. print ("usage akvoKO AkvoDataset.yaml kparams.yaml SaveString.yaml " )
  55. exit()
  56. print("Loading data")
  57. AKVO = loadAkvoData(sys.argv[1])
  58. print("Building Kernel")
  59. B_inc = AKVO.META["B_0"]["inc"]
  60. B_dec = AKVO.META["B_0"]["dec"]
  61. B0 = AKVO.META["B_0"]["intensity"]
  62. fT = AKVO.transFreq
  63. #gamma = 2.67518e8
  64. #B0 = (fL*2.*np.pi) /gamma * 1e9
  65. # read in kernel params
  66. kparams = loadAkvoData( sys.argv[2] )
  67. Kern = mrln.KernelV0()
  68. TX = []
  69. for tx in kparams['txCoils']:
  70. Coil1 = em1d.PolygonalWireAntenna.DeSerialize( tx )
  71. Coil1.SetNumberOfFrequencies(1)
  72. Coil1.SetFrequency(0, fT)
  73. Coil1.SetCurrent(1.)
  74. Kern.PushCoil( tx.split('.yml')[0], Coil1 )
  75. TX.append( tx.split('.yml')[0] )
  76. RX = []
  77. for rx in kparams['rxCoils']:
  78. if rx not in kparams['txCoils']:
  79. print("new recv")
  80. Coil1 = em1d.PolygonalWireAntenna.DeSerialize( rx )
  81. Coil1.SetNumberOfFrequencies(1)
  82. Coil1.SetFrequency(0, fT)
  83. Coil1.SetCurrent(1.)
  84. Kern.PushCoil( rx.split('.yml')[0], Coil1 )
  85. else:
  86. print("reuse tx coil")
  87. RX.append( rx.split('.yml')[0] )
  88. ## TODO
  89. # pass this in...
  90. print("Building layered earth model")
  91. lmod = em1d.LayeredEarthEM()
  92. nlay = len(kparams["sigs"])
  93. sigs = np.array(kparams["sigs"])
  94. tops = np.array(kparams["tops"])
  95. bots = np.array(kparams["bots"])
  96. if ( (len(tops)-1) != len(bots)):
  97. print("Layer mismatch")
  98. exit()
  99. thicks = bots - tops[0:-1]
  100. lmod.SetNumberOfLayers(nlay + 1)
  101. lmod.SetLayerThickness(thicks)
  102. lmod.SetLayerConductivity( np.concatenate( ( [0.0], sigs ) ))
  103. #lmod.SetNumberOfLayers(4)
  104. #lmod.SetLayerThickness([15.49, 28.18])
  105. #lmod.SetLayerConductivity([0.0, 1./16.91, 1./24.06, 1./33.23])
  106. lmod.SetMagneticFieldIncDecMag( B_inc, B_dec, B0, lc.NANOTESLA )
  107. Kern.SetLayeredEarthEM( lmod );
  108. Kern.SetIntegrationSize( (kparams["size_n"], kparams["size_e"], kparams["size_d"]) )
  109. Kern.SetIntegrationOrigin( (kparams["origin_n"], kparams["origin_e"], kparams["origin_d"]) )
  110. Kern.SetTolerance( 1e-9*kparams["branchTol"] )
  111. Kern.SetMinLevel( kparams["minLevel"] )
  112. Kern.SetMaxLevel( kparams["maxLevel"] )
  113. Kern.SetHankelTransformType( lc.FHTKEY201 )
  114. Kern.AlignWithAkvoDataset( sys.argv[1] )
  115. if str(kparams["Lspacing"]).strip() == "Geometric":
  116. thick = np.geomspace(kparams["thick1"], kparams["thickN"], num=kparams["nLay"])
  117. elif str(kparams["Lspacing"]) == "Log":
  118. thick = np.logspace(kparams["thick1"], kparams["thickN"], num=kparams["nLay"])
  119. elif str(kparams["Lspacing"]) == "Linear":
  120. thick = np.linspace(kparams["thick1"], kparams["thickN"], num=kparams["nLay"])
  121. else:
  122. print("DOOOM!, in calcAkvoKernel layer spacing was not <Geometric>, <Log>, or <Linear>")
  123. print( str(kparams["Lspacing"]) )
  124. exit()
  125. print( np.array(kparams["origin_d"]) )
  126. print( np.cumsum(thick)[0:-1] )
  127. iface = np.concatenate( (np.array( [kparams["origin_d"]] ), kparams["origin_d"]+np.cumsum(thick)[0:-1]) )
  128. Kern.SetDepthLayerInterfaces(iface)
  129. #Kern.SetDepthLayerInterfaces(np.geomspace(1, 110, num=40))
  130. #Kern.SetDepthLayerInterfaces(np.linspace(1, 110, num=50))
  131. #Kern.SetDepthLayerInterfaces(np.geomspace(1, 110, num=40))
  132. # autAkvoDataNode = YAML::LoadFile(argv[4]);
  133. # Kern->AlignWithAkvoDataset( AkvoDataNode );
  134. #Kern.CalculateK0( ["Coil 1"], ["Coil 1"], False )
  135. Kern.CalculateK0( TX, RX, False )
  136. #yml = open( 'test' + str(Kern.GetTolerance()) + '.yaml', 'w')
  137. yml = open( sys.argv[3], 'w' )
  138. print(Kern, file=yml)
  139. #
  140. K0 = Kern.GetKernel()
  141. plt.matshow(np.abs(K0))
  142. plt.show()
  143. if __name__ == "__main__":
  144. main()