Surface NMR processing and inversion GUI
Du kan inte välja fler än 25 ämnen Ämnen måste starta med en bokstav eller siffra, kan innehålla bindestreck ('-') och vara max 35 tecken långa.

invertTA.py 26KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683
  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. import multiprocessing
  10. import itertools
  11. from scipy.linalg import svd
  12. from matplotlib.backends.backend_pdf import PdfPages
  13. from matplotlib.colors import LogNorm
  14. from matplotlib.colors import LightSource
  15. from matplotlib.ticker import ScalarFormatter
  16. from matplotlib.ticker import MaxNLocator
  17. from matplotlib.ticker import AutoMinorLocator
  18. from matplotlib.ticker import LogLocator
  19. from matplotlib.ticker import FormatStrFormatter
  20. from matplotlib.colors import Normalize
  21. import cmocean
  22. from akvo.tressel.lemma_yaml import *
  23. from akvo.tressel import nonlinearinv as nl
  24. import pandas as pd
  25. # For Mihai, but consider letting fonts be user selected.
  26. #import matplotlib.pyplot as plt
  27. #plt.rcParams["font.family"] = "arial"
  28. import matplotlib.colors as colors
  29. # From https://stackoverflow.com/questions/18926031/how-to-extract-a-subset-of-a-colormap-as-a-new-colormap-in-matplotlib
  30. def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
  31. new_cmap = colors.LinearSegmentedColormap.from_list(
  32. 'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
  33. cmap(np.linspace(minval, maxval, n)))
  34. return new_cmap
  35. def buildKQT(K0,tg,T2Bins):
  36. """
  37. Constructs a QT inversion kernel from an initial amplitude one.
  38. """
  39. nlay, nq = np.shape(K0)
  40. nt2 = len(T2Bins)
  41. nt = len(tg)
  42. KQT = np.zeros( ( nq*nt,nt2*nlay), dtype=np.complex128 )
  43. for iq in range(nq):
  44. for it in range(nt):
  45. for ilay in range(nlay):
  46. for it2 in range(nt2):
  47. #KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-((10+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
  48. KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-((10+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
  49. return KQT
  50. def loadAkvoData(fnamein, chan):
  51. """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be
  52. corrected in future kernel calculations. The current was reported but not the pulse length.
  53. """
  54. fname = (os.path.splitext(fnamein)[0])
  55. with open(fnamein, 'r') as stream:
  56. try:
  57. AKVO = (yaml.load(stream, Loader=yaml.Loader))
  58. except yaml.YAMLError as exc:
  59. print(exc)
  60. exit()
  61. Z = np.zeros( (AKVO.nPulseMoments, AKVO.Gated["Pulse 1"]["abscissa"].size ) )
  62. ZS = np.zeros( (AKVO.nPulseMoments, AKVO.Gated["Pulse 1"]["abscissa"].size ) )
  63. for q in range(AKVO.nPulseMoments):
  64. Z[q] = AKVO.Gated["Pulse 1"][chan]["Q-"+str(q) + " CA"].data
  65. if chan == "Chan. 1":
  66. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  67. elif chan == "Chan. 2":
  68. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  69. elif chan == "Chan. 3":
  70. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  71. elif chan == "Chan. 4":
  72. ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
  73. else:
  74. print("DOOM!!!")
  75. exit()
  76. #Z *= 1e-9
  77. #ZS *= 1e-9
  78. J = np.array(AKVO.Pulses["Pulse 1"]["current"].data)
  79. #J = np.append(J,J[-1]+(J[-1]-J[-2])) # This was a pcolor hack, get rid of this
  80. #print("J", J)
  81. Q = AKVO.pulseLength[0]*J
  82. return Z, ZS, AKVO.Gated["Pulse 1"]["abscissa"].data, Q
  83. def catLayers(K0):
  84. K = np.zeros( (len(K0.keys()), len(K0["layer-0"].data)) , dtype=complex )
  85. for lay in range(len(K0.keys())):
  86. #print(K0["layer-"+str(lay)].data) # print (lay)
  87. K[lay] =K0["layer-"+str(lay)].data # print (lay)
  88. return 1e9*K # invert in nV
  89. def loadK0(fname):
  90. """ Loads in initial amplitude kernel
  91. """
  92. print("loading K0", fname)
  93. with open(fname) as f:
  94. K0 = yaml.load(f, Loader=yaml.Loader)
  95. K = catLayers(K0.K0)
  96. ifaces = np.array(K0.Interfaces.data)
  97. return ifaces, K
  98. #return ifaces, np.abs(K)
  99. def invertDelta(G, V_n, T2Bins, sig, alphastar):
  100. """ helper function that simply calls logBarrier, simplfies parallel execution
  101. """
  102. model = logBarrier(G, V_n, T2Bins, "Single", MAXITER=1, sigma=sig, alpha=alphastar, smooth="Smallest")
  103. return model
  104. def main():
  105. if (len (sys.argv) < 2):
  106. print ("akvoQT invertParameters.yaml")
  107. exit()
  108. with open(sys.argv[1], 'r') as stream:
  109. try:
  110. cont = (yaml.load(stream, Loader=yaml.Loader))
  111. except yaml.YAMLError as exc:
  112. print(exc)
  113. exit(1)
  114. ###############################################
  115. # Load in data
  116. ###############################################
  117. V = []
  118. VS = []
  119. QQ = []
  120. tg = 0
  121. for dat in cont['data']:
  122. for ch in cont['data'][dat]['channels']:
  123. print("dat", dat, "ch", ch)
  124. v,vs,tg,Q = loadAkvoData(dat, ch)
  125. V.append(v)
  126. VS.append(vs)
  127. QQ.append(Q)
  128. for iv in range(1, len(V)):
  129. V[0] = np.concatenate( (V[0], V[iv]) )
  130. VS[0] = np.concatenate( (VS[0], VS[iv]) )
  131. V = V[0]
  132. VS = VS[0]
  133. ###############################################
  134. # Load in kernels
  135. ###############################################
  136. K0 = []
  137. for kern in cont["K0"]:
  138. ifaces,k0 = loadK0( kern )
  139. K0.append(k0)
  140. for ik in range(1, len(K0)):
  141. K0[0] = np.concatenate( (K0[0].T, K0[ik].T) ).T
  142. K0 = K0[0]
  143. #np.save("ifaces", ifaces)
  144. #exit()
  145. #plt.matshow(np.real(K0))
  146. #plt.show()
  147. #exit()
  148. ##############################################################
  149. # VERY Simple Sensitivity based calc. of noise per layer #
  150. # minimally useful, but retained for backwards compatibility #
  151. maxq = np.argmax(np.abs(K0), axis=1)
  152. maxK = .1 * np.abs(K0)[ np.arange(0,len(ifaces)-1), maxq ] # 10% water is arbitrary
  153. SNR = maxK / (VS[0][0])
  154. ###############################################
  155. # Build full kernel
  156. ###############################################
  157. T2Bins = np.logspace( np.log10(cont["T2Bins"]["low"]), np.log10(cont["T2Bins"]["high"]), cont["T2Bins"]["number"], endpoint=True, base=10)
  158. T2Bins2 = np.append( T2Bins, T2Bins[-1] + (T2Bins[-1]-T2Bins[-2]) )
  159. NT2 = len(T2Bins)
  160. KQT = np.real(buildKQT(np.abs(K0),tg,T2Bins))
  161. ###############################################
  162. # Linear Inversion
  163. ###############################################
  164. print("Calling inversion", flush=True)
  165. inv, ibreak, errn, phim, phid, mkappa, Wd, Wm, alphastar = logBarrier(KQT, np.ravel(V), T2Bins, "lcurve", MAXITER=150, sigma=np.ravel(VS), alpha=1e7, smooth="Smallest" )
  166. ################################
  167. # Summary plots, Data Space #
  168. ################################
  169. # TODO, need to clean this up for the case of multiple channels! Each channel should be a new row. It will be ugly, but important
  170. # TODO, loop over channels
  171. ich = 0
  172. for ch in cont['data'][dat]['channels']:
  173. figx = plt.figure( figsize=(pc2in(42.0),pc2in(22.)) )
  174. ax1 = figx.add_axes([.100, .15, .200, .70])
  175. ax2 = figx.add_axes([.325, .15, .200, .70]) # shifted to make room for shared colourbar
  176. axc1= figx.add_axes([.550, .15, .025, .70]) # shifted to make room for shared colourbar
  177. ax3 = figx.add_axes([.670, .15, .200, .70])
  178. axc2= figx.add_axes([.895, .15, .025, .70]) # shifted to make room for shared colourbar
  179. ax3.set_yscale('log')
  180. ax2.set_yscale('log')
  181. ax1.set_yscale('log')
  182. ax2.yaxis.set_ticklabels([])
  183. ax3.yaxis.set_ticklabels([])
  184. ax3.set_xscale('log')
  185. ax2.set_xscale('log')
  186. ax1.set_xscale('log')
  187. ax1.set_ylabel("Q (A $\cdot$ s)")
  188. ax1.set_xlabel("time (s)")
  189. ax2.set_xlabel("time (s)")
  190. ax3.set_xlabel("time (s)")
  191. #TT, QQQ = np.meshgrid(tg, np.ravel(QQ))
  192. TT, QQQ = np.meshgrid(tg, np.ravel(QQ[ich]))
  193. nq = np.shape(QQ[ich])[0] #- 1 # -1 to account for padding in pcolor, but this was changed
  194. nt = np.shape(tg)[0]
  195. ntq = nt*nq
  196. VV = V[ich*nq:ich*nq+nq,:] # slice this channel
  197. VVS = VS[ich*nq:ich*nq+nq,:] # slice this channel
  198. mmax = np.max(np.abs(VV))
  199. mmin = np.min(VV)
  200. obs = ax1.pcolor(TT, QQQ, VV, cmap=cmocean.cm.curl_r, vmin=-mmax, vmax=mmax, shading='auto') # pcolor edge not defined
  201. ax1.set_title("observed")
  202. pre = np.dot(KQT[ich*ntq:(ich+1)*ntq,:], inv)
  203. PRE = np.reshape( pre, np.shape(VV) )
  204. prem = ax2.pcolor(TT, QQQ, PRE, cmap=cmocean.cm.curl_r, vmin=-mmax, vmax=mmax,shading='auto' )
  205. ax2.set_title("predicted")
  206. cbar = plt.colorbar(prem, axc1)
  207. axc1.set_ylim( [np.min(VV), np.max(VV)] )
  208. cbar.outline.set_edgecolor(None)
  209. cbar.set_label('$V_N$ (nV)')
  210. DIFF = (PRE-VV) / VVS
  211. md = np.max(np.abs(DIFF))
  212. dim = ax3.pcolor(TT, QQQ, DIFF, cmap=cmocean.cm.balance, vmin=-md, vmax=md, shading='auto')
  213. ax3.set_title("misfit / $\widehat{\sigma}$")
  214. cbar2 = plt.colorbar(dim, axc2)
  215. #axc1.set_ylim( [np.min(V), np.max(V)] )
  216. cbar2.outline.set_edgecolor(None)
  217. cbar2.set_label('$V_N$ (nV)')
  218. #plt.colorbar(dim, ax3)
  219. figx.suptitle(ch + " linear Inversion")
  220. plt.savefig(ch + "dataspace.pdf")
  221. ich += 1
  222. ###############################################
  223. # Non-linear refinement!
  224. ###############################################
  225. nonLinearRefinement = cont['NonLinearRefinement']
  226. if nonLinearRefinement:
  227. KQTc = buildKQT(K0, tg, T2Bins)
  228. prec = np.abs(np.dot(KQTc, inv))
  229. phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
  230. print("PHID forward linear=", errn, "PHID forward nonlinear=", phidc/len(np.ravel(V)))
  231. res = nl.nonlinearinversion(inv, Wd, KQTc, np.ravel(V), Wm, alphastar )
  232. if res.success == True:
  233. INVc = np.reshape(res.x, (len(ifaces)-1,cont["T2Bins"]["number"]) )
  234. prec = np.abs(np.dot(KQTc, res.x))
  235. phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
  236. PREc = np.reshape( prec, np.shape(V) )
  237. print("PHID linear=", errn, "PHID nonlinear=", phidc/len(np.ravel(V)))
  238. while phidc/len(np.ravel(V)) > errn:
  239. phidc_old = phidc/len(np.ravel(V))
  240. #alphastar *= .9
  241. res = nl.nonlinearinversion(res.x, Wd, KQTc, np.ravel(V), Wm, alphastar )
  242. if res.success == True:
  243. INVc = np.reshape(res.x, (len(ifaces)-1,cont["T2Bins"]["number"]) )
  244. prec = np.abs(np.dot(KQTc, res.x))
  245. phidc = np.linalg.norm(np.dot(Wd,prec-np.ravel(V)))**2
  246. PREc = np.reshape( prec, np.shape(V) )
  247. print("PHID linear=", errn, "PHID nonlinear=", phidc/len(np.ravel(V)))
  248. else:
  249. nonLinearRefinement = False
  250. print("Non-linear inversion failed, results will not be shown")
  251. break
  252. if phidc_old - phidc/len(np.ravel(V)) < 0.005:
  253. print("Not making progress reducing misfit in nonlinear refinement")
  254. break
  255. # Turn this into a nice figure w/ shared axes etc.
  256. # plt.matshow(PREc, cmap='Blues')
  257. # plt.gca().set_title("nonlinear predicted")
  258. # plt.colorbar()
  259. #
  260. # DIFFc = (PREc-V) / VS
  261. # md = np.max(np.abs(DIFF))
  262. # plt.matshow(DIFFc, cmap=cmocean.cm.balance, vmin=-md, vmax=md)
  263. # plt.gca().set_title("nonlinear misfit / $\widehat{\sigma}$")
  264. # plt.colorbar()
  265. ################################
  266. # Summary plots, Data Space #
  267. ################################
  268. ich = 0
  269. for ch in cont['data'][dat]['channels']:
  270. figx = plt.figure( figsize=(pc2in(42.0),pc2in(22.)) )
  271. ax1 = figx.add_axes([.100, .15, .200, .70])
  272. ax2 = figx.add_axes([.325, .15, .200, .70]) # shifted to make room for shared colourbar
  273. axc1= figx.add_axes([.550, .15, .025, .70]) # shifted to make room for shared colourbar
  274. ax3 = figx.add_axes([.670, .15, .200, .70])
  275. axc2= figx.add_axes([.895, .15, .025, .70]) # shifted to make room for shared colourbar
  276. ax3.set_yscale('log')
  277. ax2.set_yscale('log')
  278. ax1.set_yscale('log')
  279. ax2.yaxis.set_ticklabels([])
  280. ax3.yaxis.set_ticklabels([])
  281. ax3.set_xscale('log')
  282. ax2.set_xscale('log')
  283. ax1.set_xscale('log')
  284. ax1.set_ylabel("Q (A $\cdot$ s)")
  285. ax1.set_xlabel("time (s)")
  286. ax2.set_xlabel("time (s)")
  287. ax3.set_xlabel("time (s)")
  288. #TT, QQQ = np.meshgrid(tg, np.ravel(QQ))
  289. TT, QQQ = np.meshgrid(tg, np.ravel(QQ[ich]))
  290. nq = np.shape(QQ[ich])[0] # to account for padding in pcolor, this changed
  291. nt = np.shape(tg)[0]
  292. ntq = nt*nq
  293. VV = V[ich*nq:ich*nq+nq,:] # slice this channel
  294. VVS = VS[ich*nq:ich*nq+nq,:] # slice this channel
  295. mmax = np.max(np.abs(VV))
  296. mmin = np.min(VV)
  297. obs = ax1.pcolor(TT, QQQ, VV, cmap=cmocean.cm.curl_r, vmin=-mmax, vmax=mmax, shading='auto')
  298. ax1.set_title("observed")
  299. ## Here neds to change
  300. pre = np.abs(np.dot(KQTc[ich*ntq:(ich+1)*ntq,:], inv))
  301. PRE = np.reshape( pre, np.shape(VV) )
  302. prem = ax2.pcolor(TT, QQQ, PRE, cmap=cmocean.cm.curl_r, vmin=-mmax, vmax=mmax, shading='auto' )
  303. ax2.set_title("predicted")
  304. cbar = plt.colorbar(prem, axc1)
  305. axc1.set_ylim( [np.min(VV), np.max(VV)] )
  306. cbar.outline.set_edgecolor(None)
  307. cbar.set_label('$V_N$ (nV)')
  308. DIFF = (PRE-VV) / VVS
  309. md = np.max(np.abs(DIFF))
  310. dim = ax3.pcolor(TT, QQQ, DIFF, cmap=cmocean.cm.balance, vmin=-md, vmax=md, shading='auto')
  311. ax3.set_title("misfit / $\widehat{\sigma}$")
  312. cbar2 = plt.colorbar(dim, axc2)
  313. #axc1.set_ylim( [np.min(V), np.max(V)] )
  314. cbar2.outline.set_edgecolor(None)
  315. cbar2.set_label('$V_N$ (nV)')
  316. #plt.colorbar(dim, ax3)
  317. figx.suptitle(ch + " non-linear Inversion")
  318. plt.savefig(ch + "_NLdataspace.pdf")
  319. ich += 1
  320. ###############################################
  321. # Appraise DOI using simplified MRM
  322. ###############################################
  323. CalcDOI = cont['CalcDOI']
  324. if CalcDOI:
  325. pdf = PdfPages('resolution_analysis' + '.pdf' )
  326. MRM = np.zeros((len(ifaces)-1, len(ifaces)-1))
  327. # Build delta models
  328. DELTA = []
  329. for ilay in range(len(ifaces)-1):
  330. #for ilay in range(4):
  331. iDeltaT2 = len(T2Bins)//2
  332. deltaMod = np.zeros( (len(ifaces)-1, len(T2Bins)) )
  333. deltaMod[ilay][iDeltaT2] = 0.3
  334. dV = np.dot(KQT, np.ravel(deltaMod))
  335. #dinv, dibreak, derrn = logBarrier( KQT, dV, T2Bins, "single", MAXITER=1, sigma=np.ravel(VS), alpha=alphastar, smooth="Smallest" )
  336. #output = invertDelta(KQT, dV, T2Bins, np.ravel(VS), alphastar)
  337. DELTA.append(dV)
  338. print("Performing resolution analysis in parallel, printed output may not be inorder.", flush=True)
  339. with multiprocessing.Pool() as pool:
  340. invresults = pool.starmap(invertDelta, zip(itertools.repeat(KQT), DELTA, itertools.repeat(T2Bins), itertools.repeat(np.ravel(VS)), itertools.repeat(alphastar) ))
  341. # invresults = pool.starmap(logBarrier, zip(itertools.repeat(KQT), DELTA, itertools.repeat(T2Bins), itertools.repeat('single'), \
  342. # itertools.repeat('MAXITER=1'), itertools.repeat(np.ravel(VS)), itertools.repeat(alphastar))) #, itertools.repeat(u'smooth=\'Smallest\'')) )
  343. # This could be parallelized
  344. for ilay in range(len(ifaces)-1):
  345. # invert
  346. #dinv, dibreak, derrn = logBarrier(KQT, dV, T2Bins, "single", MAXITER=1, sigma=np.ravel(VS), alpha=alphastar, smooth="Smallest" )
  347. #print("Sum dinv from", str(ifaces[ilay]), "to", str(ifaces[ilay+1]), "=", np.sum(dinv))
  348. dinv, dibreak, derrn = invresults[ilay]
  349. DINV = np.reshape(dinv, (len(ifaces)-1,cont["T2Bins"]["number"]) )
  350. MRM[ilay,:] = np.sum(DINV, axis=1)
  351. Y,X = meshgrid( ifaces, T2Bins2 )
  352. fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  353. ax1 = fig.add_axes( [.2,.15,.6,.7] )
  354. im = ax1.pcolor(X, Y, DINV.T, cmap=cmocean.cm.tempo, shading='auto')
  355. ax1.plot( T2Bins[iDeltaT2], (ifaces[ilay]+ifaces[ilay+1])/2, 's', markersize=6, markeredgecolor='black') #, markerfacecolor=None )
  356. im.set_edgecolor('face')
  357. ax1.set_xlabel(u"$T_2^*$ (ms)")
  358. ax1.set_ylabel(u"depth (m)")
  359. ax1.set_xlim( T2Bins2[0], T2Bins2[-1] )
  360. ax1.set_ylim( ifaces[-1], ifaces[0] )
  361. ax2 = ax1.twiny()
  362. ax2.plot( np.sum(DINV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 , color='mediumblue' )
  363. ax2.set_xlabel(u"total water (m$^3$/m$^3$)", color='mediumblue')
  364. ax2.set_ylim( ifaces[-1], ifaces[0] )
  365. ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )
  366. ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
  367. pdf.savefig(facecolor=[0,0,0,0])
  368. plt.close(fig)
  369. np.save("MRM", MRM)
  370. centres = (ifaces[0:-1]+ifaces[1:])/2
  371. X,Y = np.meshgrid(ifaces,ifaces)
  372. fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  373. ax1 = fig.add_axes( [.2,.15,.6,.7] )
  374. ax1.pcolor(X,Y,MRM, cmap = cmocean.cm.ice)
  375. ax1.set_ylim(ifaces[-1], ifaces[0])
  376. maxDepth = np.argmax(MRM, axis=0)
  377. plt.plot(centres[maxDepth], centres, color='white')
  378. # Determine DOI
  379. DOIMetric = centres[maxDepth]/centres #> 0.9
  380. DOI = ifaces[ np.where(DOIMetric < 0.9 ) ][0]
  381. plt.axhline(y=DOI, color='white', linestyle='-.')
  382. ax1.set_ylim( ifaces[-1], ifaces[0] )
  383. ax1.set_xlim( ifaces[0], ifaces[-1] )
  384. ax1.set_xlabel(u"depth (m)")
  385. ax1.set_ylabel(u"depth (m)")
  386. plt.savefig("resolutionmatrix.pdf")
  387. pdf.close()
  388. INV = np.reshape(inv, (len(ifaces)-1,cont["T2Bins"]["number"]) )
  389. ############## LINEAR RESULT ##########################
  390. Y,X = meshgrid( ifaces, T2Bins2 )
  391. fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  392. ax1 = fig.add_axes( [.2,.15,.6,.7] )
  393. im = ax1.pcolor(X, Y, INV.T, cmap=cmocean.cm.tempo) #cmap='viridis')
  394. im.set_edgecolor('face')
  395. ax1.set_xlim( T2Bins[0], T2Bins[-1] )
  396. ax1.set_ylim( ifaces[-1], ifaces[0] )
  397. cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
  398. cb.locator = MaxNLocator( nbins = 4)
  399. cb.ax.yaxis.set_offset_position('left')
  400. cb.update_ticks()
  401. ax1.set_xlabel(u"$T_2^*$ (ms)")
  402. ax1.set_ylabel(u"depth (m)")
  403. ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
  404. ax1.get_yaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
  405. ax1.xaxis.set_major_locator( MaxNLocator(nbins = 4) )
  406. ax2 = ax1.twiny()
  407. ax2.plot( np.sum(INV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 , color='mediumblue' )
  408. ax2.set_xlabel(u"NMR total water (m$^3$/m$^3$)", color='mediumblue')
  409. ax2.set_ylim( ifaces[-1], ifaces[0] )
  410. ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )
  411. ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
  412. #ax2.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed' )
  413. if CalcDOI:
  414. ax2.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed' )
  415. ax2.tick_params(axis='x', colors='mediumblue')
  416. plt.setp(ax2.get_xticklabels(), color="mediumblue")
  417. plt.savefig("akvoInversion.pdf")
  418. #############
  419. # water plot#
  420. fig2 = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  421. ax = fig2.add_axes( [.2,.15,.6,.7] )
  422. # Bound water cutoff
  423. Bidx = T2Bins<33.0
  424. twater = np.sum(INV, axis=1)
  425. bwater = np.sum(INV[:,Bidx], axis=1)
  426. ax.plot( twater, (ifaces[0:-1]+ifaces[1::])/2, label="NMR total water", color='blue' )
  427. ax.plot( bwater, (ifaces[0:-1]+ifaces[1::])/2, label="NMR bound water", color='green' )
  428. ax.fill_betweenx((ifaces[0:-1]+ifaces[1::])/2 , twater, bwater, where=twater >= bwater, facecolor='blue', alpha=.5)
  429. ax.fill_betweenx((ifaces[0:-1]+ifaces[1::])/2 , bwater, 0, where=bwater >= 0, facecolor='green', alpha=.5)
  430. ax.set_xlabel(r"$\theta_N$ (m$^3$/m$^3$)")
  431. ax.set_ylabel(r"depth (m)")
  432. ax.set_ylim( ifaces[-1], ifaces[0] )
  433. ax.set_xlim( 0, ax.get_xlim()[1] )
  434. #ax.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed' )
  435. if CalcDOI:
  436. ax.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed' )
  437. # Hide the right and top spines
  438. ax.spines['right'].set_visible(False)
  439. ax.spines['top'].set_visible(False)
  440. # Only show ticks on the left and bottom spines
  441. ax.yaxis.set_ticks_position('left')
  442. ax.xaxis.set_ticks_position('bottom')
  443. plt.legend()
  444. plt.savefig("akvoInversionWC.pdf")
  445. fr = pd.DataFrame( INV, columns=T2Bins ) #[0:-1] )
  446. fr.insert(0, "layer top", ifaces[0:-1] )
  447. fr.insert(1, "layer bottom", ifaces[1::] )
  448. fr.insert(2, "NMR total water", np.sum(INV, axis=1) )
  449. fr.insert(3, "NMR bound water", bwater )
  450. fr.insert(4, "Layer SNR", SNR )
  451. if CalcDOI:
  452. fr.insert(5, "Resolution", DOIMetric )
  453. fr.to_csv("akvoInversion.csv", mode='w+')
  454. ############## NONLINEAR RESULT ##########################
  455. if nonLinearRefinement:
  456. Y,X = meshgrid( ifaces, T2Bins2 )
  457. fig = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  458. ax1 = fig.add_axes( [.2,.15,.6,.7] )
  459. im = ax1.pcolor(X, Y, INVc.T, cmap=cmocean.cm.tempo) #cmap='viridis')
  460. #im = ax1.pcolor(X[0:SNRidx,:], Y[0:SNRidx,:], INV.T[0:SNRidx,:], cmap=cmocean.cm.tempo) #cmap='viridis')
  461. #im = ax1.pcolor(X[SNRidx::,:], Y[SNRidx::,:], INV.T[SNRidx::,:], cmap=cmocean.cm.tempo, alpha=.5) #cmap='viridis')
  462. #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
  463. #im = ax1.pcolormesh(X, Y, INV.T, alpha=alphas) #, cmap=cmocean.cm.tempo) #cmap='viridis')
  464. #ax1.axhline( y=ifaces[SNRidx], xmin=T2Bins[0], xmax=T2Bins[-1], color='black' )
  465. im.set_edgecolor('face')
  466. ax1.set_xlim( T2Bins[0], T2Bins[-1] )
  467. ax1.set_ylim( ifaces[-1], ifaces[0] )
  468. cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
  469. cb.locator = MaxNLocator( nbins = 4)
  470. cb.ax.yaxis.set_offset_position('left')
  471. cb.update_ticks()
  472. ax1.set_xlabel(u"$T_2^*$ (ms)")
  473. ax1.set_ylabel(u"depth (m)")
  474. ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
  475. ax1.get_yaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
  476. ax1.xaxis.set_major_locator( MaxNLocator(nbins = 4) )
  477. #ax1.xaxis.set_label_position('top')
  478. ax2 = ax1.twiny()
  479. ax2.plot( np.sum(INVc, axis=1), (ifaces[1:]+ifaces[0:-1])/2 , color='mediumblue' )
  480. ax2.set_xlabel(u"NMR total water (m$^3$/m$^3$)", color='mediumblue')
  481. ax2.set_ylim( ifaces[-1], ifaces[0] )
  482. ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )
  483. ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
  484. #ax2.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed' )
  485. if CalcDOI:
  486. ax2.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed' )
  487. #ax2.xaxis.set_label_position('bottom')
  488. #fig.suptitle("Non linear inversion")
  489. ax2.tick_params(axis='x', colors='mediumblue')
  490. plt.setp(ax2.get_xticklabels(), color="mediumblue")
  491. plt.savefig("akvoInversionNL.pdf")
  492. #############
  493. # water plot#
  494. fig2 = plt.figure( figsize=(pc2in(20.0),pc2in(22.)) )
  495. ax = fig2.add_axes( [.2,.15,.6,.7] )
  496. # Bound water cutoff
  497. Bidx = T2Bins<33.0
  498. twater = np.sum(INVc, axis=1)
  499. bwater = np.sum(INVc[:,Bidx], axis=1)
  500. ax.plot( twater, (ifaces[0:-1]+ifaces[1::])/2, label="NMR total water", color='blue' )
  501. ax.plot( bwater, (ifaces[0:-1]+ifaces[1::])/2, label="NMR bound water", color='green' )
  502. ax.fill_betweenx((ifaces[0:-1]+ifaces[1::])/2 , twater, bwater, where=twater >= bwater, facecolor='blue', alpha=.5)
  503. ax.fill_betweenx((ifaces[0:-1]+ifaces[1::])/2 , bwater, 0, where=bwater >= 0, facecolor='green', alpha=.5)
  504. ax.set_xlabel(r"$\theta_N$ (m$^3$/m$^3$)")
  505. ax.set_ylabel(r"depth (m)")
  506. ax.set_ylim( ifaces[-1], ifaces[0] )
  507. ax.set_xlim( 0, ax.get_xlim()[1] )
  508. #ax.axhline( y=ifaces[SNRidx], xmin=0, xmax=1, color='black', linestyle='dashed' )
  509. if CalcDOI:
  510. ax.axhline( y=DOI, xmin=0, xmax=1, color='black', linestyle='dashed' )
  511. # Hide the right and top spines
  512. ax.spines['right'].set_visible(False)
  513. ax.spines['top'].set_visible(False)
  514. # Only show ticks on the left and bottom spines
  515. ax.yaxis.set_ticks_position('left')
  516. ax.xaxis.set_ticks_position('bottom')
  517. plt.savefig("akvoNLInversionWC.pdf")
  518. plt.legend()
  519. # Report results into a text file
  520. fr = pd.DataFrame( INVc, columns=T2Bins ) #[0:-1] )
  521. fr.insert(0, "layer top", ifaces[0:-1] )
  522. fr.insert(1, "layer bottom", ifaces[1::] )
  523. fr.insert(2, "NMR total water", np.sum(INVc, axis=1) )
  524. fr.insert(3, "NMR bound water", bwater )
  525. fr.insert(4, "Layer SNR", SNR )
  526. if CalcDOI:
  527. fr.insert(5, "Resolution", DOIMetric )
  528. fr.to_csv("akvoNLInversion.csv", mode='w+')
  529. #fr.to_excel("akvoInversion.xlsx")
  530. plt.show()
  531. if __name__ == "__main__":
  532. main()