diff --git a/schainpy/model/graphics/jroplot_voltage_lags.py b/schainpy/model/graphics/jroplot_voltage_lags.py index 7ebf8c9..5196da1 100644 --- a/schainpy/model/graphics/jroplot_voltage_lags.py +++ b/schainpy/model/graphics/jroplot_voltage_lags.py @@ -756,6 +756,89 @@ class EDensityPlot(Plot): plt.legend(loc='upper left',fontsize=8.5) #plt.legend(loc='lower left',fontsize=8.5) +class RelativeDenPlot(Plot): + ''' + Written by R. Flores + ''' + ''' + Plot for electron density + ''' + + CODE = 'den' + #plot_name = 'Electron Density' + plot_type = 'scatterbuffer' + + def setup(self): + + self.ncols = 1 + self.nrows = 1 + self.nplots = 1 + self.ylabel = 'Range [km]' + self.xlabel = r'$\mathrm{N_e}$ Relative Electron Density ($\mathrm{1/cm^3}$)' + self.titles = ['Electron Density'] + self.width = 3.5 + self.height = 5.5 + self.colorbar = False + self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1}) + + def update(self, dataOut): + data = {} + meta = {} + + data['den_power'] = dataOut.ph2 + data['den_error'] = dataOut.sdp2 + + meta['yrange'] = dataOut.heightList + + return data, meta + + def plot(self): + + y = self.data.yrange + + ax = self.axes[0] + + data = self.data[-1] + + DenPow = data['den_power'] + errDenPow = data['den_error'] + + if ax.firsttime: + self.autoxticks=False + ax.errorbar(DenPow, y, fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-') + + plt.legend(loc='upper left',fontsize=8.5) + #plt.legend(loc='lower left',fontsize=8.5) + ax.set_xscale("log", nonposx='clip') + grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50) + self.ystep_given=100 + ax.set_yticks(grid_y_ticks,minor=True) + locmaj = LogLocator(base=10,numticks=12) + ax.xaxis.set_major_locator(locmaj) + locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12) + ax.xaxis.set_minor_locator(locmin) + ax.xaxis.set_minor_formatter(NullFormatter()) + ax.grid(which='minor') + + else: + dataBefore = self.data[-2] + DenPowBefore = dataBefore['den_power'] + self.clear_figures() + ax.errorbar(DenPow, y, fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2,linestyle='-') + ax.errorbar(DenPowBefore, y, elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed") + + ax.set_xscale("log", nonposx='clip') + grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50) + ax.set_yticks(grid_y_ticks,minor=True) + locmaj = LogLocator(base=10,numticks=12) + ax.xaxis.set_major_locator(locmaj) + locmin = LogLocator(base=10.0,subs=(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),numticks=12) + ax.xaxis.set_minor_locator(locmin) + ax.xaxis.set_minor_formatter(NullFormatter()) + ax.grid(which='minor') + plt.legend(loc='upper left',fontsize=8.5) + #plt.legend(loc='lower left',fontsize=8.5) + class FaradayAnglePlot(Plot): ''' Written by R. Flores diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index 9c37353..d1fc98c 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -360,13 +360,14 @@ class HDFWriter(Operation): Operation.__init__(self) return - def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None): + def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None, uniqueChannel=False): self.path = path self.blocksPerFile = blocksPerFile self.metadataList = metadataList self.dataList = [s.strip() for s in dataList] self.setType = setType self.description = description + self.uniqueChannel = uniqueChannel if self.metadataList is None: self.metadataList = self.dataOut.metadata_list @@ -388,6 +389,9 @@ class HDFWriter(Operation): elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)): dsDict['nDim'] = 0 else: + if uniqueChannel: #Creates extra dimension to avoid the creation of multiple channels + dataAux = numpy.expand_dims(dataAux, axis=0) + dsDict['nDim'] = len(dataAux.shape) dsDict['shape'] = dataAux.shape dsDict['dsNumber'] = dataAux.shape[0] @@ -422,18 +426,19 @@ class HDFWriter(Operation): return False def run(self, dataOut, path, blocksPerFile=10, metadataList=None, - dataList=[], setType=None, description={}): + dataList=[], setType=None, description={}, uniqueChannel= False): self.dataOut = dataOut if not(self.isConfig): self.setup(path=path, blocksPerFile=blocksPerFile, metadataList=metadataList, dataList=dataList, - setType=setType, description=description) + setType=setType, description=description, uniqueChannel=uniqueChannel) self.isConfig = True self.setNextFile() self.putData() + return def setNextFile(self): @@ -490,7 +495,7 @@ class HDFWriter(Operation): self.writeData(self.fp) def getLabel(self, name, x=None): - + #print("x: ", x) if x is None: if 'Data' in self.description: data = self.description['Data'] @@ -558,7 +563,7 @@ class HDFWriter(Operation): dtsets = [] data = [] - + #print("self.dsList: ", self.dsList) for dsInfo in self.dsList: if dsInfo['nDim'] == 0: ds = grp.create_dataset( @@ -582,6 +587,11 @@ class HDFWriter(Operation): dtype=dsInfo['dtype']) dtsets.append(ds) data.append((dsInfo['variable'], i)) + + if self.uniqueChannel: #Deletes extra dimension created to avoid the creation of multiple channels + dataAux = getattr(self.dataOut, dsInfo['variable']) + dataAux = dataAux[0] + fp.flush() log.log('Creating file: {}'.format(fp.filename), self.name) diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 35de02c..2030c14 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -779,6 +779,73 @@ class removeInterference(Operation): return self.dataOut +class removeInterferenceAtFreq(Operation): + ''' + Written by R. Flores + ''' + """Operation to remove interfernce at a known frequency(s). + + Parameters: + ----------- + None + + Example + -------- + + op = proc_unit.addOperation(name='removeInterferenceAtFreq') + + """ + + def __init__(self): + + Operation.__init__(self) + + def run(self, dataOut, freq = None, freqList = None): + + VelRange = dataOut.getVelRange() + #print("VelRange: ", VelRange) + + freq_ids = [] + + if freq is not None: + #print("freq") + #if freq < 0: + inda = numpy.where(VelRange >= freq) + minIndex = inda[0][0] + #print(numpy.shape(dataOut.dataLag_spc)) + dataOut.data_spc[:,minIndex,:] = numpy.nan + + #inda = numpy.where(VelRange >= ymin_noise) + #indb = numpy.where(VelRange <= ymax_noise) + + #minIndex = inda[0][0] + #maxIndex = indb[0][-1] + + elif freqList is not None: + #print("freqList") + for freq in freqList: + #if freq < 0: + inda = numpy.where(VelRange >= freq) + minIndex = inda[0][0] + #print(numpy.shape(dataOut.dataLag_spc)) + if freq > 0: + #dataOut.data_spc[:,minIndex-1,:] = numpy.nan + freq_ids.append(minIndex-1) + else: + #dataOut.data_spc[:,minIndex,:] = numpy.nan + freq_ids.append(minIndex) + else: + raise ValueError("freq or freqList should be specified ...") + + #freq_ids = numpy.array(freq_ids).flatten() + + avg = numpy.mean(dataOut.data_spc[:,[t for t in range(dataOut.data_spc.shape[0]) if t not in freq_ids],:],axis=1) + + for p in list(freq_ids): + dataOut.data_spc[:,p,:] = avg#numpy.nan + + + return dataOut class IncohInt(Operation): diff --git a/schainpy/model/proc/jroproc_spectra_lags_faraday.py b/schainpy/model/proc/jroproc_spectra_lags_faraday.py index 4a02f6a..7a1f75b 100644 --- a/schainpy/model/proc/jroproc_spectra_lags_faraday.py +++ b/schainpy/model/proc/jroproc_spectra_lags_faraday.py @@ -2877,8 +2877,8 @@ class IntegrationFaradaySpectraNoLags(Operation): self.__lastdatatime = 0 self.__buffer_spc = [] - #self.__buffer_cspc = [] - self.__buffer_cspc = None + self.__buffer_cspc = [] + #self.__buffer_cspc = None self.__buffer_dc = 0 self.__profIndex = 0 @@ -3210,9 +3210,9 @@ class IntegrationFaradaySpectraNoLags(Operation): #print(self.__buffer_spc[:,1,3,20,0]) #print(self.__buffer_spc[:,1,5,37,0]) data_spc = numpy.sum(self.__buffer_spc,axis=0) - print("data_spc: ", data_spc[0,:,0]) - print("data_spc: ", data_spc[0,:,7]) - print("shape: ", numpy.shape(data_spc)) + #print("data_spc: ", data_spc[0,:,0]) + #print("data_spc: ", data_spc[0,:,7]) + #print("shape: ", numpy.shape(data_spc)) #exit(1) #data_cspc = numpy.sum(self.__buffer_cspc,axis=0) if self.__buffer_cspc is not None: @@ -3228,8 +3228,8 @@ class IntegrationFaradaySpectraNoLags(Operation): n = self.__profIndex self.__buffer_spc = [] - #self.__buffer_cspc = [] - self.__buffer_cspc = None + self.__buffer_cspc = [] + #self.__buffer_cspc = None self.__buffer_dc = 0 self.__profIndex = 0 @@ -3292,12 +3292,13 @@ class IntegrationFaradaySpectraNoLags(Operation): #print(numpy.shape(dataOut.data_spc)) #print(numpy.shape(dataOut.data_cspc)) #exit(1) - dataOut.data_cspc = None + #dataOut.data_cspc = None if not self.isConfig: self.setup(dataOut, n, timeInterval, overlapping,DPL ) self.isConfig = True if not self.ByLags: + #print("dataOut.data_cspc: ", dataOut.data_cspc) self.nProfiles=dataOut.nProfiles self.nChannels=dataOut.nChannels self.nHeights=dataOut.nHeights @@ -3320,6 +3321,7 @@ class IntegrationFaradaySpectraNoLags(Operation): dataOut.data_spc = numpy.squeeze(avgdata_spc) dataOut.data_cspc = numpy.squeeze(avgdata_cspc) + dataOut.data_cspc = numpy.expand_dims(dataOut.data_cspc, axis=0) dataOut.data_dc = avgdata_dc else: dataOut.dataLag_spc = avgdata_spc @@ -4290,7 +4292,7 @@ class SpectraDataToFaraday_MST(Operation): #MST MODE dataOut.pan = dataOut.tnoise[0] dataOut.pbn = dataOut.tnoise[1] - def noise(self,dataOut): + def noise(self,dataOut,minIndex,maxIndex): dataOut.noise_lag = numpy.zeros((dataOut.nChannels),'float32') #print("Lags") @@ -4312,9 +4314,9 @@ class SpectraDataToFaraday_MST(Operation): #MST MODE #dataOut.noise_lag[1] = dataOut.getNoise(ymin_index=80,ymax_index=106)[1] if dataOut.flagDecodeData: #dataOut.noise_lag[1] = dataOut.getNoise(ymin_index=150,ymax_index=200)[1] - dataOut.noise_lag[1] = dataOut.getNoise()[1] + dataOut.noise_lag[1] = dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[1] else: - dataOut.noise_lag[1] = dataOut.getNoise()[1] + dataOut.noise_lag[1] = dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[1] #else: #dataOut.noise_lag[1,lag] = numpy.mean(dataOut.noise_lag[1,:6]) #dataOut.noise_lag[:,lag] = dataOut.getNoise(ymin_index=33,ymax_index=46) @@ -4324,9 +4326,9 @@ class SpectraDataToFaraday_MST(Operation): #MST MODE dataOut.data_spc = dataOut.dataLag_spc[:,:,:,0] if dataOut.flagDecodeData: #dataOut.noise_lag[0] = dataOut.getNoise(ymin_index=150,ymax_index=200)[0] - dataOut.noise_lag[0] = dataOut.getNoise()[0] + dataOut.noise_lag[0] = dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[0] else: - dataOut.noise_lag[0] = dataOut.getNoise()[0] + dataOut.noise_lag[0] = dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[0] dataOut.tnoise = dataOut.noise_lag/float(dataOut.nProfiles*dataOut.nIncohInt) #dataOut.tnoise /= float(dataOut.nProfiles*dataOut.nIncohInt) @@ -4392,12 +4394,7 @@ class SpectraDataToFaraday_MST(Operation): #MST MODE input() ''' - - - - - - def run(self,dataOut): + def run(self,dataOut,ymin_noise = None,ymax_noise = None): dataOut.paramInterval=0#int(dataOut.nint*dataOut.header[7][0]*2 ) dataOut.lat=-11.95 @@ -4421,28 +4418,52 @@ class SpectraDataToFaraday_MST(Operation): #MST MODE dataOut.data_spc /= dataOut.windowOfFilter dataOut.data_cspc /= dataOut.windowOfFilter ''' - #print(dataOut.data_spc.shape) - print("*****************Sum: ", numpy.sum(dataOut.data_spc[0])) - print("*******************normFactor: *******************", dataOut.normFactor) + #print("dataOut.data_spc.shape: ", dataOut.data_spc.shape) + #print("dataOut.data_cspc.shape: ", dataOut.data_cspc.shape) + #print("*****************Sum: ", numpy.sum(dataOut.data_spc[0])) + #print("*******************normFactor: *******************", dataOut.normFactor) dataOut.dataLag_spc = numpy.stack((dataOut.data_spc, dataOut.data_spc), axis=-1) dataOut.dataLag_cspc = numpy.stack((dataOut.data_cspc, dataOut.data_cspc), axis=-1) #print(dataOut.dataLag_spc.shape) dataOut.DPL = numpy.shape(dataOut.dataLag_spc)[-1] + #exit(1) self.ConvertData(dataOut) - self.noise(dataOut) + + inda = numpy.where(dataOut.heightList >= ymin_noise) + indb = numpy.where(dataOut.heightList <= ymax_noise) + + minIndex = inda[0][0] + maxIndex = indb[0][-1] + + #print("ymin_noise: ", dataOut.heightList[minIndex]) + #print("ymax_noise: ", dataOut.heightList[maxIndex]) + + self.noise(dataOut,minIndex,maxIndex) dataOut.NAVG=16#dataOut.rnint2[0] #CHECK THIS! dataOut.MAXNRANGENDT=dataOut.NDP - ''' - if dataOut.flagDecodeData: - print(dataOut.kabxys_integrated[4][:,0,0]) + #''' + if 0: + #print(dataOut.kabxys_integrated[4][:,0,0]) + #print("dataOut.heightList: ", dataOut.heightList) + #print("dataOut.pbn: ", dataOut.pbn) + print("INSIDE") import matplotlib.pyplot as plt - plt.plot(dataOut.kabxys_integrated[4][:,0,0],dataOut.heightList) - plt.axvline(dataOut.pan) - plt.xlim(1.1*1e3,0.6*1e6) - plt.ylim(30,90) + #print("dataOut.getPower(): ", dataOut.getPower()) + plt.plot(10*numpy.log10(dataOut.kabxys_integrated[4][:,0,0]),dataOut.heightList) + #plt.plot(10**((dataOut.getPower()[1])/10),dataOut.heightList) + #plt.plot(dataOut.getPower()[0],dataOut.heightList) + #plt.plot(dataOut.dataLag_spc[:,:,:,0],dataOut.heightList) + plt.axvline(10*numpy.log10(dataOut.pan)) + #print(dataOut.nProfiles) + #plt.axvline(10*numpy.log10(1*dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[0]/dataOut.normFactor)) + #print("10*numpy.log10(dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[1]/dataOut.normFactor): ", 10*numpy.log10(dataOut.getNoise(ymin_index=minIndex,ymax_index=maxIndex)[1]/dataOut.normFactor)) + #plt.xlim(1,25000) + #plt.xlim(15,20) + #plt.ylim(30,90) + plt.grid() plt.show() - ''' + #''' dataOut.DPL = 1 return dataOut diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 1316e23..566065f 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -3115,6 +3115,8 @@ class DoublePulseACFs(Operation): dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn) ''' #print("init 2.6",pa,dataOut.pan) + #dataOut.pan = 23600/2 + #dataOut.pbn = 23600/2 dataOut.p[i,j]=pa+pb-(dataOut.pan+dataOut.pbn) #print(i,j,dataOut.p[i,j]) dataOut.sdp[i,j]=2*dataOut.rnint2[j]*((pa+pb)*(pa+pb)) @@ -3122,7 +3124,16 @@ class DoublePulseACFs(Operation): rhorp=dataOut.kabxys_integrated[8][i,j,0]+dataOut.kabxys_integrated[11][i,j,0] rhoip=dataOut.kabxys_integrated[10][i,j,0]-dataOut.kabxys_integrated[9][i,j,0] - + ''' + import matplotlib.pyplot as plt + plt.plot(numpy.abs(dataOut.kabxys_integrated[4][:,j,0]+dataOut.kabxys_integrated[5][:,j,0])+numpy.abs(dataOut.kabxys_integrated[6][:,j,0]+dataOut.kabxys_integrated[7][:,j,0]),dataOut.heightList) + plt.axvline((dataOut.pan+dataOut.pbn)) + plt.xlim(20000,30000) + #plt.plot(numpy.abs(dataOut.kabxys_integrated[4][:,j,0]+dataOut.kabxys_integrated[5][:,j,0])+numpy.abs(dataOut.kabxys_integrated[6][:,j,0]+dataOut.kabxys_integrated[7][:,j,0])-(dataOut.pan+dataOut.pbn),dataOut.heightList) + #plt.xlim(1,10000) + plt.grid() + plt.show() + ''' if ((pa>dataOut.pan)&(pb>dataOut.pbn)): ss4=numpy.abs((pa-dataOut.pan)*(pb-dataOut.pbn)) @@ -3206,6 +3217,18 @@ class DoublePulseACFs(Operation): #print("p: ",dataOut.p[33,:]) #exit(1) #''' + ''' + import matplotlib.pyplot as plt + #plt.plot(numpy.abs(dataOut.kabxys_integrated[4][:,j,0]+dataOut.kabxys_integrated[5][:,j,0])+numpy.abs(dataOut.kabxys_integrated[6][:,j,0]+dataOut.kabxys_integrated[7][:,j,0]),dataOut.heightList) + #plt.axvline((dataOut.pan+dataOut.pbn)) + print(numpy.shape(dataOut.p)) + plt.plot(dataOut.p[:,0]*dataOut.heightList*dataOut.heightList,dataOut.heightList) + + #plt.xlim(1,100000000) + plt.xlim(100,100000000) + plt.grid() + plt.show() + ''' #print(numpy.sum(dataOut.rhor)) #exit(1) return dataOut @@ -3442,9 +3465,19 @@ class FaradayAngleAndDPPower(Operation): dataOut.flagTeTiCorrection = False #print("ph2: ", numpy.sum(dataOut.ph2[:16])) #print("ph2: ", numpy.sum(dataOut.ph2[16:32])) + ''' + import matplotlib.pyplot as plt + #plt.plot(numpy.abs(dataOut.kabxys_integrated[4][:,j,0]+dataOut.kabxys_integrated[5][:,j,0])+numpy.abs(dataOut.kabxys_integrated[6][:,j,0]+dataOut.kabxys_integrated[7][:,j,0]),dataOut.heightList) + #plt.axvline((dataOut.pan+dataOut.pbn)) + #print(numpy.shape(dataOut.p)) + plt.plot(dataOut.ph2,dataOut.heightList) + plt.xlim(1000,1000000000) + #plt.ylim(50,400) + plt.grid() + plt.show() #exit(1) - + ''' return dataOut