Browse Source

Improvements to plotting of pulse moment calculation.

tags/1.6.1
Trevor Irons 5 years ago
parent
commit
719b0ba530
2 changed files with 78 additions and 20 deletions
  1. 54
    9
      akvo/terminal/plotyaml.py
  2. 24
    11
      akvo/tressel/mrsurvey.py

+ 54
- 9
akvo/terminal/plotyaml.py View File

1
-import yaml
1
+from ruamel import yaml
2
 import os, sys
2
 import os, sys
3
 import numpy as np
3
 import numpy as np
4
 
4
 
5
+import matplotlib.pyplot as plt
6
+import seaborn as sns
7
+sns.set(style="ticks")
8
+import cmocean 
9
+from SEGPlot import *
10
+from matplotlib.ticker import FormatStrFormatter
11
+import matplotlib.ticker as plticker
12
+
5
 def slicedict(d, s):
13
 def slicedict(d, s):
6
     return {k:v for k,v in d.items() if k.startswith(s)}
14
     return {k:v for k,v in d.items() if k.startswith(s)}
7
 
15
 
31
         # Converts to a dictionary with Eigen vectors represented as Numpy arrays 
39
         # Converts to a dictionary with Eigen vectors represented as Numpy arrays 
32
         return self
40
         return self
33
 
41
 
34
-def loadAkvoData(fnamein, chan):
42
+def loadAkvoData(fnamein):
35
     """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be 
43
     """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be 
36
         corrected in future kernel calculations. The current was reported but not the pulse length. 
44
         corrected in future kernel calculations. The current was reported but not the pulse length. 
37
     """
45
     """
38
     fname = (os.path.splitext(fnamein)[0])
46
     fname = (os.path.splitext(fnamein)[0])
39
     with open(fnamein, 'r') as stream:
47
     with open(fnamein, 'r') as stream:
40
         try:
48
         try:
41
-            AKVO = (yaml.load(stream))
49
+            AKVO = (yaml.load(stream, Loader=yaml.Loader))
42
         except yaml.YAMLError as exc:
50
         except yaml.YAMLError as exc:
43
             print(exc)
51
             print(exc)
44
     return AKVO 
52
     return AKVO 
45
 
53
 
46
 def plotQt( akvo ):
54
 def plotQt( akvo ):
47
-    import matplotlib.pyplot as plt
55
+ 
48
     plt.style.use('ggplot')
56
     plt.style.use('ggplot')
57
+    #plt.style.use('seaborn-white')
49
     for pulse in akvo.Gated:
58
     for pulse in akvo.Gated:
50
         if pulse[0:5] == "Pulse":
59
         if pulse[0:5] == "Pulse":
51
             #print(akvo.GATED[pulse].keys())
60
             #print(akvo.GATED[pulse].keys())
52
             nq = akvo.Pulses[pulse]["current"].size
61
             nq = akvo.Pulses[pulse]["current"].size
53
             for chan in slicedict(akvo.Gated[pulse], "Chan."):
62
             for chan in slicedict(akvo.Gated[pulse], "Chan."):
54
                 # accumulate pulse moments
63
                 # accumulate pulse moments
55
-                X = np.zeros( (nq, len( akvo.Gated[pulse]["abscissa"].data )) )
64
+                CA = np.zeros( (nq, len( akvo.Gated[pulse]["abscissa"].data )) )
65
+                RE = np.zeros( (nq, len( akvo.Gated[pulse]["abscissa"].data )) )
66
+                IM = np.zeros( (nq, len( akvo.Gated[pulse]["abscissa"].data )) )
56
                 for q in range(nq):
67
                 for q in range(nq):
57
-                    plt.plot(  akvo.Gated[pulse]["abscissa"].data,  akvo.Gated[pulse][chan]["Q-" + str(q)+" CA"].data ) 
58
-                    X[q] = akvo.Gated[pulse][chan]["Q-" + str(q)+" CA"].data
59
-    plt.matshow(X)
68
+                    #plt.plot(  akvo.Gated[pulse]["abscissa"].data,  akvo.Gated[pulse][chan]["Q-" + str(q)+" CA"].data ) 
69
+                    CA[q] = akvo.Gated[pulse][chan]["Q-" + str(q)+" CA"].data
70
+                    RE[q] = akvo.Gated[pulse][chan]["Q-" + str(q)+" RE"].data
71
+                    IM[q] = akvo.Gated[pulse][chan]["Q-" + str(q)+" IM"].data
72
+                    #X[q] = akvo.Gated[pulse][chan]["Q-" + str(q)+" RE"].data
73
+            Windows = akvo.Gated[pulse]["windows"].data
74
+            Q = np.array(akvo.Pulses[pulse]["current"].data)
75
+            print("pulse length ", akvo.pulseLength[0])
76
+            Q *= akvo.pulseLength[0]
77
+       
78
+            fig = plt.figure( figsize=( pc2in(20), pc2in(26) ) ) 
79
+            ax1 = fig.add_axes([.25,.05,.6,.9])
80
+            im = ax1.pcolormesh(Windows,Q,CA, cmap=cmocean.cm.curl_r, vmin=-np.max(np.abs(CA)), vmax=(np.max(np.abs(CA))))
81
+            cb = plt.colorbar( im, orientation='horizontal', pad=.175,  )
82
+            cb.set_label("FID (nV)", fontsize=10)
83
+            cb.ax.tick_params(labelsize=10) 
84
+            ax1.set_yscale('log')
85
+            ax1.set_ylabel("Q (A$\cdot$s)", fontsize=10)
86
+            ax1.set_xscale('log')
87
+            ax1.set_xlabel("time (s)", fontsize=10)
88
+
89
+            #loc = plticker.MultipleLocator(25) #base=50.0) # this locator puts ticks at regular intervals
90
+            #loc = plticker.MaxNLocator(5, steps=[1,2]) #base=50.0) # this locator puts ticks at regular intervals
91
+            #ax1.xaxis.set_major_locator(loc)
92
+
93
+            #ax1.xaxis.set_minor_locator(plticker.NullLocator())
94
+            ax1.xaxis.set_major_formatter(FormatStrFormatter('%2.0f'))
95
+            ax1.yaxis.set_major_formatter(FormatStrFormatter('%2.1f'))
96
+
97
+            #plt.figure()
98
+            #sns.kdeplot(Windows, Q, CA) #, kind="hex", color="#4CB391")
99
+            #sns.heatmap(CA, annot=False, center=0)
100
+
101
+    #plt.matshow(RE)
102
+    #plt.matshow(IM)
103
+
104
+
60
 
105
 
61
     plt.show()
106
     plt.show()
62
 if __name__ == "__main__":
107
 if __name__ == "__main__":
63
-    akvo = loadAkvoData( sys.argv[1] , "Chan. 1")
108
+    akvo = loadAkvoData( sys.argv[1] ) #, "Chan. 1")
64
     plotQt(akvo)
109
     plotQt(akvo)

+ 24
- 11
akvo/tressel/mrsurvey.py View File

363
                 ax1.set_xlabel(r"time (ms)", fontsize=8)
363
                 ax1.set_xlabel(r"time (ms)", fontsize=8)
364
 
364
 
365
                 if ichan == 0:
365
                 if ichan == 0:
366
-                    ax1.set_ylabel(r"signal [nV]", fontsize=8)
366
+                    ax1.set_ylabel(r"signal (nV)", fontsize=8)
367
                 else:
367
                 else:
368
                     plt.setp(ax1.get_yticklabels(), visible=False)
368
                     plt.setp(ax1.get_yticklabels(), visible=False)
369
                     plt.setp(ax1.get_yaxis().get_offset_text(), visible=False) 
369
                     plt.setp(ax1.get_yaxis().get_offset_text(), visible=False) 
499
                 #im2.append(ax2.matshow( db, aspect='auto', vmin=vvmin, vmax=vvmax))
499
                 #im2.append(ax2.matshow( db, aspect='auto', vmin=vvmin, vmax=vvmax))
500
                 #im2 = ax2.matshow( db, aspect='auto', cmap=cmocean.cm.ice, vmin=vvmin, vmax=vvmax)
500
                 #im2 = ax2.matshow( db, aspect='auto', cmap=cmocean.cm.ice, vmin=vvmin, vmax=vvmax)
501
                 if ichan == 0:
501
                 if ichan == 0:
502
-                    ax2.set_ylabel(r"$q$ (A $\cdot$ s)", fontsize=8)
502
+                    #ax2.set_ylabel(r"$q$ (A $\cdot$ s)", fontsize=8)
503
+                    ax2.set_ylabel(r"pulse index", fontsize=8)
504
+                    #ax1.set_ylabel(r"FID (nV)", fontsize=8)
503
                 else:
505
                 else:
504
                     #ax2.yaxis.set_ticklabels([])
506
                     #ax2.yaxis.set_ticklabels([])
505
                     plt.setp(ax2.get_yticklabels(), visible=False)
507
                     plt.setp(ax2.get_yticklabels(), visible=False)
1084
                 ax1 = axes[2*ichan  ]
1086
                 ax1 = axes[2*ichan  ]
1085
                 ax2 = axes[2*ichan+1] # TODO fix hard coded number
1087
                 ax2 = axes[2*ichan+1] # TODO fix hard coded number
1086
                 if phase == 0: # Re Im 
1088
                 if phase == 0: # Re Im 
1087
-                    print("plot dog", np.shape(QQ), np.shape(self.DATADICT["RE"][pulse][chan]))
1088
-                    print("QQ", QQ) 
1089
+                    #print("plot dog", np.shape(QQ), np.shape(self.DATADICT["RE"][pulse][chan]))
1090
+                    #print("QQ", QQ) 
1089
                     im1 = ax1.pcolormesh( time_sp, QQ[iQ], self.DATADICT["RE"][pulse][chan][iQ], cmap=dcmap, \
1091
                     im1 = ax1.pcolormesh( time_sp, QQ[iQ], self.DATADICT["RE"][pulse][chan][iQ], cmap=dcmap, \
1090
                          vmin=-self.DATADICT["REmax"][pulse] , vmax=self.DATADICT["REmax"][pulse] )
1092
                          vmin=-self.DATADICT["REmax"][pulse] , vmax=self.DATADICT["REmax"][pulse] )
1091
                     im2 = ax2.pcolormesh( time_sp, QQ[iQ], self.DATADICT["IM"][pulse][chan][iQ], cmap=dcmap, \
1093
                     im2 = ax2.pcolormesh( time_sp, QQ[iQ], self.DATADICT["IM"][pulse][chan][iQ], cmap=dcmap, \
1575
         
1577
         
1576
         canvas.reAxH(2)
1578
         canvas.reAxH(2)
1577
         nstack = len(self.DATADICT["stacks"]) 
1579
         nstack = len(self.DATADICT["stacks"]) 
1578
-        canvas.ax1.set_yscale('log')
1580
+        #canvas.ax1.set_yscale('log')
1579
 
1581
 
1580
         for pulse in self.DATADICT["PULSES"]:
1582
         for pulse in self.DATADICT["PULSES"]:
1581
             self.DATADICT[pulse]["qeff"] = {}
1583
             self.DATADICT[pulse]["qeff"] = {}
1600
                     v = np.fft.fftfreq(len(x), self.dt)
1602
                     v = np.fft.fftfreq(len(x), self.dt)
1601
                     v = v[0:len(X)]
1603
                     v = v[0:len(X)]
1602
                     v[-1] = np.abs(v[-1])
1604
                     v[-1] = np.abs(v[-1])
1605
+
1603
                     # calculate effective current/moment
1606
                     # calculate effective current/moment
1604
                     I0 = np.abs(X)/len(X) 
1607
                     I0 = np.abs(X)/len(X) 
1605
-                    qeff = I0 * (self.DATADICT[pulse]["PULSE_TIMES"][-1]-self.DATADICT[pulse]["PULSE_TIMES"][0])
1606
-                    canvas.ax1.set_title(r"pulse moment index " +str(ipm), fontsize=10)
1607
-                    canvas.ax1.set_xlabel(r"$\nu$ [Hz]", fontsize=8)
1608
-                    canvas.ax1.set_ylabel(r"$q_{eff}$ [A$\cdot$sec]", fontsize=8)
1609
-                    canvas.ax1.plot(v, qeff, color=scolours[iistack] ) # eff current
1608
+                    qeff = I0 #* (self.DATADICT[pulse]["PULSE_TIMES"][-1]-self.DATADICT[pulse]["PULSE_TIMES"][0])
1609
+
1610
+                    # frequency plot
1611
+                    #canvas.ax1.set_title(r"pulse moment index " +str(ipm), fontsize=10)
1612
+                    #canvas.ax1.set_xlabel(r"$\nu$ [Hz]", fontsize=8)
1613
+                    #canvas.ax1.set_ylabel(r"$q_{eff}$ [A$\cdot$sec]", fontsize=8)
1614
+                    #canvas.ax1.plot(v, qeff, color=scolours[iistack] ) # eff current
1615
+                    
1616
+                    # time plot
1617
+                    canvas.ax1.plot(1e2*(self.DATADICT[pulse]["PULSE_TIMES"]-self.DATADICT[pulse]["PULSE_TIMES"][0]), x, color=scolours[iistack])
1618
+
1610
                     self.DATADICT[pulse]["qeff"][ipm][istack] = qeff
1619
                     self.DATADICT[pulse]["qeff"][ipm][istack] = qeff
1611
                     self.DATADICT[pulse]["q_nu"][ipm][istack] = v
1620
                     self.DATADICT[pulse]["q_nu"][ipm][istack] = v
1612
                     iistack += 1
1621
                     iistack += 1
1622
+
1623
+                canvas.ax1.set_xlabel("time (ms)", fontsize=8)
1624
+                canvas.ax1.set_ylabel("current (A)", fontsize=8)
1613
                 canvas.draw()
1625
                 canvas.draw()
1626
+
1614
                         
1627
                         
1615
                 percent = int(1e2* (float)((istack)+ipm*self.DATADICT["nPulseMoments"]) / 
1628
                 percent = int(1e2* (float)((istack)+ipm*self.DATADICT["nPulseMoments"]) / 
1616
                                    (float)(len(self.DATADICT["PULSES"])*self.DATADICT["nPulseMoments"]*nstack))
1629
                                    (float)(len(self.DATADICT["PULSES"])*self.DATADICT["nPulseMoments"]*nstack))
1640
                 istack = 0
1653
                 istack = 0
1641
                 for stack in self.DATADICT["stacks"]:
1654
                 for stack in self.DATADICT["stacks"]:
1642
                     # find index 
1655
                     # find index 
1643
-                    icv = int (round(cv / self.DATADICT[pulse]["q_nu"][ipm][stack][1]))
1656
+                    icv = int(round(cv / self.DATADICT[pulse]["q_nu"][ipm][stack][1]))
1644
                     self.DATADICT[pulse]["Q"][ipm,istack] = self.DATADICT[pulse]["qeff"][ipm][stack][icv]
1657
                     self.DATADICT[pulse]["Q"][ipm,istack] = self.DATADICT[pulse]["qeff"][ipm][stack][icv]
1645
                     if ilabel:
1658
                     if ilabel:
1646
                         ax.scatter(ipm, self.DATADICT[pulse]["qeff"][ipm][stack][icv], facecolors='none', edgecolors=scolours[istack], label=(str(pulse)))
1659
                         ax.scatter(ipm, self.DATADICT[pulse]["qeff"][ipm][stack][icv], facecolors='none', edgecolors=scolours[istack], label=(str(pulse)))

Loading…
Cancel
Save