##// END OF EJS Templates
schain
RSS
cambios para amisr ISR
joabAM -
r1465:65c0d2b45bc1
parent child
Show More
Side by Side Unified
Context file:
r1465:65c0d2b45bc1 Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,132
1 import numpy
2
3 def setConstants(dataOut):
4 dictionary = {}
5 dictionary["M"] = dataOut.normFactor
6 dictionary["N"] = dataOut.nFFTPoints
7 dictionary["ippSeconds"] = dataOut.ippSeconds
8 dictionary["K"] = dataOut.nIncohInt
9
10 return dictionary
11
12 def initialValuesFunction(data_spc, constants):
13 #Constants
14 M = constants["M"]
15 N = constants["N"]
16 ippSeconds = constants["ippSeconds"]
17
18 S1 = data_spc[0,:]/(N*M)
19 S2 = data_spc[1,:]/(N*M)
20
21 Nl=min(S1)
22 A=sum(S1-Nl)/N
23 #x = dataOut.getVelRange() #below matches Madrigal data better
24 x=numpy.linspace(-(N/2)/(N*ippSeconds),(N/2-1)/(N*ippSeconds),N)*-(6.0/2)
25 v=sum(x*(S1-Nl))/sum(S1-Nl)
26 al1=numpy.sqrt(sum(x**2*(S1-Nl))/sum(S2-Nl)-v**2)
27 p0=[al1,A,A,v,min(S1),min(S2)]#first guess(width,amplitude,velocity,noise)
28 return p0
29
30 def modelFunction(p, constants):
31 ippSeconds = constants["ippSeconds"]
32 N = constants["N"]
33
34 fm_c = ACFtoSPC(p, constants)
35 fm = numpy.hstack((fm_c[0],fm_c[1]))
36 return fm
37
38 def errorFunction(p, constants, LT):
39
40 J=makeJacobian(p, constants)
41 J =numpy.dot(LT,J)
42 covm =numpy.linalg.inv(numpy.dot(J.T ,J))
43 #calculate error as the square root of the covariance matrix diagonal
44 #multiplying by 1.96 would give 95% confidence interval
45 err =numpy.sqrt(numpy.diag(covm))
46 return err
47
48 #-----------------------------------------------------------------------------------
49
50 def ACFw(alpha,A1,A2,vd,x,N,ippSeconds):
51 #creates weighted autocorrelation function based on the operational model
52 #x is n or N-n
53 k=2*numpy.pi/3.0
54 pdt=x*ippSeconds
55 #both correlated channels ACFs are created at the sametime
56 R1=A1*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2)
57 R2=A2*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2)
58 # T is the triangle weigthing function
59 T=1-abs(x)/N
60 Rp1=T*R1
61 Rp2=T*R2
62 return [Rp1,Rp2]
63
64 def ACFtoSPC(p, constants):
65 #calls the create ACF function and transforms the ACF to spectra
66 N = constants["N"]
67 ippSeconds = constants["ippSeconds"]
68
69 n=numpy.linspace(0,(N-1),N)
70 Nn=N-n
71 R = ACFw(p[0],p[1],p[2],p[3],n,N,ippSeconds)
72 RN = ACFw(p[0],p[1],p[2],p[3],Nn,N,ippSeconds)
73 Rf1=R[0]+numpy.conjugate(RN[0])
74 Rf2=R[1]+numpy.conjugate(RN[1])
75 sw1=numpy.fft.fft(Rf1,n=N)
76 sw2=numpy.fft.fft(Rf2,n=N)
77 #the fft needs to be shifted, noise added, and takes only the real part
78 sw0=numpy.real(numpy.fft.fftshift(sw1))+abs(p[4])
79 sw1=numpy.real(numpy.fft.fftshift(sw2))+abs(p[5])
80 return [sw0,sw1]
81
82 def makeJacobian(p, constants):
83 #create Jacobian matrix
84 N = constants["N"]
85 IPPt = constants["ippSeconds"]
86
87 n=numpy.linspace(0,(N-1),N)
88 Nn=N-n
89 k=2*numpy.pi/3.0
90 #created weighted ACF
91 R=ACFw(p[0],p[1],p[2],p[3],n,N,IPPt)
92 RN=ACFw(p[0],p[1],p[2],p[3],Nn,N,IPPt)
93 #take derivatives with respect to the fit parameters
94 Jalpha1=R[0]*-1*(k*n*IPPt)**2*p[0]/(1+(p[0]*k*n*IPPt)**2)+numpy.conjugate(RN[0]*-1*(k*Nn*IPPt)**2*p[0]/(1+(p[0]*k*Nn*IPPt)**2))
95 Jalpha2=R[1]*-1*(k*n*IPPt)**2*p[0]/(1+(p[0]*k*n*IPPt)**2)+numpy.conjugate(RN[1]*-1*(k*Nn*IPPt)**2*p[0]/(1+(p[0]*k*Nn*IPPt)**2))
96 JA1=R[0]/p[1]+numpy.conjugate(RN[0]/p[1])
97 JA2=R[1]/p[2]+numpy.conjugate(RN[1]/p[2])
98 Jvd1=R[0]*-1j*k*n*IPPt+numpy.conjugate(RN[0]*-1j*k*Nn*IPPt)
99 Jvd2=R[1]*-1j*k*n*IPPt+numpy.conjugate(RN[1]*-1j*k*Nn*IPPt)
100 #fft
101 sJalp1=numpy.fft.fft(Jalpha1,n=N)
102 sJalp2=numpy.fft.fft(Jalpha2,n=N)
103 sJA1=numpy.fft.fft(JA1,n=N)
104 sJA2=numpy.fft.fft(JA2,n=N)
105 sJvd1=numpy.fft.fft(Jvd1,n=N)
106 sJvd2=numpy.fft.fft(Jvd2,n=N)
107 sJalp1=numpy.real(numpy.fft.fftshift(sJalp1))
108 sJalp2=numpy.real(numpy.fft.fftshift(sJalp2))
109 sJA1=numpy.real(numpy.fft.fftshift(sJA1))
110 sJA2=numpy.real(numpy.fft.fftshift(sJA2))
111 sJvd1=numpy.real(numpy.fft.fftshift(sJvd1))
112 sJvd2=numpy.real(numpy.fft.fftshift(sJvd2))
113 sJnoise=numpy.ones(numpy.shape(sJvd1))
114 #combine arrays
115 za=numpy.zeros([N])
116 sJalp=zip(sJalp1,sJalp2)
117 sJA1=zip(sJA1,za)
118 sJA2=zip(za,sJA2)
119 sJvd=zip(sJvd1,sJvd2)
120 sJn1=zip(sJnoise, za)
121 sJn2=zip(za, sJnoise)
122 #reshape from 2D to 1D
123 sJalp=numpy.reshape(list(sJalp), [2*N])
124 sJA1=numpy.reshape(list(sJA1), [2*N])
125 sJA2=numpy.reshape(list(sJA2), [2*N])
126 sJvd=numpy.reshape(list(sJvd), [2*N])
127 sJn1=numpy.reshape(list(sJn1), [2*N])
128 sJn2=numpy.reshape(list(sJn2), [2*N])
129 #combine into matrix and transpose
130 J=numpy.array([sJalp,sJA1,sJA2,sJvd,sJn1,sJn2])
131 J=J.T
132 return J
Show file before | Show file after | Hide comments
@@ -257,6 +257,8 class Plot(Operation):
257 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
257 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
258 self.tmin = kwargs.get('tmin', None)
258 self.tmin = kwargs.get('tmin', None)
259 self.t_units = kwargs.get('t_units', "h_m")
259 self.t_units = kwargs.get('t_units', "h_m")
260 self.selectedHeight = kwargs.get('selectedHeight', None)
261
260
262
261 if self.server:
263 if self.server:
262 if not self.server.startswith('tcp://'):
264 if not self.server.startswith('tcp://'):
Show file before | Show file after | Hide comments
@@ -14,7 +14,7 import matplotlib.pyplot as plt
14 class SpectraHeisPlot(Plot):
14 class SpectraHeisPlot(Plot):
15
15
16 CODE = 'spc_heis'
16 CODE = 'spc_heis'
17
17 channelList = []
18 def setup(self):
18 def setup(self):
19
19
20 self.nplots = len(self.data.channels)
20 self.nplots = len(self.data.channels)
@@ -28,6 +28,8 class SpectraHeisPlot(Plot):
28 self.colorbar = False
28 self.colorbar = False
29
29
30 def update(self, dataOut):
30 def update(self, dataOut):
31 if len(self.channelList) == 0:
32 self.channelList = dataOut.channelList
31
33
32 data = {}
34 data = {}
33 meta = {}
35 meta = {}
@@ -51,9 +53,12 class SpectraHeisPlot(Plot):
51 self.xmin = min(x) if self.xmin is None else self.xmin
53 self.xmin = min(x) if self.xmin is None else self.xmin
52 self.xmax = max(x) if self.xmax is None else self.xmax
54 self.xmax = max(x) if self.xmax is None else self.xmax
53 ax.plt = ax.plot(x, ychannel, lw=1, color='b')[0]
55 ax.plt = ax.plot(x, ychannel, lw=1, color='b')[0]
56 ax.set_ylim(ymin=self.zmin, ymax=self.zmax)
57 ax.set_xlim(xmin=self.xmin, xmax=self.xmax)
54 else:
58 else:
55 ax.plt.set_data(x, ychannel)
59 ax.plt.set_data(x, ychannel)
56
60 ax.set_ylim(ymin=self.zmin, ymax=self.zmax)
61 ax.set_xlim(xmin=self.xmin, xmax=self.xmax)
57 self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel)))
62 self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel)))
58 plt.suptitle(Maintitle)
63 plt.suptitle(Maintitle)
59
64
Show file before | Show file after | Hide comments
@@ -48,12 +48,13 class SnrPlot(RTIPlot):
48 colormap = 'jet'
48 colormap = 'jet'
49
49
50 def update(self, dataOut):
50 def update(self, dataOut):
51 self.update_list(dataOut)
51 if len(self.channelList) == 0:
52 data = {
52 self.update_list(dataOut)
53 'snr': 10*numpy.log10(dataOut.data_snr)
53 data = {}
54 }
54 meta = {}
55
55 data['snr'] = 10*numpy.log10(dataOut.data_snr)
56 return data, {}
56
57 return data, meta
57
58
58 class DopplerPlot(RTIPlot):
59 class DopplerPlot(RTIPlot):
59 '''
60 '''
@@ -84,6 +85,7 class PowerPlot(RTIPlot):
84 data = {
85 data = {
85 'pow': 10*numpy.log10(dataOut.data_pow)
86 'pow': 10*numpy.log10(dataOut.data_pow)
86 }
87 }
88 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
87 return data, {}
89 return data, {}
88
90
89 class SpectralWidthPlot(RTIPlot):
91 class SpectralWidthPlot(RTIPlot):
@@ -99,7 +101,7 class SpectralWidthPlot(RTIPlot):
99 data = {
101 data = {
100 'width': dataOut.data_width
102 'width': dataOut.data_width
101 }
103 }
102
104 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
103 return data, {}
105 return data, {}
104
106
105 class SkyMapPlot(Plot):
107 class SkyMapPlot(Plot):
Show file before | Show file after | Hide comments
@@ -25,6 +25,7 class SpectraPlot(Plot):
25 channelList = []
25 channelList = []
26
26
27 def setup(self):
27 def setup(self):
28
28 self.nplots = len(self.data.channels)
29 self.nplots = len(self.data.channels)
29 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
30 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
30 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
31 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
@@ -37,16 +38,18 class SpectraPlot(Plot):
37 self.width = 3.5 * self.ncols
38 self.width = 3.5 * self.ncols
38 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
39 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
39 self.ylabel = 'Range [km]'
40 self.ylabel = 'Range [km]'
41
42
40 def update_list(self,dataOut):
43 def update_list(self,dataOut):
41 if len(self.channelList) == 0:
44 if len(self.channelList) == 0:
42 self.channelList = dataOut.channelList
45 self.channelList = dataOut.channelList
43
46
44 def update(self, dataOut):
47 def update(self, dataOut):
48
45 self.update_list(dataOut)
49 self.update_list(dataOut)
46 data = {}
50 data = {}
47 meta = {}
51 meta = {}
48 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
52 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
49
50 data['spc'] = spc
53 data['spc'] = spc
51 data['rti'] = dataOut.getPower()
54 data['rti'] = dataOut.getPower()
52 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
55 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
@@ -72,7 +75,6 class SpectraPlot(Plot):
72 self.xlabel = "Velocity (m/s)"
75 self.xlabel = "Velocity (m/s)"
73
76
74 self.titles = []
77 self.titles = []
75
76 y = self.data.yrange
78 y = self.data.yrange
77 self.y = y
79 self.y = y
78
80
@@ -192,8 +194,8 class CrossSpectraPlot(Plot):
192
194
193 if ax.firsttime:
195 if ax.firsttime:
194 ax.plt = ax.pcolormesh(x, y, coh.T,
196 ax.plt = ax.pcolormesh(x, y, coh.T,
195 vmin=0,
197 vmin=self.zmin_coh,
196 vmax=1,
198 vmax=self.zmax_coh,
197 cmap=plt.get_cmap(self.colormap_coh)
199 cmap=plt.get_cmap(self.colormap_coh)
198 )
200 )
199 else:
201 else:
@@ -210,6 +212,7 class CrossSpectraPlot(Plot):
210 )
212 )
211 else:
213 else:
212 ax.plt.set_array(phase.T.ravel())
214 ax.plt.set_array(phase.T.ravel())
215
213 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
216 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
214
217
215
218
@@ -258,6 +261,7 class RTIPlot(Plot):
258 self.z = self.data[self.CODE]
261 self.z = self.data[self.CODE]
259 self.z = numpy.array(self.z, dtype=float)
262 self.z = numpy.array(self.z, dtype=float)
260 self.z = numpy.ma.masked_invalid(self.z)
263 self.z = numpy.ma.masked_invalid(self.z)
264
261 try:
265 try:
262 if self.channelList != None:
266 if self.channelList != None:
263 self.titles = ['{} Channel {}'.format(
267 self.titles = ['{} Channel {}'.format(
@@ -282,9 +286,10 class RTIPlot(Plot):
282 cmap=plt.get_cmap(self.colormap)
286 cmap=plt.get_cmap(self.colormap)
283 )
287 )
284 if self.showprofile:
288 if self.showprofile:
285 ax.plot_profile = self.pf_axes[n].plot(
289 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
286 data['rti'][n], self.y)[0]
290
287 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
291 if "noise" in self.data:
292 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
288 color="k", linestyle="dashed", lw=1)[0]
293 color="k", linestyle="dashed", lw=1)[0]
289 else:
294 else:
290 ax.collections.remove(ax.collections[0])
295 ax.collections.remove(ax.collections[0])
@@ -294,8 +299,9 class RTIPlot(Plot):
294 cmap=plt.get_cmap(self.colormap)
299 cmap=plt.get_cmap(self.colormap)
295 )
300 )
296 if self.showprofile:
301 if self.showprofile:
297 ax.plot_profile.set_data(data['rti'][n], self.y)
302 ax.plot_profile.set_data(data[self.CODE][n], self.y)
298 ax.plot_noise.set_data(numpy.repeat(
303 if "noise" in self.data:
304 ax.plot_noise.set_data(numpy.repeat(
299 data['noise'][n], len(self.y)), self.y)
305 data['noise'][n], len(self.y)), self.y)
300
306
301
307
@@ -447,24 +453,40 class SpectraCutPlot(Plot):
447 CODE = 'spc_cut'
453 CODE = 'spc_cut'
448 plot_type = 'scatter'
454 plot_type = 'scatter'
449 buffering = False
455 buffering = False
456 heights = []
457 channelList = []
458 maintitle = "Spectra Cuts"
450
459
451 def setup(self):
460 def setup(self):
452
461
453 self.nplots = len(self.data.channels)
462 self.nplots = len(self.data.channels)
454 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
463 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
455 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
464 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
456 self.width = 3.4 * self.ncols + 1.5
465 self.width = 3.6 * self.ncols + 1.5
457 self.height = 3 * self.nrows
466 self.height = 3 * self.nrows
458 self.ylabel = 'Power [dB]'
467 self.ylabel = 'Power [dB]'
459 self.colorbar = False
468 self.colorbar = False
460 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
469 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
470 if self.selectedHeight:
471 self.maintitle = "Spectra Cut for %d km " %(int(self.selectedHeight))
461
472
462 def update(self, dataOut):
473 def update(self, dataOut):
474 if len(self.channelList) == 0:
475 self.channelList = dataOut.channelList
463
476
477 self.heights = dataOut.heightList
478 if self.selectedHeight:
479 index_list = numpy.where(self.heights >= self.selectedHeight)
480 self.height_index = index_list[0]
481 self.height_index = self.height_index[0]
482 #print(self.height_index)
464 data = {}
483 data = {}
465 meta = {}
484 meta = {}
466 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
485 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
467 data['spc'] = spc
486 if self.selectedHeight:
487 data['spc'] = spc[:,:,self.height_index]
488 else:
489 data['spc'] = spc
468 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
490 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
469
491
470 return data, meta
492 return data, meta
@@ -484,7 +506,7 class SpectraCutPlot(Plot):
484
506
485 y = self.data.yrange
507 y = self.data.yrange
486 z = self.data[-1]['spc']
508 z = self.data[-1]['spc']
487
509 #print(z.shape)
488 if self.height_index:
510 if self.height_index:
489 index = numpy.array(self.height_index)
511 index = numpy.array(self.height_index)
490 else:
512 else:
@@ -496,14 +518,21 class SpectraCutPlot(Plot):
496 self.xmin = self.xmin if self.xmin else -self.xmax
518 self.xmin = self.xmin if self.xmin else -self.xmax
497 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
519 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
498 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
520 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
499 ax.plt = ax.plot(x, z[n, :, index].T)
521 if self.selectedHeight:
500 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
522 ax.plt = ax.plot(x, z[n,:])
501 self.figures[0].legend(ax.plt, labels, loc='center right')
523
524 else:
525 ax.plt = ax.plot(x, z[n, :, index].T)
526 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
527 self.figures[0].legend(ax.plt, labels, loc='center right')
502 else:
528 else:
503 for i, line in enumerate(ax.plt):
529 for i, line in enumerate(ax.plt):
504 line.set_data(x, z[n, :, index[i]])
530 if self.selectedHeight:
505 self.titles.append('CH {}'.format(n))
531 line.set_data(x, z[n, :])
506
532 else:
533 line.set_data(x, z[n, :, index[i]])
534 self.titles.append('CH {}'.format(self.channelList[n]))
535 plt.suptitle(self.maintitle)
507
536
508 class BeaconPhase(Plot):
537 class BeaconPhase(Plot):
509
538
@@ -735,3 +764,198 class BeaconPhase(Plot):
735 update_figfile=update_figfile)
764 update_figfile=update_figfile)
736
765
737 return dataOut
766 return dataOut
767
768 class NoiselessSpectraPlot(Plot):
769 '''
770 Plot for Spectra data, subtracting
771 the noise in all channels, using for
772 amisr-14 data
773 '''
774
775 CODE = 'noiseless_spc'
776 colormap = 'nipy_spectral'
777 plot_type = 'pcolor'
778 buffering = False
779 channelList = []
780
781 def setup(self):
782
783 self.nplots = len(self.data.channels)
784 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
785 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
786 self.height = 2.6 * self.nrows
787
788 self.cb_label = 'dB'
789 if self.showprofile:
790 self.width = 4 * self.ncols
791 else:
792 self.width = 3.5 * self.ncols
793 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
794 self.ylabel = 'Range [km]'
795
796
797 def update_list(self,dataOut):
798 if len(self.channelList) == 0:
799 self.channelList = dataOut.channelList
800
801 def update(self, dataOut):
802
803 self.update_list(dataOut)
804 data = {}
805 meta = {}
806 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
807 (nch, nff, nh) = dataOut.data_spc.shape
808 n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
809 noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh))
810 #print(noise.shape, "noise", noise)
811
812 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise
813
814 data['spc'] = spc
815 data['rti'] = dataOut.getPower() - n1
816
817 data['noise'] = n0
818 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
819
820 return data, meta
821
822 def plot(self):
823 if self.xaxis == "frequency":
824 x = self.data.xrange[0]
825 self.xlabel = "Frequency (kHz)"
826 elif self.xaxis == "time":
827 x = self.data.xrange[1]
828 self.xlabel = "Time (ms)"
829 else:
830 x = self.data.xrange[2]
831 self.xlabel = "Velocity (m/s)"
832
833 self.titles = []
834 y = self.data.yrange
835 self.y = y
836
837 data = self.data[-1]
838 z = data['spc']
839
840 for n, ax in enumerate(self.axes):
841 noise = data['noise'][n]
842
843 if ax.firsttime:
844 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
845 self.xmin = self.xmin if self.xmin else -self.xmax
846 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
847 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
848 ax.plt = ax.pcolormesh(x, y, z[n].T,
849 vmin=self.zmin,
850 vmax=self.zmax,
851 cmap=plt.get_cmap(self.colormap)
852 )
853
854 if self.showprofile:
855 ax.plt_profile = self.pf_axes[n].plot(
856 data['rti'][n], y)[0]
857 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
858 color="k", linestyle="dashed", lw=1)[0]
859
860 else:
861 ax.plt.set_array(z[n].T.ravel())
862 if self.showprofile:
863 ax.plt_profile.set_data(data['rti'][n], y)
864 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
865
866 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
867
868
869 class NoiselessRTIPlot(Plot):
870 '''
871 Plot for RTI data
872 '''
873
874 CODE = 'noiseless_rti'
875 colormap = 'jet'
876 plot_type = 'pcolorbuffer'
877 titles = None
878 channelList = []
879
880 def setup(self):
881 self.xaxis = 'time'
882 self.ncols = 1
883 #print("dataChannels ",self.data.channels)
884 self.nrows = len(self.data.channels)
885 self.nplots = len(self.data.channels)
886 self.ylabel = 'Range [km]'
887 self.xlabel = 'Time'
888 self.cb_label = 'dB'
889 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
890 self.titles = ['{} Channel {}'.format(
891 self.CODE.upper(), x) for x in range(self.nplots)]
892
893 def update_list(self,dataOut):
894
895 self.channelList = dataOut.channelList
896
897
898 def update(self, dataOut):
899 if len(self.channelList) == 0:
900 self.update_list(dataOut)
901 data = {}
902 meta = {}
903
904 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
905 (nch, nff, nh) = dataOut.data_spc.shape
906 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
907
908
909 data['noiseless_rti'] = dataOut.getPower() - noise
910 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
911 return data, meta
912
913 def plot(self):
914
915 self.x = self.data.times
916 self.y = self.data.yrange
917 self.z = self.data['noiseless_rti']
918 self.z = numpy.array(self.z, dtype=float)
919 self.z = numpy.ma.masked_invalid(self.z)
920
921 try:
922 if self.channelList != None:
923 self.titles = ['{} Channel {}'.format(
924 self.CODE.upper(), x) for x in self.channelList]
925 except:
926 if self.channelList.any() != None:
927 self.titles = ['{} Channel {}'.format(
928 self.CODE.upper(), x) for x in self.channelList]
929 if self.decimation is None:
930 x, y, z = self.fill_gaps(self.x, self.y, self.z)
931 else:
932 x, y, z = self.fill_gaps(*self.decimate())
933 dummy_var = self.axes #ExtraΓ±amente esto actualiza el valor axes
934 for n, ax in enumerate(self.axes):
935 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
936 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
937 data = self.data[-1]
938 if ax.firsttime:
939 ax.plt = ax.pcolormesh(x, y, z[n].T,
940 vmin=self.zmin,
941 vmax=self.zmax,
942 cmap=plt.get_cmap(self.colormap)
943 )
944 if self.showprofile:
945 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
946
947 if "noise" in self.data:
948 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
949 color="k", linestyle="dashed", lw=1)[0]
950 else:
951 ax.collections.remove(ax.collections[0])
952 ax.plt = ax.pcolormesh(x, y, z[n].T,
953 vmin=self.zmin,
954 vmax=self.zmax,
955 cmap=plt.get_cmap(self.colormap)
956 )
957 if self.showprofile:
958 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
959 if "noise" in self.data:
960 ax.plot_noise.set_data(numpy.repeat(
961 data['noise'][n], len(self.y)), self.y)
Show file before | Show file after | Hide comments
@@ -681,7 +681,10 class Reader(object):
681
681
682 self.filename = filename
682 self.filename = filename
683 self.fileSize = os.path.getsize(filename)
683 self.fileSize = os.path.getsize(filename)
684 self.fp = self.open_file(filename, self.open_mode)
684 try:
685 self.fp = self.open_file(filename, self.open_mode)
686 except Exception as e:
687 raise schainpy.admin.SchainError("[Reading] Error in {} file, unable to open".format(filename))
685 self.flagIsNewFile = 1
688 self.flagIsNewFile = 1
686
689
687 return 1
690 return 1
@@ -691,7 +694,7 class Reader(object):
691 """Check if the given datetime is in range"""
694 """Check if the given datetime is in range"""
692 startDateTime= datetime.datetime.combine(startDate,startTime)
695 startDateTime= datetime.datetime.combine(startDate,startTime)
693 endDateTime = datetime.datetime.combine(endDate,endTime)
696 endDateTime = datetime.datetime.combine(endDate,endTime)
694
697
695 if startDateTime <= dt <= endDateTime:
698 if startDateTime <= dt <= endDateTime:
696 return True
699 return True
697 return False
700 return False
Show file before | Show file after | Hide comments
@@ -81,7 +81,7 class AMISRReader(ProcessingUnit):
81 #Is really necessary create the output object in the initializer
81 #Is really necessary create the output object in the initializer
82 self.dataOut = Voltage()
82 self.dataOut = Voltage()
83 self.dataOut.error=False
83 self.dataOut.error=False
84
84 self.margin_days = 1
85
85
86 def setup(self,path=None,
86 def setup(self,path=None,
87 startDate=None,
87 startDate=None,
@@ -95,7 +95,8 class AMISRReader(ProcessingUnit):
95 nCode = 0,
95 nCode = 0,
96 nBaud = 0,
96 nBaud = 0,
97 online=False,
97 online=False,
98 old_beams=False):
98 old_beams=False,
99 margin_days=1):
99
100
100
101
101
102
@@ -106,7 +107,7 class AMISRReader(ProcessingUnit):
106 self.code = code
107 self.code = code
107 self.nCode = int(nCode)
108 self.nCode = int(nCode)
108 self.nBaud = int(nBaud)
109 self.nBaud = int(nBaud)
109
110 self.margin_days = margin_days
110
111
111
112
112 #self.findFiles()
113 #self.findFiles()
@@ -195,7 +196,7 class AMISRReader(ProcessingUnit):
195 self.__nSamples = self.nsa
196 self.__nSamples = self.nsa
196 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
197 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
197 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
198 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
198
199 #print("amisr-ipp:",self.ippSeconds, self.__ippKm)
199 #for now until understand why the code saved is different (code included even though code not in tuf file)
200 #for now until understand why the code saved is different (code included even though code not in tuf file)
200 #self.__codeType = 0
201 #self.__codeType = 0
201 # self.__nCode = None
202 # self.__nCode = None
@@ -248,7 +249,7 class AMISRReader(ProcessingUnit):
248 dom = int(amisr_dirname_format[6:8])
249 dom = int(amisr_dirname_format[6:8])
249 thisDate = datetime.date(year,month,dom)
250 thisDate = datetime.date(year,month,dom)
250 #margen de un dΓ­a extra, igual luego se filtra for fecha y hora
251 #margen de un dΓ­a extra, igual luego se filtra for fecha y hora
251 if (thisDate>=(self.startDate - datetime.timedelta(days=1)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)):
252 if (thisDate>=(self.startDate - datetime.timedelta(days=self.margin_days)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)):
252 return amisr_dirname_format
253 return amisr_dirname_format
253 except:
254 except:
254 return None
255 return None
@@ -502,7 +503,6 class AMISRReader(ProcessingUnit):
502 nCode=self.__nCode, nBaud=self.__nBaud,
503 nCode=self.__nCode, nBaud=self.__nBaud,
503 code = self.__code,
504 code = self.__code,
504 fClock=1)
505 fClock=1)
505
506 #fill system header
506 #fill system header
507 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
507 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
508 nProfiles=self.newProfiles,
508 nProfiles=self.newProfiles,
Show file before | Show file after | Hide comments
@@ -415,7 +415,7 class HDFWrite(Operation):
415 dataAux = getattr(self.dataOut, self.dataList[i])
415 dataAux = getattr(self.dataOut, self.dataList[i])
416 dsDict['variable'] = self.dataList[i]
416 dsDict['variable'] = self.dataList[i]
417 else:
417 else:
418 log.warning('Attribute {} not found in dataOut', self.name)
418 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name)
419 continue
419 continue
420
420
421 if dataAux is None:
421 if dataAux is None:
Show file before | Show file after | Hide comments
@@ -57,15 +57,14 class ProcessingUnit(object):
57 '''
57 '''
58
58
59 try:
59 try:
60 # if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
60
61 # return self.dataIn.isReady()
61 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
62 #dataIn=None es unidades de Lectura, segunda parte unidades de procesamiento
62 return self.dataIn.isReady()
63 if self.dataIn is None or (not self.dataIn.error and not self.dataIn.flagNoData):
63 elif self.dataIn is None or not self.dataIn.error:
64 self.run(**kwargs)
64 self.run(**kwargs)
65 elif self.dataIn.error:
65 elif self.dataIn.error:
66 self.dataOut.error = self.dataIn.error
66 self.dataOut.error = self.dataIn.error
67 self.dataOut.flagNoData = True
67 self.dataOut.flagNoData = True
68
69 except:
68 except:
70
69
71 err = traceback.format_exc()
70 err = traceback.format_exc()
Show file before | Show file after | Hide comments
@@ -99,7 +99,7 class ParametersProc(ProcessingUnit):
99 self.dataOut.elevationList = self.dataIn.elevationList
99 self.dataOut.elevationList = self.dataIn.elevationList
100
100
101 def run(self):
101 def run(self):
102
102 # print("run parameter proc")
103 #---------------------- Voltage Data ---------------------------
103 #---------------------- Voltage Data ---------------------------
104
104
105 if self.dataIn.type == "Voltage":
105 if self.dataIn.type == "Voltage":
@@ -136,11 +136,13 class ParametersProc(ProcessingUnit):
136 self.dataOut.nIncohInt = self.dataIn.nIncohInt
136 self.dataOut.nIncohInt = self.dataIn.nIncohInt
137 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
137 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
138 self.dataOut.ippFactor = self.dataIn.ippFactor
138 self.dataOut.ippFactor = self.dataIn.ippFactor
139 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
140 self.dataOut.ipp = self.dataIn.ipp
139 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
141 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
140 self.dataOut.spc_noise = self.dataIn.getNoise()
142 self.dataOut.spc_noise = self.dataIn.getNoise()
141 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
143 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
142 # self.dataOut.normFactor = self.dataIn.normFactor
144 # self.dataOut.normFactor = self.dataIn.normFactor
143 if hasattr(self.dataIn, 'channelList'):
145 if hasattr(self.dataIn, 'channelList'):
144 self.dataOut.channelList = self.dataIn.channelList
146 self.dataOut.channelList = self.dataIn.channelList
145 if hasattr(self.dataIn, 'pairsList'):
147 if hasattr(self.dataIn, 'pairsList'):
146 self.dataOut.pairsList = self.dataIn.pairsList
148 self.dataOut.pairsList = self.dataIn.pairsList
@@ -308,7 +310,7 class SpectralFilters(Operation):
308 Operation.__init__(self)
310 Operation.__init__(self)
309 self.i = 0
311 self.i = 0
310
312
311 def run(self, dataOut, ):
313 def run(self, dataOut ):
312
314
313 self.spc = dataOut.data_pre[0].copy()
315 self.spc = dataOut.data_pre[0].copy()
314 self.Num_Chn = self.spc.shape[0]
316 self.Num_Chn = self.spc.shape[0]
@@ -1470,7 +1472,9 class SpectralFitting(Operation):
1470
1472
1471 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1473 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1472
1474
1473
1475 print("run SpectralFitting")
1476 self.dataOut = dataOut.copy()
1477 print(self.dataOut.nIncohInt)
1474 if path != None:
1478 if path != None:
1475 sys.path.append(path)
1479 sys.path.append(path)
1476 self.dataOut.library = importlib.import_module(file)
1480 self.dataOut.library = importlib.import_module(file)
@@ -1481,20 +1485,20 class SpectralFitting(Operation):
1481 self.dataOut.groupList = groupArray
1485 self.dataOut.groupList = groupArray
1482
1486
1483 nGroups = groupArray.shape[0]
1487 nGroups = groupArray.shape[0]
1484 nChannels = self.dataIn.nChannels
1488 nChannels = self.dataOut.nChannels
1485 nHeights=self.dataIn.heightList.size
1489 nHeights=self.dataOut.heightList.size
1486
1490
1487 #Parameters Array
1491 #Parameters Array
1488 self.dataOut.data_param = None
1492 self.dataOut.data_param = None
1489
1493
1490 #Set constants
1494 #Set constants
1491 constants = self.dataOut.library.setConstants(self.dataIn)
1495 constants = self.dataOut.library.setConstants(self.dataOut)
1492 self.dataOut.constants = constants
1496 self.dataOut.constants = constants
1493 M = self.dataIn.normFactor
1497 M = self.dataOut.normFactor
1494 N = self.dataIn.nFFTPoints
1498 N = self.dataOut.nFFTPoints
1495 ippSeconds = self.dataIn.ippSeconds
1499 ippSeconds = self.dataOut.ippSeconds
1496 K = self.dataIn.nIncohInt
1500 K = self.dataOut.nIncohInt
1497 pairsArray = numpy.array(self.dataIn.pairsList)
1501 pairsArray = numpy.array(self.dataOut.pairsList)
1498
1502
1499 #List of possible combinations
1503 #List of possible combinations
1500 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1504 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
@@ -1503,14 +1507,14 class SpectralFitting(Operation):
1503 if getSNR:
1507 if getSNR:
1504 listChannels = groupArray.reshape((groupArray.size))
1508 listChannels = groupArray.reshape((groupArray.size))
1505 listChannels.sort()
1509 listChannels.sort()
1506 noise = self.dataIn.getNoise()
1510 noise = self.dataOut.getNoise()
1507 self.dataOut.data_snr = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1511 self.dataOut.data_snr = self.__getSNR(self.dataOut.data_spc[listChannels,:,:], noise[listChannels])
1508
1512
1509 for i in range(nGroups):
1513 for i in range(nGroups):
1510 coord = groupArray[i,:]
1514 coord = groupArray[i,:]
1511
1515
1512 #Input data array
1516 #Input data array
1513 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1517 data = self.dataOut.data_spc[coord,:,:]/(M*N)
1514 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1518 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1515
1519
1516 #Cross Spectra data array for Covariance Matrixes
1520 #Cross Spectra data array for Covariance Matrixes
@@ -1519,7 +1523,7 class SpectralFitting(Operation):
1519 pairsSel = numpy.array([coord[x],coord[y]])
1523 pairsSel = numpy.array([coord[x],coord[y]])
1520 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1524 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1521 ind += 1
1525 ind += 1
1522 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1526 dataCross = self.dataOut.data_cspc[indCross,:,:]/(M*N)
1523 dataCross = dataCross**2/K
1527 dataCross = dataCross**2/K
1524
1528
1525 for h in range(nHeights):
1529 for h in range(nHeights):
@@ -1548,13 +1552,13 class SpectralFitting(Operation):
1548 dp = numpy.dot(LT,d)
1552 dp = numpy.dot(LT,d)
1549
1553
1550 #Initial values
1554 #Initial values
1551 data_spc = self.dataIn.data_spc[coord,:,h]
1555 data_spc = self.dataOut.data_spc[coord,:,h]
1552
1556
1553 if (h>0)and(error1[3]<5):
1557 if (h>0)and(error1[3]<5):
1554 p0 = self.dataOut.data_param[i,:,h-1]
1558 p0 = self.dataOut.data_param[i,:,h-1]
1555 else:
1559 else:
1556 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1560 #p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1557
1561 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants))
1558 try:
1562 try:
1559 #Least Squares
1563 #Least Squares
1560 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1564 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
@@ -1563,7 +1567,8 class SpectralFitting(Operation):
1563 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1567 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1564 #Error with Jacobian
1568 #Error with Jacobian
1565 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1569 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1566 except:
1570 except Exception as e:
1571 print("Param error: ",e)
1567 minp = p0*numpy.nan
1572 minp = p0*numpy.nan
1568 error0 = numpy.nan
1573 error0 = numpy.nan
1569 error1 = p0*numpy.nan
1574 error1 = p0*numpy.nan
@@ -1575,7 +1580,11 class SpectralFitting(Operation):
1575
1580
1576 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1581 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1577 self.dataOut.data_param[i,:,h] = minp
1582 self.dataOut.data_param[i,:,h] = minp
1578 return
1583
1584
1585 print("SpectralFitting Done")
1586 dataOut = self.dataOut
1587 return dataOut
1579
1588
1580 def __residFunction(self, p, dp, LT, constants):
1589 def __residFunction(self, p, dp, LT, constants):
1581
1590
@@ -2320,8 +2329,8 class WindProfiler(Operation):
2320
2329
2321 class EWDriftsEstimation(Operation):
2330 class EWDriftsEstimation(Operation):
2322
2331
2323 def __init__(self):
2332 # def __init__(self):
2324 Operation.__init__(self)
2333 # Operation.__init__(self)
2325
2334
2326 def __correctValues(self, heiRang, phi, velRadial, SNR):
2335 def __correctValues(self, heiRang, phi, velRadial, SNR):
2327 listPhi = phi.tolist()
2336 listPhi = phi.tolist()
@@ -2338,25 +2347,27 class EWDriftsEstimation(Operation):
2338
2347
2339 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2348 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2340 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2349 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2341
2350 try:
2342 for i in rango:
2351 for i in rango:
2343 x = heiRang*math.cos(phi[i])
2352 x = heiRang*math.cos(phi[i])
2344 y1 = velRadial[i,:]
2353 y1 = velRadial[i,:]
2345 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2354 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2346
2355
2347 x1 = heiRang1
2356 x1 = heiRang1
2348 y11 = f1(x1)
2357 y11 = f1(x1)
2349
2358
2350 y2 = SNR[i,:]
2359 y2 = SNR[i,:]
2351 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2360 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2352 y21 = f2(x1)
2361 y21 = f2(x1)
2353
2362
2354 velRadial1[i,:] = y11
2363 velRadial1[i,:] = y11
2355 SNR1[i,:] = y21
2364 SNR1[i,:] = y21
2356
2365 except Exception as e:
2366 print("wrong values!",e)
2357 return heiRang1, velRadial1, SNR1
2367 return heiRang1, velRadial1, SNR1
2358
2368
2359 def run(self, dataOut, zenith, zenithCorrection):
2369 def run(self, dataOut, zenith, zenithCorrection):
2370 print("run EWDriftsEstimation")
2360 heiRang = dataOut.heightList
2371 heiRang = dataOut.heightList
2361 velRadial = dataOut.data_param[:,3,:]
2372 velRadial = dataOut.data_param[:,3,:]
2362 SNR = dataOut.data_snr
2373 SNR = dataOut.data_snr
@@ -2366,7 +2377,7 class EWDriftsEstimation(Operation):
2366 zenith *= numpy.pi/180
2377 zenith *= numpy.pi/180
2367
2378
2368 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2379 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2369
2380 print("__correctValues Done ")
2370 alp = zenith[0]
2381 alp = zenith[0]
2371 bet = zenith[1]
2382 bet = zenith[1]
2372
2383
@@ -2384,6 +2395,8 class EWDriftsEstimation(Operation):
2384
2395
2385 dataOut.utctimeInit = dataOut.utctime
2396 dataOut.utctimeInit = dataOut.utctime
2386 dataOut.outputInterval = dataOut.timeInterval
2397 dataOut.outputInterval = dataOut.timeInterval
2398
2399 print("EWDriftsEstimation Done ")
2387 return
2400 return
2388
2401
2389 #--------------- Non Specular Meteor ----------------
2402 #--------------- Non Specular Meteor ----------------
Show file before | Show file after | Hide comments
@@ -68,7 +68,6 class SpectraProc(ProcessingUnit):
68 self.dataOut.elevationList = self.dataIn.elevationList
68 self.dataOut.elevationList = self.dataIn.elevationList
69
69
70
70
71
72 def __getFft(self):
71 def __getFft(self):
73 """
72 """
74 Convierte valores de Voltaje a Spectra
73 Convierte valores de Voltaje a Spectra
Show file before | Show file after | Hide comments
@@ -178,7 +178,7 class selectHeights(Operation):
178 1 si el metodo se ejecuto con exito caso contrario devuelve 0
178 1 si el metodo se ejecuto con exito caso contrario devuelve 0
179 """
179 """
180
180
181 dataOut = dataOut
181 self.dataOut = dataOut
182
182
183 if minHei and maxHei:
183 if minHei and maxHei:
184
184
@@ -207,7 +207,7 class selectHeights(Operation):
207
207
208 self.selectHeightsByIndex(minIndex, maxIndex)
208 self.selectHeightsByIndex(minIndex, maxIndex)
209
209
210 return self.dataOut
210 return dataOut
211
211
212 def selectHeightsByIndex(self, minIndex, maxIndex):
212 def selectHeightsByIndex(self, minIndex, maxIndex):
213 """
213 """
@@ -1636,3 +1636,123 class PulsePairVoltage(Operation):
1636 # self.__startIndex += self.__newNSamples
1636 # self.__startIndex += self.__newNSamples
1637 #
1637 #
1638 # return
1638 # return
1639 class SSheightProfiles(Operation):
1640
1641 step = None
1642 nsamples = None
1643 bufferShape = None
1644 profileShape = None
1645 sshProfiles = None
1646 profileIndex = None
1647
1648 def __init__(self, **kwargs):
1649
1650 Operation.__init__(self, **kwargs)
1651 self.isConfig = False
1652
1653 def setup(self,dataOut ,step = None , nsamples = None):
1654
1655 if step == None and nsamples == None:
1656 raise ValueError("step or nheights should be specified ...")
1657
1658 self.step = step
1659 self.nsamples = nsamples
1660 self.__nChannels = dataOut.nChannels
1661 self.__nProfiles = dataOut.nProfiles
1662 self.__nHeis = dataOut.nHeights
1663 shape = dataOut.data.shape #nchannels, nprofiles, nsamples
1664
1665 residue = (shape[1] - self.nsamples) % self.step
1666 if residue != 0:
1667 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
1668
1669 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1670 numberProfile = self.nsamples
1671 numberSamples = (shape[1] - self.nsamples)/self.step
1672
1673 self.bufferShape = int(shape[0]), int(numberSamples), int(numberProfile) # nchannels, nsamples , nprofiles
1674 self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples) # nchannels, nprofiles, nsamples
1675
1676 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1677 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1678
1679 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1680 dataOut.flagNoData = True
1681
1682 profileIndex = None
1683 #print(dataOut.getFreqRange(1)/1000.)
1684 #exit(1)
1685 if dataOut.flagDataAsBlock:
1686 dataOut.data = numpy.average(dataOut.data,axis=1)
1687 #print("jee")
1688 dataOut.flagDataAsBlock = False
1689 if not self.isConfig:
1690 self.setup(dataOut, step=step , nsamples=nsamples)
1691 #print("Setup done")
1692 self.isConfig = True
1693
1694 #DC_Hae = numpy.array([0.398+0.588j, -0.926+0.306j, -0.536-0.682j, -0.072+0.53j, 0.368-0.356j, 0.996+0.362j])
1695 #DC_Hae = numpy.array([ 0.001025 +0.0516375j, 0.03485 +0.20923125j, -0.168 -0.02720625j,
1696 #-0.1105375 +0.0707125j, -0.20309375-0.09670625j, 0.189775 +0.02716875j])*(-3.5)
1697
1698 #DC_Hae = numpy.array([ -32.26 +8.66j, -32.26 +8.66j])
1699
1700 #DC_Hae = numpy.array([-2.78500000e-01 -1.39175j, -6.63237294e+02+210.4268625j])
1701
1702 #print(dataOut.data[0,13:15])
1703 #dataOut.data = dataOut.data - DC_Hae[:,None]
1704 #print(dataOut.data[0,13:15])
1705 #exit(1)
1706 if code is not None:
1707 code = numpy.array(code)
1708 code_block = code
1709 '''
1710 roll = 0
1711 code = numpy.roll(code,roll,axis=0)
1712 code = numpy.reshape(code,(5,100,64))
1713 block = dataOut.CurrentBlock%5
1714
1715 day_dif = 0 #day_19_Oct_2021: 3
1716 code_block = code[block-1-day_dif,:,:]
1717 '''
1718 if repeat is not None:
1719 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1720 #print(code_block.shape)
1721 for i in range(self.buffer.shape[1]):
1722
1723 if code is not None:
1724 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block
1725
1726 else:
1727
1728 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1729
1730 #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
1731
1732 for j in range(self.buffer.shape[0]):
1733 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1734
1735 profileIndex = self.nsamples
1736 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1737 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1738 #print("ippSeconds, dH: ",ippSeconds,deltaHeight)
1739 try:
1740 if dataOut.concat_m is not None:
1741 ippSeconds= ippSeconds/float(dataOut.concat_m)
1742 #print "Profile concat %d"%dataOut.concat_m
1743 except:
1744 pass
1745
1746 dataOut.data = self.sshProfiles
1747 dataOut.flagNoData = False
1748 dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
1749 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
1750
1751 dataOut.profileIndex = profileIndex
1752 dataOut.flagDataAsBlock = True
1753 dataOut.ippSeconds = ippSeconds
1754 dataOut.step = self.step
1755 #print(numpy.shape(dataOut.data))
1756 #exit(1)
1757
1758 return dataOut
General Comments 0
You need to be logged in to leave comments. Login now