Surface NMR processing and inversion GUI
Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290
  1. from akvo.tressel.SlidesPlot import *
  2. import numpy as np
  3. import sys
  4. import matplotlib.pyplot as plt
  5. import cmocean
  6. from pylab import meshgrid
  7. from akvo.tressel.logbarrier import *
  8. import yaml,os
  9. from matplotlib.colors import LogNorm
  10. from matplotlib.colors import LightSource
  11. from matplotlib.ticker import ScalarFormatter
  12. from matplotlib.ticker import MaxNLocator
  13. from matplotlib.ticker import AutoMinorLocator
  14. from matplotlib.ticker import LogLocator
  15. from matplotlib.ticker import FormatStrFormatter
  16. from matplotlib.colors import Normalize
  17. import cmocean
  18. from akvo.tressel.lemma_yaml import *
  19. import pandas as pd
  20. def buildKQT(K0,tg,T2Bins):
  21. """
  22. Constructs a QT inversion kernel from an initial amplitude one.
  23. """
  24. nlay, nq = np.shape(K0)
  25. nt2 = len(T2Bins)
  26. nt = len(tg)
  27. KQT = np.zeros( ( nq*nt,nt2*nlay) )
  28. for iq in range(nq):
  29. for it in range(nt):
  30. for ilay in range(nlay):
  31. for it2 in range(nt2):
  32. #KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-((10+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
  33. KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-((10+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
  34. return KQT
  35. def loadAkvoData(fnamein, chan):
  36. """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be
  37. corrected in future kernel calculations. The current was reported but not the pulse length.
  38. """
  39. fname = (os.path.splitext(fnamein)[0])
  40. with open(fnamein, 'r') as stream:
  41. try:
  42. AKVO = (yaml.load(stream, Loader=yaml.Loader))
  43. except yaml.YAMLError as exc:
  44. print(exc)
  45. exit()
  46. Z = np.zeros( (AKVO.nPulseMoments, AKVO.Gated["Pulse 1"]["abscissa"].size ) )
  47. ZS = np.zeros( (AKVO.nPulseMoments, AKVO.Gated["Pulse 1"]["abscissa"].size ) )
  48. for q in range(AKVO.nPulseMoments):
  49. Z[q] = AKVO.Gated["Pulse 1"][chan]["Q-"+str(q) + " CA"].data
  50. if chan == "Chan. 1":
  51. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  52. elif chan == "Chan. 2":
  53. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  54. elif chan == "Chan. 3":
  55. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  56. elif chan == "Chan. 4":
  57. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  58. else:
  59. print("DOOM!!!")
  60. exit()
  61. #Z *= 1e-9
  62. #ZS *= 1e-9
  63. J = AKVO.Pulses["Pulse 1"]["current"].data
  64. J = np.append(J,J[-1]+(J[-1]-J[-2]))
  65. Q = AKVO.pulseLength[0]*J
  66. return Z, ZS, AKVO.Gated["Pulse 1"]["abscissa"].data #, Q
  67. def catLayers(K0):
  68. K = np.zeros( (len(K0.keys()), len(K0["layer-0"].data)) , dtype=complex )
  69. for lay in range(len(K0.keys())):
  70. #print(K0["layer-"+str(lay)].data) # print (lay)
  71. K[lay] =K0["layer-"+str(lay)].data # print (lay)
  72. return 1e9*K # invert in nV
  73. def loadK0(fname):
  74. """ Loads in initial amplitude kernel
  75. """
  76. print("loading K0", fname)
  77. with open(fname) as f:
  78. K0 = yaml.load(f, Loader=yaml.Loader)
  79. K = catLayers(K0.K0)
  80. ifaces = np.array(K0.Interfaces.data)
  81. return ifaces, np.abs(K)
  82. def main():
  83. if (len (sys.argv) < 2):
  84. print ("akvoQT invertParameters.yaml")
  85. exit()
  86. with open(sys.argv[1], 'r') as stream:
  87. try:
  88. cont = (yaml.load(stream, Loader=yaml.Loader))
  89. except yaml.YAMLError as exc:
  90. print(exc)
  91. ###############################################
  92. # Load in data
  93. ###############################################
  94. V = []
  95. VS = []
  96. tg = 0
  97. for dat in cont['data']:
  98. for ch in cont['data'][dat]['channels']:
  99. print("dat", dat, "ch", ch)
  100. v,vs,tg = loadAkvoData(dat, ch)
  101. V.append(v)
  102. VS.append(vs)
  103. for iv in range(1, len(V)):
  104. V[0] = np.concatenate( (V[0], V[iv]) )
  105. VS[0] = np.concatenate( (VS[0], VS[iv]) )
  106. V = V[0]
  107. VS = VS[0]
  108. ###############################################
  109. # Load in kernels
  110. ###############################################
  111. K0 = []
  112. for kern in cont["K0"]:
  113. ifaces,k0 = loadK0( kern )
  114. K0.append(k0)
  115. for ik in range(1, len(K0)):
  116. K0[0] = np.concatenate( (K0[0].T, K0[ik].T) ).T
  117. K0 = K0[0]
  118. #plt.matshow(K0)
  119. ###################
  120. # VERY Simple DOI #
  121. maxq = np.argmax(K0, axis=1)
  122. maxK = .1 * K0[ np.arange(0,len(ifaces)-1), maxq ] # 10% water is arbitrary
  123. SNR = maxK / (VS[0][0])
  124. #SNR[SNR>1] = 1
  125. SNRidx = len(ifaces)-3
  126. while SNR[SNRidx] > SNR[SNRidx+1] and SNRidx > 2:
  127. SNRidx -= 1
  128. #print("IDX", SNRidx)
  129. #plt.plot(ifaces[0:-1], SNR)
  130. #plt.plot(ifaces[0:-1][SNRidx], SNR[SNRidx], '.',markersize=12)
  131. #plt.gca().axhline(y=VS[0][0], xmin=0, xmax=ifaces[-1], color='r')
  132. #plt.gca().axhline(y=1, xmin=0, xmax=ifaces[-1], color='r')
  133. #K0T = np.dot(K0, K0.T)
  134. #K0T = np.dot(K0, np.dot( VS[0][0]* np.eye(len(ifaces)-1,1), K0.T) )
  135. #K0T = np.dot(K0, np.dot( 1/(VS[0][0]**2) * np.eye(np.shape(K0)[1]), K0.T) )
  136. #plt.matshow(0.05 * np.sqrt(K0T))
  137. #plt.colorbar()
  138. #plt.plot(ifaces[0:-1], np.diag( 0.01* np.sqrt(K0T)))
  139. #print(np.shape(K0T), len(ifaces)-1)
  140. #plt.show()
  141. #exit()
  142. ###############################################
  143. # Build full kernel
  144. ###############################################
  145. T2Bins = np.logspace( np.log10(cont["T2Bins"]["low"]), np.log10(cont["T2Bins"]["high"]), cont["T2Bins"]["number"], endpoint=True, base=10)
  146. KQT = buildKQT(K0,tg,T2Bins)
  147. ###############################################
  148. # Invert
  149. ###############################################
  150. print("Calling inversion", flush=True)
  151. inv, ibreak, errn, phim, phid, mkappa = logBarrier(KQT, np.ravel(V), T2Bins, "lcurve", MAXITER=150, sigma=np.ravel(VS), alpha=1e6, smooth="Smallest" )
  152. ###############################################
  153. # Appraise
  154. ###############################################
  155. pre = np.dot(KQT,inv)
  156. PRE = np.reshape( pre, np.shape(V) )
  157. plt.matshow(PRE, cmap='Blues')
  158. plt.gca().set_title("predicted")
  159. plt.colorbar()
  160. DIFF = (PRE-V) / VS
  161. md = np.max(np.abs(DIFF))
  162. plt.matshow(DIFF, cmap=cmocean.cm.balance, vmin=-md, vmax=md)
  163. plt.gca().set_title("misfit / $\widehat{\sigma}$")
  164. plt.colorbar()
  165. plt.matshow(V, cmap='Blues')
  166. plt.gca().set_title("observed")
  167. plt.colorbar()
  168. T2Bins = np.append( T2Bins, T2Bins[-1] + (T2Bins[-1]-T2Bins[-2]) )
  169. INV = np.reshape(inv, (len(ifaces)-1,cont["T2Bins"]["number"]) )
  170. #alphas = np.tile(SNR, (len(T2Bins)-1,1))
  171. #colors = Normalize(1e-6, np.max(INV.T), clip=True)(INV.T)
  172. #colors = cmocean.cm.tempo(colors)
  173. ##colors[..., -1] = alphas
  174. #print(np.shape(colors))
  175. #print(np.shape(INV.T))
  176. #greys = np.full((*(INV.T).shape, 3), 70, dtype=np.uint8)
  177. Y,X = meshgrid( ifaces, T2Bins )
  178. fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  179. ax1 = fig.add_axes( [.2,.15,.6,.7] )
  180. im = ax1.pcolor(X, Y, INV.T, cmap=cmocean.cm.tempo) #cmap='viridis')
  181. #im = ax1.pcolor(X[0:SNRidx,:], Y[0:SNRidx,:], INV.T[0:SNRidx,:], cmap=cmocean.cm.tempo) #cmap='viridis')
  182. #im = ax1.pcolor(X[SNRidx::,:], Y[SNRidx::,:], INV.T[SNRidx::,:], cmap=cmocean.cm.tempo, alpha=.5) #cmap='viridis')
  183. #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
  184. #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
  185. #ax1.axhline( y=ifaces[SNRidx], xmin=T2Bins[0], xmax=T2Bins[-1], color='black' )
  186. im.set_edgecolor('face')
  187. ax1.set_xlim( T2Bins[0], T2Bins[-1] )
  188. ax1.set_ylim( ifaces[-1], ifaces[0] )
  189. cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
  190. cb.locator = MaxNLocator( nbins = 4)
  191. cb.ax.yaxis.set_offset_position('left')
  192. cb.update_ticks()
  193. ax1.set_xlabel(u"$T_2^*$ (ms)")
  194. ax1.set_ylabel(u"depth (m)")
  195. ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
  196. ax1.get_yaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
  197. ax1.xaxis.set_major_locator( MaxNLocator(nbins = 4) )
  198. #ax1.xaxis.set_label_position('top')
  199. ax2 = ax1.twiny()
  200. ax2.plot( np.sum(INV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 , color='red' )
  201. ax2.set_xlabel(u"total water (m$^3$/m$^3$)")
  202. ax2.set_ylim( ifaces[-1], ifaces[0] )
  203. ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )
  204. ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
  205. ax2.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed' )
  206. #ax2.xaxis.set_label_position('bottom')
  207. plt.savefig("akvoInversion.pdf")
  208. #############
  209. # water plot#
  210. fig2 = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  211. ax = fig2.add_axes( [.2,.15,.6,.7] )
  212. # Bound water cutoff
  213. Bidx = T2Bins[0:-1]<33.0
  214. twater = np.sum(INV, axis=1)
  215. bwater = np.sum(INV[:,Bidx], axis=1)
  216. ax.plot( twater, (ifaces[0:-1]+ifaces[1::])/2, label="NMR total water", color='blue' )
  217. ax.plot( bwater, (ifaces[0:-1]+ifaces[1::])/2, label="NMR bound water", color='green' )
  218. ax.fill_betweenx((ifaces[0:-1]+ifaces[1::])/2 , twater, bwater, where=twater >= bwater, facecolor='blue', alpha=.5)
  219. ax.fill_betweenx((ifaces[0:-1]+ifaces[1::])/2 , bwater, 0, where=bwater >= 0, facecolor='green', alpha=.5)
  220. ax.set_xlabel(r"$\theta_N$ (m$^3$/m$^3$)")
  221. ax.set_ylabel(r"depth (m)")
  222. ax.set_ylim( ifaces[-1], ifaces[0] )
  223. ax.set_xlim( 0, ax.get_xlim()[1] )
  224. ax.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed' )
  225. plt.savefig("akvoInversionWC.pdf")
  226. plt.legend()
  227. # Report results into a text file
  228. fr = pd.DataFrame( INV, columns=T2Bins[0:-1] )
  229. fr.insert(0, "layer top", ifaces[0:-1] )
  230. fr.insert(1, "layer bottom", ifaces[1::] )
  231. fr.insert(2, "NMR total water", np.sum(INV, axis=1) )
  232. fr.insert(3, "NMR bound water", bwater )
  233. fr.insert(4, "Layer SNR", SNR )
  234. fr.to_csv("akvoInversion.csv")
  235. #fr.to_excel("akvoInversion.xlsx")
  236. plt.show()
  237. if __name__ == "__main__":
  238. main()