diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index 517062e..db0fd85 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -257,6 +257,8 @@ class Plot(Operation): self.data = PlotterData(self.CODE, self.exp_code, self.localtime) self.tmin = kwargs.get('tmin', None) self.t_units = kwargs.get('t_units', "h_m") + self.selectedHeight = kwargs.get('selectedHeight', None) + if self.server: if not self.server.startswith('tcp://'): diff --git a/schainpy/model/graphics/jroplot_heispectra.py b/schainpy/model/graphics/jroplot_heispectra.py index da391ea..8486c26 100644 --- a/schainpy/model/graphics/jroplot_heispectra.py +++ b/schainpy/model/graphics/jroplot_heispectra.py @@ -14,7 +14,7 @@ import matplotlib.pyplot as plt class SpectraHeisPlot(Plot): CODE = 'spc_heis' - + channelList = [] def setup(self): self.nplots = len(self.data.channels) @@ -28,6 +28,8 @@ class SpectraHeisPlot(Plot): self.colorbar = False def update(self, dataOut): + if len(self.channelList) == 0: + self.channelList = dataOut.channelList data = {} meta = {} @@ -51,9 +53,12 @@ class SpectraHeisPlot(Plot): self.xmin = min(x) if self.xmin is None else self.xmin self.xmax = max(x) if self.xmax is None else self.xmax ax.plt = ax.plot(x, ychannel, lw=1, color='b')[0] + ax.set_ylim(ymin=self.zmin, ymax=self.zmax) + ax.set_xlim(xmin=self.xmin, xmax=self.xmax) else: ax.plt.set_data(x, ychannel) - + ax.set_ylim(ymin=self.zmin, ymax=self.zmax) + ax.set_xlim(xmin=self.xmin, xmax=self.xmax) self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel))) plt.suptitle(Maintitle) diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index 732060d..b4c4216 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -48,12 +48,13 @@ class SnrPlot(RTIPlot): colormap = 'jet' def update(self, dataOut): - self.update_list(dataOut) - data = { - 'snr': 10*numpy.log10(dataOut.data_snr) - } - - return data, {} + if len(self.channelList) == 0: + self.update_list(dataOut) + data = {} + meta = {} + data['snr'] = 10*numpy.log10(dataOut.data_snr) + + return data, meta class DopplerPlot(RTIPlot): ''' @@ -84,6 +85,7 @@ class PowerPlot(RTIPlot): data = { 'pow': 10*numpy.log10(dataOut.data_pow) } + data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) return data, {} class SpectralWidthPlot(RTIPlot): @@ -99,7 +101,7 @@ class SpectralWidthPlot(RTIPlot): data = { 'width': dataOut.data_width } - + data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) return data, {} class SkyMapPlot(Plot): diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index ed93c21..1f2a84a 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -25,6 +25,7 @@ class SpectraPlot(Plot): channelList = [] def setup(self): + self.nplots = len(self.data.channels) self.ncols = int(numpy.sqrt(self.nplots) + 0.9) self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9) @@ -37,16 +38,18 @@ class SpectraPlot(Plot): self.width = 3.5 * self.ncols self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08}) self.ylabel = 'Range [km]' + + def update_list(self,dataOut): if len(self.channelList) == 0: self.channelList = dataOut.channelList def update(self, dataOut): + self.update_list(dataOut) data = {} meta = {} spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - data['spc'] = spc data['rti'] = dataOut.getPower() data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) @@ -72,7 +75,6 @@ class SpectraPlot(Plot): self.xlabel = "Velocity (m/s)" self.titles = [] - y = self.data.yrange self.y = y @@ -192,8 +194,8 @@ class CrossSpectraPlot(Plot): if ax.firsttime: ax.plt = ax.pcolormesh(x, y, coh.T, - vmin=0, - vmax=1, + vmin=self.zmin_coh, + vmax=self.zmax_coh, cmap=plt.get_cmap(self.colormap_coh) ) else: @@ -210,6 +212,7 @@ class CrossSpectraPlot(Plot): ) else: ax.plt.set_array(phase.T.ravel()) + self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1])) @@ -258,6 +261,7 @@ class RTIPlot(Plot): self.z = self.data[self.CODE] self.z = numpy.array(self.z, dtype=float) self.z = numpy.ma.masked_invalid(self.z) + try: if self.channelList != None: self.titles = ['{} Channel {}'.format( @@ -282,9 +286,10 @@ class RTIPlot(Plot): cmap=plt.get_cmap(self.colormap) ) if self.showprofile: - ax.plot_profile = self.pf_axes[n].plot( - data['rti'][n], self.y)[0] - ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y, + ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0] + + if "noise" in self.data: + ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y, color="k", linestyle="dashed", lw=1)[0] else: ax.collections.remove(ax.collections[0]) @@ -294,8 +299,9 @@ class RTIPlot(Plot): cmap=plt.get_cmap(self.colormap) ) if self.showprofile: - ax.plot_profile.set_data(data['rti'][n], self.y) - ax.plot_noise.set_data(numpy.repeat( + ax.plot_profile.set_data(data[self.CODE][n], self.y) + if "noise" in self.data: + ax.plot_noise.set_data(numpy.repeat( data['noise'][n], len(self.y)), self.y) @@ -447,24 +453,40 @@ class SpectraCutPlot(Plot): CODE = 'spc_cut' plot_type = 'scatter' buffering = False + heights = [] + channelList = [] + maintitle = "Spectra Cuts" def setup(self): self.nplots = len(self.data.channels) self.ncols = int(numpy.sqrt(self.nplots) + 0.9) self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9) - self.width = 3.4 * self.ncols + 1.5 + self.width = 3.6 * self.ncols + 1.5 self.height = 3 * self.nrows self.ylabel = 'Power [dB]' self.colorbar = False self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08}) + if self.selectedHeight: + self.maintitle = "Spectra Cut for %d km " %(int(self.selectedHeight)) def update(self, dataOut): + if len(self.channelList) == 0: + self.channelList = dataOut.channelList + self.heights = dataOut.heightList + if self.selectedHeight: + index_list = numpy.where(self.heights >= self.selectedHeight) + self.height_index = index_list[0] + self.height_index = self.height_index[0] + #print(self.height_index) data = {} meta = {} spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - data['spc'] = spc + if self.selectedHeight: + data['spc'] = spc[:,:,self.height_index] + else: + data['spc'] = spc meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) return data, meta @@ -484,7 +506,7 @@ class SpectraCutPlot(Plot): y = self.data.yrange z = self.data[-1]['spc'] - + #print(z.shape) if self.height_index: index = numpy.array(self.height_index) else: @@ -496,14 +518,21 @@ class SpectraCutPlot(Plot): self.xmin = self.xmin if self.xmin else -self.xmax self.ymin = self.ymin if self.ymin else numpy.nanmin(z) self.ymax = self.ymax if self.ymax else numpy.nanmax(z) - ax.plt = ax.plot(x, z[n, :, index].T) - labels = ['Range = {:2.1f}km'.format(y[i]) for i in index] - self.figures[0].legend(ax.plt, labels, loc='center right') + if self.selectedHeight: + ax.plt = ax.plot(x, z[n,:]) + + else: + ax.plt = ax.plot(x, z[n, :, index].T) + labels = ['Range = {:2.1f}km'.format(y[i]) for i in index] + self.figures[0].legend(ax.plt, labels, loc='center right') else: for i, line in enumerate(ax.plt): - line.set_data(x, z[n, :, index[i]]) - self.titles.append('CH {}'.format(n)) - + if self.selectedHeight: + line.set_data(x, z[n, :]) + else: + line.set_data(x, z[n, :, index[i]]) + self.titles.append('CH {}'.format(self.channelList[n])) + plt.suptitle(self.maintitle) class BeaconPhase(Plot): @@ -735,3 +764,198 @@ class BeaconPhase(Plot): update_figfile=update_figfile) return dataOut + +class NoiselessSpectraPlot(Plot): + ''' + Plot for Spectra data, subtracting + the noise in all channels, using for + amisr-14 data + ''' + + CODE = 'noiseless_spc' + colormap = 'nipy_spectral' + plot_type = 'pcolor' + buffering = False + channelList = [] + + def setup(self): + + self.nplots = len(self.data.channels) + self.ncols = int(numpy.sqrt(self.nplots) + 0.9) + self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9) + self.height = 2.6 * self.nrows + + self.cb_label = 'dB' + if self.showprofile: + self.width = 4 * self.ncols + else: + self.width = 3.5 * self.ncols + self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08}) + self.ylabel = 'Range [km]' + + + def update_list(self,dataOut): + if len(self.channelList) == 0: + self.channelList = dataOut.channelList + + def update(self, dataOut): + + self.update_list(dataOut) + data = {} + meta = {} + n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + (nch, nff, nh) = dataOut.data_spc.shape + n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh)) + noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh)) + #print(noise.shape, "noise", noise) + + spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise + + data['spc'] = spc + data['rti'] = dataOut.getPower() - n1 + + data['noise'] = n0 + meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) + + return data, meta + + def plot(self): + if self.xaxis == "frequency": + x = self.data.xrange[0] + self.xlabel = "Frequency (kHz)" + elif self.xaxis == "time": + x = self.data.xrange[1] + self.xlabel = "Time (ms)" + else: + x = self.data.xrange[2] + self.xlabel = "Velocity (m/s)" + + self.titles = [] + y = self.data.yrange + self.y = y + + data = self.data[-1] + z = data['spc'] + + for n, ax in enumerate(self.axes): + noise = data['noise'][n] + + if ax.firsttime: + self.xmax = self.xmax if self.xmax else numpy.nanmax(x) + self.xmin = self.xmin if self.xmin else -self.xmax + self.zmin = self.zmin if self.zmin else numpy.nanmin(z) + self.zmax = self.zmax if self.zmax else numpy.nanmax(z) + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) + + if self.showprofile: + ax.plt_profile = self.pf_axes[n].plot( + data['rti'][n], y)[0] + ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y, + color="k", linestyle="dashed", lw=1)[0] + + else: + ax.plt.set_array(z[n].T.ravel()) + if self.showprofile: + ax.plt_profile.set_data(data['rti'][n], y) + ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y) + + self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise)) + + +class NoiselessRTIPlot(Plot): + ''' + Plot for RTI data + ''' + + CODE = 'noiseless_rti' + colormap = 'jet' + plot_type = 'pcolorbuffer' + titles = None + channelList = [] + + def setup(self): + self.xaxis = 'time' + self.ncols = 1 + #print("dataChannels ",self.data.channels) + self.nrows = len(self.data.channels) + self.nplots = len(self.data.channels) + self.ylabel = 'Range [km]' + self.xlabel = 'Time' + self.cb_label = 'dB' + self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95}) + self.titles = ['{} Channel {}'.format( + self.CODE.upper(), x) for x in range(self.nplots)] + + def update_list(self,dataOut): + + self.channelList = dataOut.channelList + + + def update(self, dataOut): + if len(self.channelList) == 0: + self.update_list(dataOut) + data = {} + meta = {} + + n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + (nch, nff, nh) = dataOut.data_spc.shape + noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh)) + + + data['noiseless_rti'] = dataOut.getPower() - noise + data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + return data, meta + + def plot(self): + + self.x = self.data.times + self.y = self.data.yrange + self.z = self.data['noiseless_rti'] + self.z = numpy.array(self.z, dtype=float) + self.z = numpy.ma.masked_invalid(self.z) + + try: + if self.channelList != None: + self.titles = ['{} Channel {}'.format( + self.CODE.upper(), x) for x in self.channelList] + except: + if self.channelList.any() != None: + self.titles = ['{} Channel {}'.format( + self.CODE.upper(), x) for x in self.channelList] + if self.decimation is None: + x, y, z = self.fill_gaps(self.x, self.y, self.z) + else: + x, y, z = self.fill_gaps(*self.decimate()) + dummy_var = self.axes #Extrañamente esto actualiza el valor axes + for n, ax in enumerate(self.axes): + self.zmin = self.zmin if self.zmin else numpy.min(self.z) + self.zmax = self.zmax if self.zmax else numpy.max(self.z) + data = self.data[-1] + if ax.firsttime: + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) + if self.showprofile: + ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0] + + if "noise" in self.data: + ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y, + color="k", linestyle="dashed", lw=1)[0] + else: + ax.collections.remove(ax.collections[0]) + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) + if self.showprofile: + ax.plot_profile.set_data(data['noiseless_rti'][n], self.y) + if "noise" in self.data: + ax.plot_noise.set_data(numpy.repeat( + data['noise'][n], len(self.y)), self.y) diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py index 10cee2c..1ff1b65 100644 --- a/schainpy/model/io/jroIO_base.py +++ b/schainpy/model/io/jroIO_base.py @@ -681,7 +681,10 @@ class Reader(object): self.filename = filename self.fileSize = os.path.getsize(filename) - self.fp = self.open_file(filename, self.open_mode) + try: + self.fp = self.open_file(filename, self.open_mode) + except Exception as e: + raise schainpy.admin.SchainError("[Reading] Error in {} file, unable to open".format(filename)) self.flagIsNewFile = 1 return 1 @@ -691,7 +694,7 @@ class Reader(object): """Check if the given datetime is in range""" startDateTime= datetime.datetime.combine(startDate,startTime) endDateTime = datetime.datetime.combine(endDate,endTime) - + if startDateTime <= dt <= endDateTime: return True return False diff --git a/schainpy/model/io/jroIO_kamisr.py b/schainpy/model/io/jroIO_kamisr.py index 2e209f0..db52c38 100644 --- a/schainpy/model/io/jroIO_kamisr.py +++ b/schainpy/model/io/jroIO_kamisr.py @@ -81,7 +81,7 @@ class AMISRReader(ProcessingUnit): #Is really necessary create the output object in the initializer self.dataOut = Voltage() self.dataOut.error=False - + self.margin_days = 1 def setup(self,path=None, startDate=None, @@ -95,7 +95,8 @@ class AMISRReader(ProcessingUnit): nCode = 0, nBaud = 0, online=False, - old_beams=False): + old_beams=False, + margin_days=1): @@ -106,7 +107,7 @@ class AMISRReader(ProcessingUnit): self.code = code self.nCode = int(nCode) self.nBaud = int(nBaud) - + self.margin_days = margin_days #self.findFiles() @@ -195,7 +196,7 @@ class AMISRReader(ProcessingUnit): self.__nSamples = self.nsa self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000 - + #print("amisr-ipp:",self.ippSeconds, self.__ippKm) #for now until understand why the code saved is different (code included even though code not in tuf file) #self.__codeType = 0 # self.__nCode = None @@ -248,7 +249,7 @@ class AMISRReader(ProcessingUnit): dom = int(amisr_dirname_format[6:8]) thisDate = datetime.date(year,month,dom) #margen de un día extra, igual luego se filtra for fecha y hora - if (thisDate>=(self.startDate - datetime.timedelta(days=1)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)): + if (thisDate>=(self.startDate - datetime.timedelta(days=self.margin_days)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)): return amisr_dirname_format except: return None @@ -502,7 +503,6 @@ class AMISRReader(ProcessingUnit): nCode=self.__nCode, nBaud=self.__nBaud, code = self.__code, fClock=1) - #fill system header self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples, nProfiles=self.newProfiles, diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index e74e7ac..0d2d818 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -415,7 +415,7 @@ class HDFWrite(Operation): dataAux = getattr(self.dataOut, self.dataList[i]) dsDict['variable'] = self.dataList[i] else: - log.warning('Attribute {} not found in dataOut', self.name) + log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name) continue if dataAux is None: diff --git a/schainpy/model/proc/jroproc_base.py b/schainpy/model/proc/jroproc_base.py index f98e4fa..a4aa194 100644 --- a/schainpy/model/proc/jroproc_base.py +++ b/schainpy/model/proc/jroproc_base.py @@ -57,15 +57,14 @@ class ProcessingUnit(object): ''' try: - # if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error: - # return self.dataIn.isReady() - #dataIn=None es unidades de Lectura, segunda parte unidades de procesamiento - if self.dataIn is None or (not self.dataIn.error and not self.dataIn.flagNoData): + + if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error: + return self.dataIn.isReady() + elif self.dataIn is None or not self.dataIn.error: self.run(**kwargs) elif self.dataIn.error: self.dataOut.error = self.dataIn.error self.dataOut.flagNoData = True - except: err = traceback.format_exc() diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 94764ca..4d8f0af 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -99,7 +99,7 @@ class ParametersProc(ProcessingUnit): self.dataOut.elevationList = self.dataIn.elevationList def run(self): - + # print("run parameter proc") #---------------------- Voltage Data --------------------------- if self.dataIn.type == "Voltage": @@ -136,11 +136,13 @@ class ParametersProc(ProcessingUnit): self.dataOut.nIncohInt = self.dataIn.nIncohInt self.dataOut.nFFTPoints = self.dataIn.nFFTPoints self.dataOut.ippFactor = self.dataIn.ippFactor + self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() + self.dataOut.ipp = self.dataIn.ipp self.dataOut.abscissaList = self.dataIn.getVelRange(1) self.dataOut.spc_noise = self.dataIn.getNoise() self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1)) # self.dataOut.normFactor = self.dataIn.normFactor - if hasattr(self.dataIn, 'channelList'): + if hasattr(self.dataIn, 'channelList'): self.dataOut.channelList = self.dataIn.channelList if hasattr(self.dataIn, 'pairsList'): self.dataOut.pairsList = self.dataIn.pairsList @@ -308,7 +310,7 @@ class SpectralFilters(Operation): Operation.__init__(self) self.i = 0 - def run(self, dataOut, ): + def run(self, dataOut ): self.spc = dataOut.data_pre[0].copy() self.Num_Chn = self.spc.shape[0] @@ -1470,7 +1472,9 @@ class SpectralFitting(Operation): def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): - + print("run SpectralFitting") + self.dataOut = dataOut.copy() + print(self.dataOut.nIncohInt) if path != None: sys.path.append(path) self.dataOut.library = importlib.import_module(file) @@ -1481,20 +1485,20 @@ class SpectralFitting(Operation): self.dataOut.groupList = groupArray nGroups = groupArray.shape[0] - nChannels = self.dataIn.nChannels - nHeights=self.dataIn.heightList.size + nChannels = self.dataOut.nChannels + nHeights=self.dataOut.heightList.size #Parameters Array self.dataOut.data_param = None #Set constants - constants = self.dataOut.library.setConstants(self.dataIn) + constants = self.dataOut.library.setConstants(self.dataOut) self.dataOut.constants = constants - M = self.dataIn.normFactor - N = self.dataIn.nFFTPoints - ippSeconds = self.dataIn.ippSeconds - K = self.dataIn.nIncohInt - pairsArray = numpy.array(self.dataIn.pairsList) + M = self.dataOut.normFactor + N = self.dataOut.nFFTPoints + ippSeconds = self.dataOut.ippSeconds + K = self.dataOut.nIncohInt + pairsArray = numpy.array(self.dataOut.pairsList) #List of possible combinations listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) @@ -1503,14 +1507,14 @@ class SpectralFitting(Operation): if getSNR: listChannels = groupArray.reshape((groupArray.size)) listChannels.sort() - noise = self.dataIn.getNoise() - self.dataOut.data_snr = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) + noise = self.dataOut.getNoise() + self.dataOut.data_snr = self.__getSNR(self.dataOut.data_spc[listChannels,:,:], noise[listChannels]) for i in range(nGroups): coord = groupArray[i,:] #Input data array - data = self.dataIn.data_spc[coord,:,:]/(M*N) + data = self.dataOut.data_spc[coord,:,:]/(M*N) data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) #Cross Spectra data array for Covariance Matrixes @@ -1519,7 +1523,7 @@ class SpectralFitting(Operation): pairsSel = numpy.array([coord[x],coord[y]]) indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) ind += 1 - dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N) + dataCross = self.dataOut.data_cspc[indCross,:,:]/(M*N) dataCross = dataCross**2/K for h in range(nHeights): @@ -1548,13 +1552,13 @@ class SpectralFitting(Operation): dp = numpy.dot(LT,d) #Initial values - data_spc = self.dataIn.data_spc[coord,:,h] + data_spc = self.dataOut.data_spc[coord,:,h] if (h>0)and(error1[3]<5): p0 = self.dataOut.data_param[i,:,h-1] else: - p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i)) - + #p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i)) + p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants)) try: #Least Squares minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) @@ -1563,7 +1567,8 @@ class SpectralFitting(Operation): error0 = numpy.sum(infodict['fvec']**2)/(2*N) #Error with Jacobian error1 = self.dataOut.library.errorFunction(minp,constants,LT) - except: + except Exception as e: + print("Param error: ",e) minp = p0*numpy.nan error0 = numpy.nan error1 = p0*numpy.nan @@ -1575,7 +1580,11 @@ class SpectralFitting(Operation): self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) self.dataOut.data_param[i,:,h] = minp - return + + + print("SpectralFitting Done") + dataOut = self.dataOut + return dataOut def __residFunction(self, p, dp, LT, constants): @@ -2320,8 +2329,8 @@ class WindProfiler(Operation): class EWDriftsEstimation(Operation): - def __init__(self): - Operation.__init__(self) + # def __init__(self): + # Operation.__init__(self) def __correctValues(self, heiRang, phi, velRadial, SNR): listPhi = phi.tolist() @@ -2338,25 +2347,27 @@ class EWDriftsEstimation(Operation): velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) SNR1 = numpy.zeros([len(phi),len(heiRang1)]) - - for i in rango: - x = heiRang*math.cos(phi[i]) - y1 = velRadial[i,:] - f1 = interpolate.interp1d(x,y1,kind = 'cubic') - - x1 = heiRang1 - y11 = f1(x1) - - y2 = SNR[i,:] - f2 = interpolate.interp1d(x,y2,kind = 'cubic') - y21 = f2(x1) - - velRadial1[i,:] = y11 - SNR1[i,:] = y21 - + try: + for i in rango: + x = heiRang*math.cos(phi[i]) + y1 = velRadial[i,:] + f1 = interpolate.interp1d(x,y1,kind = 'cubic') + + x1 = heiRang1 + y11 = f1(x1) + + y2 = SNR[i,:] + f2 = interpolate.interp1d(x,y2,kind = 'cubic') + y21 = f2(x1) + + velRadial1[i,:] = y11 + SNR1[i,:] = y21 + except Exception as e: + print("wrong values!",e) return heiRang1, velRadial1, SNR1 def run(self, dataOut, zenith, zenithCorrection): + print("run EWDriftsEstimation") heiRang = dataOut.heightList velRadial = dataOut.data_param[:,3,:] SNR = dataOut.data_snr @@ -2366,7 +2377,7 @@ class EWDriftsEstimation(Operation): zenith *= numpy.pi/180 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) - + print("__correctValues Done ") alp = zenith[0] bet = zenith[1] @@ -2384,6 +2395,8 @@ class EWDriftsEstimation(Operation): dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval + + print("EWDriftsEstimation Done ") return #--------------- Non Specular Meteor ---------------- diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 24a56de..ff2e730 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -68,7 +68,6 @@ class SpectraProc(ProcessingUnit): self.dataOut.elevationList = self.dataIn.elevationList - def __getFft(self): """ Convierte valores de Voltaje a Spectra diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index acea65e..7cd3b7a 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -178,7 +178,7 @@ class selectHeights(Operation): 1 si el metodo se ejecuto con exito caso contrario devuelve 0 """ - dataOut = dataOut + self.dataOut = dataOut if minHei and maxHei: @@ -207,7 +207,7 @@ class selectHeights(Operation): self.selectHeightsByIndex(minIndex, maxIndex) - return self.dataOut + return dataOut def selectHeightsByIndex(self, minIndex, maxIndex): """ @@ -1636,3 +1636,123 @@ class PulsePairVoltage(Operation): # self.__startIndex += self.__newNSamples # # return +class SSheightProfiles(Operation): + + step = None + nsamples = None + bufferShape = None + profileShape = None + sshProfiles = None + profileIndex = None + + def __init__(self, **kwargs): + + Operation.__init__(self, **kwargs) + self.isConfig = False + + def setup(self,dataOut ,step = None , nsamples = None): + + if step == None and nsamples == None: + raise ValueError("step or nheights should be specified ...") + + self.step = step + self.nsamples = nsamples + self.__nChannels = dataOut.nChannels + self.__nProfiles = dataOut.nProfiles + self.__nHeis = dataOut.nHeights + shape = dataOut.data.shape #nchannels, nprofiles, nsamples + + residue = (shape[1] - self.nsamples) % self.step + if residue != 0: + 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)) + + deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + numberProfile = self.nsamples + numberSamples = (shape[1] - self.nsamples)/self.step + + self.bufferShape = int(shape[0]), int(numberSamples), int(numberProfile) # nchannels, nsamples , nprofiles + self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples) # nchannels, nprofiles, nsamples + + self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex) + self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex) + + def run(self, dataOut, step, nsamples, code = None, repeat = None): + dataOut.flagNoData = True + + profileIndex = None + #print(dataOut.getFreqRange(1)/1000.) + #exit(1) + if dataOut.flagDataAsBlock: + dataOut.data = numpy.average(dataOut.data,axis=1) + #print("jee") + dataOut.flagDataAsBlock = False + if not self.isConfig: + self.setup(dataOut, step=step , nsamples=nsamples) + #print("Setup done") + self.isConfig = True + + #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]) + #DC_Hae = numpy.array([ 0.001025 +0.0516375j, 0.03485 +0.20923125j, -0.168 -0.02720625j, + #-0.1105375 +0.0707125j, -0.20309375-0.09670625j, 0.189775 +0.02716875j])*(-3.5) + + #DC_Hae = numpy.array([ -32.26 +8.66j, -32.26 +8.66j]) + + #DC_Hae = numpy.array([-2.78500000e-01 -1.39175j, -6.63237294e+02+210.4268625j]) + + #print(dataOut.data[0,13:15]) + #dataOut.data = dataOut.data - DC_Hae[:,None] + #print(dataOut.data[0,13:15]) + #exit(1) + if code is not None: + code = numpy.array(code) + code_block = code + ''' + roll = 0 + code = numpy.roll(code,roll,axis=0) + code = numpy.reshape(code,(5,100,64)) + block = dataOut.CurrentBlock%5 + + day_dif = 0 #day_19_Oct_2021: 3 + code_block = code[block-1-day_dif,:,:] + ''' + if repeat is not None: + code_block = numpy.repeat(code_block, repeats=repeat, axis=1) + #print(code_block.shape) + for i in range(self.buffer.shape[1]): + + if code is not None: + self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block + + else: + + self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:] + + #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]) + + for j in range(self.buffer.shape[0]): + self.sshProfiles[j] = numpy.transpose(self.buffer[j]) + + profileIndex = self.nsamples + deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + ippSeconds = (deltaHeight*1.0e-6)/(0.15) + #print("ippSeconds, dH: ",ippSeconds,deltaHeight) + try: + if dataOut.concat_m is not None: + ippSeconds= ippSeconds/float(dataOut.concat_m) + #print "Profile concat %d"%dataOut.concat_m + except: + pass + + dataOut.data = self.sshProfiles + dataOut.flagNoData = False + dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0] + dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples) + + dataOut.profileIndex = profileIndex + dataOut.flagDataAsBlock = True + dataOut.ippSeconds = ippSeconds + dataOut.step = self.step + #print(numpy.shape(dataOut.data)) + #exit(1) + + return dataOut diff --git a/schainpy/model/proc/modelSpectralFitting.py b/schainpy/model/proc/modelSpectralFitting.py new file mode 100644 index 0000000..2460d47 --- /dev/null +++ b/schainpy/model/proc/modelSpectralFitting.py @@ -0,0 +1,132 @@ +import numpy + +def setConstants(dataOut): + dictionary = {} + dictionary["M"] = dataOut.normFactor + dictionary["N"] = dataOut.nFFTPoints + dictionary["ippSeconds"] = dataOut.ippSeconds + dictionary["K"] = dataOut.nIncohInt + + return dictionary + +def initialValuesFunction(data_spc, constants): + #Constants + M = constants["M"] + N = constants["N"] + ippSeconds = constants["ippSeconds"] + + S1 = data_spc[0,:]/(N*M) + S2 = data_spc[1,:]/(N*M) + + Nl=min(S1) + A=sum(S1-Nl)/N + #x = dataOut.getVelRange() #below matches Madrigal data better + x=numpy.linspace(-(N/2)/(N*ippSeconds),(N/2-1)/(N*ippSeconds),N)*-(6.0/2) + v=sum(x*(S1-Nl))/sum(S1-Nl) + al1=numpy.sqrt(sum(x**2*(S1-Nl))/sum(S2-Nl)-v**2) + p0=[al1,A,A,v,min(S1),min(S2)]#first guess(width,amplitude,velocity,noise) + return p0 + +def modelFunction(p, constants): + ippSeconds = constants["ippSeconds"] + N = constants["N"] + + fm_c = ACFtoSPC(p, constants) + fm = numpy.hstack((fm_c[0],fm_c[1])) + return fm + +def errorFunction(p, constants, LT): + + J=makeJacobian(p, constants) + J =numpy.dot(LT,J) + covm =numpy.linalg.inv(numpy.dot(J.T ,J)) + #calculate error as the square root of the covariance matrix diagonal + #multiplying by 1.96 would give 95% confidence interval + err =numpy.sqrt(numpy.diag(covm)) + return err + +#----------------------------------------------------------------------------------- + +def ACFw(alpha,A1,A2,vd,x,N,ippSeconds): + #creates weighted autocorrelation function based on the operational model + #x is n or N-n + k=2*numpy.pi/3.0 + pdt=x*ippSeconds + #both correlated channels ACFs are created at the sametime + R1=A1*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2) + R2=A2*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2) + # T is the triangle weigthing function + T=1-abs(x)/N + Rp1=T*R1 + Rp2=T*R2 + return [Rp1,Rp2] + +def ACFtoSPC(p, constants): + #calls the create ACF function and transforms the ACF to spectra + N = constants["N"] + ippSeconds = constants["ippSeconds"] + + n=numpy.linspace(0,(N-1),N) + Nn=N-n + R = ACFw(p[0],p[1],p[2],p[3],n,N,ippSeconds) + RN = ACFw(p[0],p[1],p[2],p[3],Nn,N,ippSeconds) + Rf1=R[0]+numpy.conjugate(RN[0]) + Rf2=R[1]+numpy.conjugate(RN[1]) + sw1=numpy.fft.fft(Rf1,n=N) + sw2=numpy.fft.fft(Rf2,n=N) + #the fft needs to be shifted, noise added, and takes only the real part + sw0=numpy.real(numpy.fft.fftshift(sw1))+abs(p[4]) + sw1=numpy.real(numpy.fft.fftshift(sw2))+abs(p[5]) + return [sw0,sw1] + +def makeJacobian(p, constants): + #create Jacobian matrix + N = constants["N"] + IPPt = constants["ippSeconds"] + + n=numpy.linspace(0,(N-1),N) + Nn=N-n + k=2*numpy.pi/3.0 + #created weighted ACF + R=ACFw(p[0],p[1],p[2],p[3],n,N,IPPt) + RN=ACFw(p[0],p[1],p[2],p[3],Nn,N,IPPt) + #take derivatives with respect to the fit parameters + 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)) + 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)) + JA1=R[0]/p[1]+numpy.conjugate(RN[0]/p[1]) + JA2=R[1]/p[2]+numpy.conjugate(RN[1]/p[2]) + Jvd1=R[0]*-1j*k*n*IPPt+numpy.conjugate(RN[0]*-1j*k*Nn*IPPt) + Jvd2=R[1]*-1j*k*n*IPPt+numpy.conjugate(RN[1]*-1j*k*Nn*IPPt) + #fft + sJalp1=numpy.fft.fft(Jalpha1,n=N) + sJalp2=numpy.fft.fft(Jalpha2,n=N) + sJA1=numpy.fft.fft(JA1,n=N) + sJA2=numpy.fft.fft(JA2,n=N) + sJvd1=numpy.fft.fft(Jvd1,n=N) + sJvd2=numpy.fft.fft(Jvd2,n=N) + sJalp1=numpy.real(numpy.fft.fftshift(sJalp1)) + sJalp2=numpy.real(numpy.fft.fftshift(sJalp2)) + sJA1=numpy.real(numpy.fft.fftshift(sJA1)) + sJA2=numpy.real(numpy.fft.fftshift(sJA2)) + sJvd1=numpy.real(numpy.fft.fftshift(sJvd1)) + sJvd2=numpy.real(numpy.fft.fftshift(sJvd2)) + sJnoise=numpy.ones(numpy.shape(sJvd1)) + #combine arrays + za=numpy.zeros([N]) + sJalp=zip(sJalp1,sJalp2) + sJA1=zip(sJA1,za) + sJA2=zip(za,sJA2) + sJvd=zip(sJvd1,sJvd2) + sJn1=zip(sJnoise, za) + sJn2=zip(za, sJnoise) + #reshape from 2D to 1D + sJalp=numpy.reshape(list(sJalp), [2*N]) + sJA1=numpy.reshape(list(sJA1), [2*N]) + sJA2=numpy.reshape(list(sJA2), [2*N]) + sJvd=numpy.reshape(list(sJvd), [2*N]) + sJn1=numpy.reshape(list(sJn1), [2*N]) + sJn2=numpy.reshape(list(sJn2), [2*N]) + #combine into matrix and transpose + J=numpy.array([sJalp,sJA1,sJA2,sJvd,sJn1,sJn2]) + J=J.T + return J