##// 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 257 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
258 258 self.tmin = kwargs.get('tmin', None)
259 259 self.t_units = kwargs.get('t_units', "h_m")
260 self.selectedHeight = kwargs.get('selectedHeight', None)
261
260 262
261 263 if self.server:
262 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 14 class SpectraHeisPlot(Plot):
15 15
16 16 CODE = 'spc_heis'
17
17 channelList = []
18 18 def setup(self):
19 19
20 20 self.nplots = len(self.data.channels)
@@ -28,6 +28,8 class SpectraHeisPlot(Plot):
28 28 self.colorbar = False
29 29
30 30 def update(self, dataOut):
31 if len(self.channelList) == 0:
32 self.channelList = dataOut.channelList
31 33
32 34 data = {}
33 35 meta = {}
@@ -51,9 +53,12 class SpectraHeisPlot(Plot):
51 53 self.xmin = min(x) if self.xmin is None else self.xmin
52 54 self.xmax = max(x) if self.xmax is None else self.xmax
53 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 58 else:
55 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 62 self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel)))
58 63 plt.suptitle(Maintitle)
59 64
Show file before | Show file after | Hide comments
@@ -48,12 +48,13 class SnrPlot(RTIPlot):
48 48 colormap = 'jet'
49 49
50 50 def update(self, dataOut):
51 self.update_list(dataOut)
52 data = {
53 'snr': 10*numpy.log10(dataOut.data_snr)
54 }
55
56 return data, {}
51 if len(self.channelList) == 0:
52 self.update_list(dataOut)
53 data = {}
54 meta = {}
55 data['snr'] = 10*numpy.log10(dataOut.data_snr)
56
57 return data, meta
57 58
58 59 class DopplerPlot(RTIPlot):
59 60 '''
@@ -84,6 +85,7 class PowerPlot(RTIPlot):
84 85 data = {
85 86 'pow': 10*numpy.log10(dataOut.data_pow)
86 87 }
88 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
87 89 return data, {}
88 90
89 91 class SpectralWidthPlot(RTIPlot):
@@ -99,7 +101,7 class SpectralWidthPlot(RTIPlot):
99 101 data = {
100 102 'width': dataOut.data_width
101 103 }
102
104 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
103 105 return data, {}
104 106
105 107 class SkyMapPlot(Plot):
Show file before | Show file after | Hide comments
@@ -25,6 +25,7 class SpectraPlot(Plot):
25 25 channelList = []
26 26
27 27 def setup(self):
28
28 29 self.nplots = len(self.data.channels)
29 30 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
30 31 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
@@ -37,16 +38,18 class SpectraPlot(Plot):
37 38 self.width = 3.5 * self.ncols
38 39 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
39 40 self.ylabel = 'Range [km]'
41
42
40 43 def update_list(self,dataOut):
41 44 if len(self.channelList) == 0:
42 45 self.channelList = dataOut.channelList
43 46
44 47 def update(self, dataOut):
48
45 49 self.update_list(dataOut)
46 50 data = {}
47 51 meta = {}
48 52 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
49
50 53 data['spc'] = spc
51 54 data['rti'] = dataOut.getPower()
52 55 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
@@ -72,7 +75,6 class SpectraPlot(Plot):
72 75 self.xlabel = "Velocity (m/s)"
73 76
74 77 self.titles = []
75
76 78 y = self.data.yrange
77 79 self.y = y
78 80
@@ -192,8 +194,8 class CrossSpectraPlot(Plot):
192 194
193 195 if ax.firsttime:
194 196 ax.plt = ax.pcolormesh(x, y, coh.T,
195 vmin=0,
196 vmax=1,
197 vmin=self.zmin_coh,
198 vmax=self.zmax_coh,
197 199 cmap=plt.get_cmap(self.colormap_coh)
198 200 )
199 201 else:
@@ -210,6 +212,7 class CrossSpectraPlot(Plot):
210 212 )
211 213 else:
212 214 ax.plt.set_array(phase.T.ravel())
215
213 216 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
214 217
215 218
@@ -258,6 +261,7 class RTIPlot(Plot):
258 261 self.z = self.data[self.CODE]
259 262 self.z = numpy.array(self.z, dtype=float)
260 263 self.z = numpy.ma.masked_invalid(self.z)
264
261 265 try:
262 266 if self.channelList != None:
263 267 self.titles = ['{} Channel {}'.format(
@@ -282,9 +286,10 class RTIPlot(Plot):
282 286 cmap=plt.get_cmap(self.colormap)
283 287 )
284 288 if self.showprofile:
285 ax.plot_profile = self.pf_axes[n].plot(
286 data['rti'][n], self.y)[0]
287 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
289 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
290
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 293 color="k", linestyle="dashed", lw=1)[0]
289 294 else:
290 295 ax.collections.remove(ax.collections[0])
@@ -294,8 +299,9 class RTIPlot(Plot):
294 299 cmap=plt.get_cmap(self.colormap)
295 300 )
296 301 if self.showprofile:
297 ax.plot_profile.set_data(data['rti'][n], self.y)
298 ax.plot_noise.set_data(numpy.repeat(
302 ax.plot_profile.set_data(data[self.CODE][n], self.y)
303 if "noise" in self.data:
304 ax.plot_noise.set_data(numpy.repeat(
299 305 data['noise'][n], len(self.y)), self.y)
300 306
301 307
@@ -447,24 +453,40 class SpectraCutPlot(Plot):
447 453 CODE = 'spc_cut'
448 454 plot_type = 'scatter'
449 455 buffering = False
456 heights = []
457 channelList = []
458 maintitle = "Spectra Cuts"
450 459
451 460 def setup(self):
452 461
453 462 self.nplots = len(self.data.channels)
454 463 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
455 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 466 self.height = 3 * self.nrows
458 467 self.ylabel = 'Power [dB]'
459 468 self.colorbar = False
460 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 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 483 data = {}
465 484 meta = {}
466 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 490 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
469 491
470 492 return data, meta
@@ -484,7 +506,7 class SpectraCutPlot(Plot):
484 506
485 507 y = self.data.yrange
486 508 z = self.data[-1]['spc']
487
509 #print(z.shape)
488 510 if self.height_index:
489 511 index = numpy.array(self.height_index)
490 512 else:
@@ -496,14 +518,21 class SpectraCutPlot(Plot):
496 518 self.xmin = self.xmin if self.xmin else -self.xmax
497 519 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
498 520 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
499 ax.plt = ax.plot(x, z[n, :, index].T)
500 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
501 self.figures[0].legend(ax.plt, labels, loc='center right')
521 if self.selectedHeight:
522 ax.plt = ax.plot(x, z[n,:])
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 528 else:
503 529 for i, line in enumerate(ax.plt):
504 line.set_data(x, z[n, :, index[i]])
505 self.titles.append('CH {}'.format(n))
506
530 if self.selectedHeight:
531 line.set_data(x, z[n, :])
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 537 class BeaconPhase(Plot):
509 538
@@ -735,3 +764,198 class BeaconPhase(Plot):
735 764 update_figfile=update_figfile)
736 765
737 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 682 self.filename = filename
683 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 688 self.flagIsNewFile = 1
686 689
687 690 return 1
@@ -691,7 +694,7 class Reader(object):
691 694 """Check if the given datetime is in range"""
692 695 startDateTime= datetime.datetime.combine(startDate,startTime)
693 696 endDateTime = datetime.datetime.combine(endDate,endTime)
694
697
695 698 if startDateTime <= dt <= endDateTime:
696 699 return True
697 700 return False
Show file before | Show file after | Hide comments
@@ -81,7 +81,7 class AMISRReader(ProcessingUnit):
81 81 #Is really necessary create the output object in the initializer
82 82 self.dataOut = Voltage()
83 83 self.dataOut.error=False
84
84 self.margin_days = 1
85 85
86 86 def setup(self,path=None,
87 87 startDate=None,
@@ -95,7 +95,8 class AMISRReader(ProcessingUnit):
95 95 nCode = 0,
96 96 nBaud = 0,
97 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 107 self.code = code
107 108 self.nCode = int(nCode)
108 109 self.nBaud = int(nBaud)
109
110 self.margin_days = margin_days
110 111
111 112
112 113 #self.findFiles()
@@ -195,7 +196,7 class AMISRReader(ProcessingUnit):
195 196 self.__nSamples = self.nsa
196 197 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
197 198 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
198
199 #print("amisr-ipp:",self.ippSeconds, self.__ippKm)
199 200 #for now until understand why the code saved is different (code included even though code not in tuf file)
200 201 #self.__codeType = 0
201 202 # self.__nCode = None
@@ -248,7 +249,7 class AMISRReader(ProcessingUnit):
248 249 dom = int(amisr_dirname_format[6:8])
249 250 thisDate = datetime.date(year,month,dom)
250 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 253 return amisr_dirname_format
253 254 except:
254 255 return None
@@ -502,7 +503,6 class AMISRReader(ProcessingUnit):
502 503 nCode=self.__nCode, nBaud=self.__nBaud,
503 504 code = self.__code,
504 505 fClock=1)
505
506 506 #fill system header
507 507 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
508 508 nProfiles=self.newProfiles,
Show file before | Show file after | Hide comments
@@ -415,7 +415,7 class HDFWrite(Operation):
415 415 dataAux = getattr(self.dataOut, self.dataList[i])
416 416 dsDict['variable'] = self.dataList[i]
417 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 419 continue
420 420
421 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 59 try:
60 # if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
61 # return self.dataIn.isReady()
62 #dataIn=None es unidades de Lectura, segunda parte unidades de procesamiento
63 if self.dataIn is None or (not self.dataIn.error and not self.dataIn.flagNoData):
60
61 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
62 return self.dataIn.isReady()
63 elif self.dataIn is None or not self.dataIn.error:
64 64 self.run(**kwargs)
65 65 elif self.dataIn.error:
66 66 self.dataOut.error = self.dataIn.error
67 67 self.dataOut.flagNoData = True
68
69 68 except:
70 69
71 70 err = traceback.format_exc()
Show file before | Show file after | Hide comments
@@ -99,7 +99,7 class ParametersProc(ProcessingUnit):
99 99 self.dataOut.elevationList = self.dataIn.elevationList
100 100
101 101 def run(self):
102
102 # print("run parameter proc")
103 103 #---------------------- Voltage Data ---------------------------
104 104
105 105 if self.dataIn.type == "Voltage":
@@ -136,11 +136,13 class ParametersProc(ProcessingUnit):
136 136 self.dataOut.nIncohInt = self.dataIn.nIncohInt
137 137 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
138 138 self.dataOut.ippFactor = self.dataIn.ippFactor
139 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
140 self.dataOut.ipp = self.dataIn.ipp
139 141 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
140 142 self.dataOut.spc_noise = self.dataIn.getNoise()
141 143 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
142 144 # self.dataOut.normFactor = self.dataIn.normFactor
143 if hasattr(self.dataIn, 'channelList'):
145 if hasattr(self.dataIn, 'channelList'):
144 146 self.dataOut.channelList = self.dataIn.channelList
145 147 if hasattr(self.dataIn, 'pairsList'):
146 148 self.dataOut.pairsList = self.dataIn.pairsList
@@ -308,7 +310,7 class SpectralFilters(Operation):
308 310 Operation.__init__(self)
309 311 self.i = 0
310 312
311 def run(self, dataOut, ):
313 def run(self, dataOut ):
312 314
313 315 self.spc = dataOut.data_pre[0].copy()
314 316 self.Num_Chn = self.spc.shape[0]
@@ -1470,7 +1472,9 class SpectralFitting(Operation):
1470 1472
1471 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 1478 if path != None:
1475 1479 sys.path.append(path)
1476 1480 self.dataOut.library = importlib.import_module(file)
@@ -1481,20 +1485,20 class SpectralFitting(Operation):
1481 1485 self.dataOut.groupList = groupArray
1482 1486
1483 1487 nGroups = groupArray.shape[0]
1484 nChannels = self.dataIn.nChannels
1485 nHeights=self.dataIn.heightList.size
1488 nChannels = self.dataOut.nChannels
1489 nHeights=self.dataOut.heightList.size
1486 1490
1487 1491 #Parameters Array
1488 1492 self.dataOut.data_param = None
1489 1493
1490 1494 #Set constants
1491 constants = self.dataOut.library.setConstants(self.dataIn)
1495 constants = self.dataOut.library.setConstants(self.dataOut)
1492 1496 self.dataOut.constants = constants
1493 M = self.dataIn.normFactor
1494 N = self.dataIn.nFFTPoints
1495 ippSeconds = self.dataIn.ippSeconds
1496 K = self.dataIn.nIncohInt
1497 pairsArray = numpy.array(self.dataIn.pairsList)
1497 M = self.dataOut.normFactor
1498 N = self.dataOut.nFFTPoints
1499 ippSeconds = self.dataOut.ippSeconds
1500 K = self.dataOut.nIncohInt
1501 pairsArray = numpy.array(self.dataOut.pairsList)
1498 1502
1499 1503 #List of possible combinations
1500 1504 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
@@ -1503,14 +1507,14 class SpectralFitting(Operation):
1503 1507 if getSNR:
1504 1508 listChannels = groupArray.reshape((groupArray.size))
1505 1509 listChannels.sort()
1506 noise = self.dataIn.getNoise()
1507 self.dataOut.data_snr = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1510 noise = self.dataOut.getNoise()
1511 self.dataOut.data_snr = self.__getSNR(self.dataOut.data_spc[listChannels,:,:], noise[listChannels])
1508 1512
1509 1513 for i in range(nGroups):
1510 1514 coord = groupArray[i,:]
1511 1515
1512 1516 #Input data array
1513 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1517 data = self.dataOut.data_spc[coord,:,:]/(M*N)
1514 1518 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1515 1519
1516 1520 #Cross Spectra data array for Covariance Matrixes
@@ -1519,7 +1523,7 class SpectralFitting(Operation):
1519 1523 pairsSel = numpy.array([coord[x],coord[y]])
1520 1524 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1521 1525 ind += 1
1522 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1526 dataCross = self.dataOut.data_cspc[indCross,:,:]/(M*N)
1523 1527 dataCross = dataCross**2/K
1524 1528
1525 1529 for h in range(nHeights):
@@ -1548,13 +1552,13 class SpectralFitting(Operation):
1548 1552 dp = numpy.dot(LT,d)
1549 1553
1550 1554 #Initial values
1551 data_spc = self.dataIn.data_spc[coord,:,h]
1555 data_spc = self.dataOut.data_spc[coord,:,h]
1552 1556
1553 1557 if (h>0)and(error1[3]<5):
1554 1558 p0 = self.dataOut.data_param[i,:,h-1]
1555 1559 else:
1556 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1557
1560 #p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1561 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants))
1558 1562 try:
1559 1563 #Least Squares
1560 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 1567 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1564 1568 #Error with Jacobian
1565 1569 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1566 except:
1570 except Exception as e:
1571 print("Param error: ",e)
1567 1572 minp = p0*numpy.nan
1568 1573 error0 = numpy.nan
1569 1574 error1 = p0*numpy.nan
@@ -1575,7 +1580,11 class SpectralFitting(Operation):
1575 1580
1576 1581 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1577 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 1589 def __residFunction(self, p, dp, LT, constants):
1581 1590
@@ -2320,8 +2329,8 class WindProfiler(Operation):
2320 2329
2321 2330 class EWDriftsEstimation(Operation):
2322 2331
2323 def __init__(self):
2324 Operation.__init__(self)
2332 # def __init__(self):
2333 # Operation.__init__(self)
2325 2334
2326 2335 def __correctValues(self, heiRang, phi, velRadial, SNR):
2327 2336 listPhi = phi.tolist()
@@ -2338,25 +2347,27 class EWDriftsEstimation(Operation):
2338 2347
2339 2348 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2340 2349 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2341
2342 for i in rango:
2343 x = heiRang*math.cos(phi[i])
2344 y1 = velRadial[i,:]
2345 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2346
2347 x1 = heiRang1
2348 y11 = f1(x1)
2349
2350 y2 = SNR[i,:]
2351 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2352 y21 = f2(x1)
2353
2354 velRadial1[i,:] = y11
2355 SNR1[i,:] = y21
2356
2350 try:
2351 for i in rango:
2352 x = heiRang*math.cos(phi[i])
2353 y1 = velRadial[i,:]
2354 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2355
2356 x1 = heiRang1
2357 y11 = f1(x1)
2358
2359 y2 = SNR[i,:]
2360 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2361 y21 = f2(x1)
2362
2363 velRadial1[i,:] = y11
2364 SNR1[i,:] = y21
2365 except Exception as e:
2366 print("wrong values!",e)
2357 2367 return heiRang1, velRadial1, SNR1
2358 2368
2359 2369 def run(self, dataOut, zenith, zenithCorrection):
2370 print("run EWDriftsEstimation")
2360 2371 heiRang = dataOut.heightList
2361 2372 velRadial = dataOut.data_param[:,3,:]
2362 2373 SNR = dataOut.data_snr
@@ -2366,7 +2377,7 class EWDriftsEstimation(Operation):
2366 2377 zenith *= numpy.pi/180
2367 2378
2368 2379 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2369
2380 print("__correctValues Done ")
2370 2381 alp = zenith[0]
2371 2382 bet = zenith[1]
2372 2383
@@ -2384,6 +2395,8 class EWDriftsEstimation(Operation):
2384 2395
2385 2396 dataOut.utctimeInit = dataOut.utctime
2386 2397 dataOut.outputInterval = dataOut.timeInterval
2398
2399 print("EWDriftsEstimation Done ")
2387 2400 return
2388 2401
2389 2402 #--------------- Non Specular Meteor ----------------
Show file before | Show file after | Hide comments
@@ -68,7 +68,6 class SpectraProc(ProcessingUnit):
68 68 self.dataOut.elevationList = self.dataIn.elevationList
69 69
70 70
71
72 71 def __getFft(self):
73 72 """
74 73 Convierte valores de Voltaje a Spectra
Show file before | Show file after | Hide comments
@@ -178,7 +178,7 class selectHeights(Operation):
178 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 183 if minHei and maxHei:
184 184
@@ -207,7 +207,7 class selectHeights(Operation):
207 207
208 208 self.selectHeightsByIndex(minIndex, maxIndex)
209 209
210 return self.dataOut
210 return dataOut
211 211
212 212 def selectHeightsByIndex(self, minIndex, maxIndex):
213 213 """
@@ -1636,3 +1636,123 class PulsePairVoltage(Operation):
1636 1636 # self.__startIndex += self.__newNSamples
1637 1637 #
1638 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