##// END OF EJS Templates
Se modifico la forma de calcular la phase de la coherencia
Miguel Valdez -
r264:338ee2ca09fe
parent child
Show More
@@ -2,17 +2,16 import numpy
2 2 import datetime
3 3 import sys
4 4 import matplotlib
5 if sys.platform == 'linux':
6 matplotlib.use("GTKAgg")
7 if sys.platform == 'darwin':
5
6 if 'linux' in sys.platform:
7 matplotlib.use("TKAgg")
8
9 if 'darwin' in sys.platform:
8 10 matplotlib.use("TKAgg")
11
9 12 import matplotlib.pyplot
10 import matplotlib.dates
11 #import scitools.numpyutils
12 from mpl_toolkits.axes_grid1 import make_axes_locatable
13 13
14 from matplotlib.dates import DayLocator, HourLocator, MinuteLocator, SecondLocator, DateFormatter
15 from matplotlib.ticker import FuncFormatter
14 from mpl_toolkits.axes_grid1 import make_axes_locatable
16 15 from matplotlib.ticker import *
17 16
18 17 ###########################################
@@ -360,10 +360,10 class JRODataReader(JRODataIO, ProcessingUnit):
360 360 thisDate += datetime.timedelta(1)
361 361
362 362 if pathList == []:
363 print "Any folder found into date range %s-%s" %(startDate, endDate)
363 print "Any folder was found for the date range: %s-%s" %(startDate, endDate)
364 364 return None, None
365 365
366 print "%d folder(s) found [%s, ...]" %(len(pathList), pathList[0])
366 print "%d folder(s) was(were) found for the date range: %s-%s" %(len(pathList), startDate, endDate)
367 367
368 368 filenameList = []
369 369 for thisPath in pathList:
@@ -379,9 +379,11 class JRODataReader(JRODataIO, ProcessingUnit):
379 379 filenameList.append(filename)
380 380
381 381 if not(filenameList):
382 print "Any file found into time range %s-%s" %(startTime, endTime)
382 print "Any file was found for the time range %s - %s" %(startTime, endTime)
383 383 return None, None
384 384
385 print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime)
386
385 387 self.filenameList = filenameList
386 388
387 389 return pathList, filenameList
@@ -2077,10 +2079,19 class SpectraReader(JRODataReader):
2077 2079
2078 2080 self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
2079 2081
2080 self.dataOut.flagDecodeData = True #asumo q la data no esta decodificada
2082 self.dataOut.flagDecodeData = False #asumo q la data no esta decodificada
2081 2083
2082 2084 self.dataOut.flagDeflipData = True #asumo q la data no esta sin flip
2083 2085
2086 if self.processingHeaderObj.code != None:
2087
2088 self.dataOut.nCode = self.processingHeaderObj.nCode
2089
2090 self.dataOut.nBaud = self.processingHeaderObj.nBaud
2091
2092 self.dataOut.code = self.processingHeaderObj.code
2093
2094 self.dataOut.flagDecodeData = True
2084 2095
2085 2096 return self.dataOut.data_spc
2086 2097
@@ -54,7 +54,8 class CrossSpectraPlot(Figure):
54 54
55 55 def run(self, dataOut, idfigure, wintitle="", pairsList=None, showprofile='True',
56 56 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
57 save=False, figpath='./', figfile=None):
57 save=False, figpath='./', figfile=None,
58 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r'):
58 59
59 60 """
60 61
@@ -90,8 +91,8 class CrossSpectraPlot(Figure):
90 91 x = dataOut.getVelRange(1)
91 92 y = dataOut.getHeiRange()
92 93 z = dataOut.data_spc[:,:,:]/factor
93 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
94 avg = numpy.average(numpy.abs(z), axis=1)
94 # z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
95 avg = numpy.abs(numpy.average(z, axis=1))
95 96 noise = dataOut.getNoise()/factor
96 97
97 98 zdB = 10*numpy.log10(z)
@@ -133,7 +134,7 class CrossSpectraPlot(Figure):
133 134 axes0.pcolor(x, y, zdB,
134 135 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
135 136 xlabel=xlabel, ylabel=ylabel, title=title,
136 ticksize=9, cblabel='')
137 ticksize=9, colormap=power_cmap, cblabel='')
137 138
138 139 title = "Channel %d: %4.2fdB" %(pair[1], noisedB[pair[1]])
139 140 zdB = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:]/factor)
@@ -141,26 +142,26 class CrossSpectraPlot(Figure):
141 142 axes0.pcolor(x, y, zdB,
142 143 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
143 144 xlabel=xlabel, ylabel=ylabel, title=title,
144 ticksize=9, cblabel='')
145 ticksize=9, colormap=power_cmap, cblabel='')
145 146
146 147 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
147 148 coherence = numpy.abs(coherenceComplex)
148 phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
149
149 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
150 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
150 151
151 152 title = "Coherence %d%d" %(pair[0], pair[1])
152 153 axes0 = self.axesList[i*self.__nsubplots+2]
153 154 axes0.pcolor(x, y, coherence,
154 155 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1,
155 156 xlabel=xlabel, ylabel=ylabel, title=title,
156 ticksize=9, cblabel='')
157 ticksize=9, colormap=coherence_cmap, cblabel='')
157 158
158 159 title = "Phase %d%d" %(pair[0], pair[1])
159 160 axes0 = self.axesList[i*self.__nsubplots+3]
160 161 axes0.pcolor(x, y, phase,
161 162 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
162 163 xlabel=xlabel, ylabel=ylabel, title=title,
163 ticksize=9, cblabel='', colormap='RdBu_r')
164 ticksize=9, colormap=phase_cmap, cblabel='')
164 165
165 166
166 167
@@ -719,7 +720,7 class CoherenceMap(Figure):
719 720
720 721 WIDTHPROF = None
721 722 HEIGHTPROF = None
722 PREFIX = 'coherencemap'
723 PREFIX = 'cmap'
723 724
724 725 def __init__(self):
725 726 self.timerange = 2*60*60
@@ -766,7 +767,8 class CoherenceMap(Figure):
766 767 def run(self, dataOut, idfigure, wintitle="", pairsList=None, showprofile='True',
767 768 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
768 769 timerange=None,
769 save=False, figpath='./', figfile=None):
770 save=False, figpath='./', figfile=None,
771 coherence_cmap='jet', phase_cmap='RdBu_r'):
770 772
771 773 if pairsList == None:
772 774 pairsIndexList = dataOut.pairsIndexList
@@ -817,9 +819,12 class CoherenceMap(Figure):
817 819
818 820 pair = dataOut.pairsList[pairsIndexList[i]]
819 821 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
820 coherence = numpy.abs(coherenceComplex)
821 avg = numpy.average(coherence, axis=0)
822 z = avg.reshape((1,-1))
822 avgcoherenceComplex = numpy.average(coherenceComplex, axis=0)
823 coherence = numpy.abs(avgcoherenceComplex)
824 # coherence = numpy.abs(coherenceComplex)
825 # avg = numpy.average(coherence, axis=0)
826
827 z = coherence.reshape((1,-1))
823 828
824 829 counter = 0
825 830
@@ -828,33 +833,34 class CoherenceMap(Figure):
828 833 axes.pcolor(x, y, z,
829 834 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1,
830 835 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
831 ticksize=9, cblabel='', cbsize="1%")
836 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
832 837
833 838 if self.__showprofile:
834 839 counter += 1
835 840 axes = self.axesList[i*self.__nsubplots*2 + counter]
836 axes.pline(avg, y,
841 axes.pline(coherence, y,
837 842 xmin=0, xmax=1, ymin=ymin, ymax=ymax,
838 843 xlabel='', ylabel='', title='', ticksize=7,
839 844 ytick_visible=False, nxticks=5,
840 845 grid='x')
841 846
842 847 counter += 1
843 phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
844 avg = numpy.average(phase, axis=0)
845 z = avg.reshape((1,-1))
848 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
849 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
850 # avg = numpy.average(phase, axis=0)
851 z = phase.reshape((1,-1))
846 852
847 853 title = "Phase %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
848 854 axes = self.axesList[i*self.__nsubplots*2 + counter]
849 855 axes.pcolor(x, y, z,
850 856 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
851 857 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
852 ticksize=9, cblabel='', colormap='RdBu', cbsize="1%")
858 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
853 859
854 860 if self.__showprofile:
855 861 counter += 1
856 862 axes = self.axesList[i*self.__nsubplots*2 + counter]
857 axes.pline(avg, y,
863 axes.pline(phase, y,
858 864 xmin=-180, xmax=180, ymin=ymin, ymax=ymax,
859 865 xlabel='', ylabel='', title='', ticksize=7,
860 866 ytick_visible=False, nxticks=4,
General Comments 0
You need to be logged in to leave comments. Login now