schain Commit - r1377:5453f456345b · Jicamarca Repository

isr commit

rflores -

r1377:5453f456345b

parent child

Context file:

r1377:5453f456345b

Expand all files

The requested changes are too big and content was truncated. Show full diff

schainpy/model/graphics/jroplot_voltage_lags.py created 644 +1156 0

This diff has been collapsed as it changes many lines, (1156 lines changed) Show them Hide them
		@@ -0,0 +1,1156
	1
	2	import os
	3	import datetime
	4	import numpy
	5	from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
	6
	7	from .jroplot_spectra import RTIPlot, NoisePlot
	8
	9	from schainpy.utils import log
	10	from .plotting_codes import *
	11
	12	from schainpy.model.graphics.jroplot_base import Plot, plt
	13
	14	import matplotlib.pyplot as plt
	15	import matplotlib.colors as colors
	16
	17	import time
	18	import math
	19
	20
	21	from matplotlib.ticker import MultipleLocator
	22
	23
	24
	25	class RTIDPPlot(RTIPlot):
	26
	27	'''
	28	Plot for RTI Double Pulse Experiment
	29	'''
	30
	31	CODE = 'RTIDP'
	32	colormap = 'jro'
	33	plot_name = 'RTI'
	34
	35	#cb_label = 'Ne Electron Density (1/cm3)'
	36
	37	def setup(self):
	38	self.xaxis = 'time'
	39	self.ncols = 1
	40	self.nrows = 3
	41	self.nplots = self.nrows
	42	#self.height=10
	43	if self.showSNR:
	44	self.nrows += 1
	45	self.nplots += 1
	46
	47	self.ylabel = 'Height [km]'
	48	self.xlabel = 'Time (LT)'
	49
	50	self.cb_label = 'Intensity (dB)'
	51
	52
	53	#self.cb_label = cb_label
	54
	55	self.titles = ['{} Channel {}'.format(
	56	self.plot_name.upper(), '0x1'),'{} Channel {}'.format(
	57	self.plot_name.upper(), '0'),'{} Channel {}'.format(
	58	self.plot_name.upper(), '1')]
	59
	60
	61	def plot(self):
	62
	63	self.data.normalize_heights()
	64	self.x = self.data.times
	65	self.y = self.data.heights[0:self.data.NDP]
	66
	67	if self.showSNR:
	68	self.z = numpy.concatenate(
	69	(self.data[self.CODE], self.data['snr'])
	70	)
	71	else:
	72
	73	self.z = self.data[self.CODE]
	74	#print(numpy.max(self.z[0,0:]))
	75
	76	self.z = numpy.ma.masked_invalid(self.z)
	77
	78	if self.decimation is None:
	79	x, y, z = self.fill_gaps(self.x, self.y, self.z)
	80	else:
	81	x, y, z = self.fill_gaps(*self.decimate())
	82
	83	for n, ax in enumerate(self.axes):
	84
	85
	86	self.zmax = self.zmax if self.zmax is not None else numpy.max(
	87	self.z[1][0,12:40])
	88	self.zmin = self.zmin if self.zmin is not None else numpy.min(
	89	self.z[1][0,12:40])
	90
	91
	92
	93	if ax.firsttime:
	94
	95	if self.zlimits is not None:
	96	self.zmin, self.zmax = self.zlimits[n]
	97
	98
	99	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	100	vmin=self.zmin,
	101	vmax=self.zmax,
	102	cmap=self.cmaps[n]
	103	)
	104	#plt.tight_layout()
	105	else:
	106	if self.zlimits is not None:
	107	self.zmin, self.zmax = self.zlimits[n]
	108	ax.collections.remove(ax.collections[0])
	109	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	110	vmin=self.zmin,
	111	vmax=self.zmax,
	112	cmap=self.cmaps[n]
	113	)
	114	#plt.tight_layout()
	115
	116
	117	class RTILPPlot(RTIPlot):
	118
	119	'''
	120	Plot for RTI Long Pulse
	121	'''
	122
	123	CODE = 'RTILP'
	124	colormap = 'jro'
	125	plot_name = 'RTI LP'
	126
	127	#cb_label = 'Ne Electron Density (1/cm3)'
	128
	129	def setup(self):
	130	self.xaxis = 'time'
	131	self.ncols = 1
	132	self.nrows = 4
	133	self.nplots = self.nrows
	134	if self.showSNR:
	135	self.nrows += 1
	136	self.nplots += 1
	137
	138	self.ylabel = 'Height [km]'
	139	self.xlabel = 'Time (LT)'
	140
	141	self.cb_label = 'Intensity (dB)'
	142
	143
	144
	145	#self.cb_label = cb_label
	146
	147	self.titles = ['{} Channel {}'.format(
	148	self.plot_name.upper(), '0'),'{} Channel {}'.format(
	149	self.plot_name.upper(), '1'),'{} Channel {}'.format(
	150	self.plot_name.upper(), '2'),'{} Channel {}'.format(
	151	self.plot_name.upper(), '3')]
	152
	153
	154	def plot(self):
	155
	156	self.data.normalize_heights()
	157	self.x = self.data.times
	158	self.y = self.data.heights[0:self.data.NRANGE]
	159
	160	if self.showSNR:
	161	self.z = numpy.concatenate(
	162	(self.data[self.CODE], self.data['snr'])
	163	)
	164	else:
	165
	166	self.z = self.data[self.CODE]
	167	#print(numpy.max(self.z[0,0:]))
	168
	169	self.z = numpy.ma.masked_invalid(self.z)
	170
	171	if self.decimation is None:
	172	x, y, z = self.fill_gaps(self.x, self.y, self.z)
	173	else:
	174	x, y, z = self.fill_gaps(*self.decimate())
	175
	176	for n, ax in enumerate(self.axes):
	177
	178
	179	self.zmax = self.zmax if self.zmax is not None else numpy.max(
	180	self.z[1][0,12:40])
	181	self.zmin = self.zmin if self.zmin is not None else numpy.min(
	182	self.z[1][0,12:40])
	183
	184	if ax.firsttime:
	185
	186	if self.zlimits is not None:
	187	self.zmin, self.zmax = self.zlimits[n]
	188
	189
	190	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	191	vmin=self.zmin,
	192	vmax=self.zmax,
	193	cmap=self.cmaps[n]
	194	)
	195	#plt.tight_layout()
	196	else:
	197	if self.zlimits is not None:
	198	self.zmin, self.zmax = self.zlimits[n]
	199	ax.collections.remove(ax.collections[0])
	200	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	201	vmin=self.zmin,
	202	vmax=self.zmax,
	203	cmap=self.cmaps[n]
	204	)
	205	#plt.tight_layout()
	206
	207
	208	class DenRTIPlot(RTIPlot):
	209
	210	'''
	211	Plot for Den
	212	'''
	213
	214	CODE = 'denrti'
	215	colormap = 'jro'
	216	plot_name = 'Electron Density'
	217
	218	#cb_label = 'Ne Electron Density (1/cm3)'
	219
	220	def setup(self):
	221	self.xaxis = 'time'
	222	self.ncols = 1
	223	self.nrows = self.data.shape(self.CODE)[0]
	224	self.nplots = self.nrows
	225	if self.showSNR:
	226	self.nrows += 1
	227	self.nplots += 1
	228
	229	self.ylabel = 'Height [km]'
	230	self.xlabel = 'Time (LT)'
	231
	232	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
	233
	234	if self.CODE == 'denrti' or self.CODE=='denrtiLP':
	235	self.cb_label = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
	236
	237	#self.cb_label = cb_label
	238	if not self.titles:
	239	self.titles = self.data.parameters \
	240	if self.data.parameters else ['{}'.format(self.plot_name)]
	241	if self.showSNR:
	242	self.titles.append('SNR')
	243
	244	def plot(self):
	245
	246	self.data.normalize_heights()
	247	self.x = self.data.times
	248	self.y = self.data.heights
	249
	250
	251
	252	if self.showSNR:
	253	self.z = numpy.concatenate(
	254	(self.data[self.CODE], self.data['snr'])
	255	)
	256	else:
	257	self.z = self.data[self.CODE]
	258
	259	self.z = numpy.ma.masked_invalid(self.z)
	260
	261	if self.decimation is None:
	262	x, y, z = self.fill_gaps(self.x, self.y, self.z)
	263	else:
	264	x, y, z = self.fill_gaps(*self.decimate())
	265
	266	for n, ax in enumerate(self.axes):
	267
	268	self.zmax = self.zmax if self.zmax is not None else numpy.max(
	269	self.z[n])
	270	self.zmin = self.zmin if self.zmin is not None else numpy.min(
	271	self.z[n])
	272
	273	if ax.firsttime:
	274
	275	if self.zlimits is not None:
	276	self.zmin, self.zmax = self.zlimits[n]
	277	if numpy.log10(self.zmin)<0:
	278	self.zmin=1
	279	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	280	vmin=self.zmin,
	281	vmax=self.zmax,
	282	cmap=self.cmaps[n],
	283	norm=colors.LogNorm()
	284	)
	285	#plt.tight_layout()
	286
	287	else:
	288	if self.zlimits is not None:
	289	self.zmin, self.zmax = self.zlimits[n]
	290	ax.collections.remove(ax.collections[0])
	291	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	292	vmin=self.zmin,
	293	vmax=self.zmax,
	294	cmap=self.cmaps[n],
	295	norm=colors.LogNorm()
	296	)
	297	#plt.tight_layout()
	298
	299
	300
	301	class DenRTILPPlot(DenRTIPlot):
	302
	303	'''
	304	Plot for Electron Temperature
	305	'''
	306
	307	CODE = 'denrtiLP'
	308	colormap = 'jro'
	309	plot_name = 'Electron Density'
	310
	311
	312	class ETempRTIPlot(RTIPlot):
	313
	314	'''
	315	Plot for Electron Temperature
	316	'''
	317
	318	CODE = 'ETemp'
	319	colormap = 'jet'
	320	plot_name = 'Electron Temperature'
	321
	322	#cb_label = 'Ne Electron Density (1/cm3)'
	323
	324	def setup(self):
	325	self.xaxis = 'time'
	326	self.ncols = 1
	327	self.nrows = self.data.shape(self.CODE)[0]
	328	self.nplots = self.nrows
	329	if self.showSNR:
	330	self.nrows += 1
	331	self.nplots += 1
	332
	333	self.ylabel = 'Height [km]'
	334	self.xlabel = 'Time (LT)'
	335	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
	336	if self.CODE == 'ETemp' or self.CODE == 'ETempLP':
	337	self.cb_label = 'Electron Temperature (K)'
	338	if self.CODE == 'ITemp' or self.CODE == 'ITempLP':
	339	self.cb_label = 'Ion Temperature (K)'
	340
	341
	342	if not self.titles:
	343	self.titles = self.data.parameters \
	344	if self.data.parameters else ['{}'.format(self.plot_name)]
	345	if self.showSNR:
	346	self.titles.append('SNR')
	347
	348	def plot(self):
	349
	350	self.data.normalize_heights()
	351	self.x = self.data.times
	352	self.y = self.data.heights
	353
	354	if self.showSNR:
	355	self.z = numpy.concatenate(
	356	(self.data[self.CODE], self.data['snr'])
	357	)
	358	else:
	359	self.z = self.data[self.CODE]
	360
	361	self.z = numpy.ma.masked_invalid(self.z)
	362
	363	if self.decimation is None:
	364	x, y, z = self.fill_gaps(self.x, self.y, self.z)
	365	else:
	366	x, y, z = self.fill_gaps(*self.decimate())
	367
	368	for n, ax in enumerate(self.axes):
	369
	370	self.zmax = self.zmax if self.zmax is not None else numpy.max(
	371	self.z[n])
	372	self.zmin = self.zmin if self.zmin is not None else numpy.min(
	373	self.z[n])
	374
	375	if ax.firsttime:
	376
	377	if self.zlimits is not None:
	378	self.zmin, self.zmax = self.zlimits[n]
	379
	380	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	381	vmin=self.zmin,
	382	vmax=self.zmax,
	383	cmap=self.cmaps[n]
	384	)
	385	#plt.tight_layout()
	386
	387	else:
	388	if self.zlimits is not None:
	389	self.zmin, self.zmax = self.zlimits[n]
	390	ax.collections.remove(ax.collections[0])
	391	ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
	392	vmin=self.zmin,
	393	vmax=self.zmax,
	394	cmap=self.cmaps[n]
	395	)
	396	#plt.tight_layout()
	397
	398
	399
	400	class ITempRTIPlot(ETempRTIPlot):
	401
	402	'''
	403	Plot for Ion Temperature
	404	'''
	405
	406	CODE = 'ITemp'
	407	colormap = 'jet'
	408	plot_name = 'Ion Temperature'
	409
	410
	411	class ElectronTempLPPlot(ETempRTIPlot):
	412
	413	'''
	414	Plot for Electron Temperature LP
	415	'''
	416
	417	CODE = 'ETempLP'
	418	colormap = 'jet'
	419	plot_name = 'Electron Temperature'
	420
	421
	422	class IonTempLPPlot(ETempRTIPlot):
	423
	424	'''
	425	Plot for Ion Temperature LP
	426	'''
	427
	428	CODE = 'ITempLP'
	429	colormap = 'jet'
	430	plot_name = 'Ion Temperature'
	431
	432
	433	class HFracRTIPlot(ETempRTIPlot):
	434
	435	'''
	436	Plot for H+ LP
	437	'''
	438
	439	CODE = 'HFracLP'
	440	colormap = 'jet'
	441	plot_name = 'H+ Frac'
	442
	443
	444	class HeFracRTIPlot(ETempRTIPlot):
	445
	446	'''
	447	Plot for He+ LP
	448	'''
	449
	450	CODE = 'HeFracLP'
	451	colormap = 'jet'
	452	plot_name = 'He+ Frac'
	453
	454
	455	class TempsDPPlot(Plot):
	456	'''
	457	Plot for Electron - Ion Temperatures
	458	'''
	459
	460	CODE = 'tempsDP'
	461	plot_name = 'Temperatures'
	462	plot_type = 'scatterbuffer'
	463
	464
	465	def setup(self):
	466
	467	self.ncols = 1
	468	self.nrows = 1
	469	self.nplots = 1
	470	self.ylabel = 'Range [km]'
	471	self.xlabel = 'Temperature (K)'
	472	self.width = 3.5
	473	self.height = 5.5
	474	self.colorbar = False
	475	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
	476	if not self.titles:
	477	self.titles = self.data.parameters \
	478	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	479
	480	def plot(self):
	481
	482	self.x = self.data['tempsDP'][:,-1]
	483	self.y = self.data.heights[0:self.data.NSHTS]
	484
	485	self.xmin = -100
	486	self.xmax = 5000
	487	ax = self.axes[0]
	488
	489	if ax.firsttime:
	490
	491	ax.errorbar(self.x, self.y, xerr=self.data.ete2, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
	492	ax.errorbar(self.data.ti2, self.y, fmt='k^', xerr=self.data.eti2,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
	493	plt.legend(loc='lower right')
	494	self.ystep_given = 50
	495	ax.yaxis.set_minor_locator(MultipleLocator(15))
	496	ax.grid(which='minor')
	497	#plt.tight_layout()
	498
	499
	500	else:
	501	self.clear_figures()
	502	ax.errorbar(self.x, self.y, xerr=self.data.ete2, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
	503	ax.errorbar(self.data.ti2, self.y, fmt='k^', xerr=self.data.eti2,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
	504	plt.legend(loc='lower right')
	505	ax.yaxis.set_minor_locator(MultipleLocator(15))
	506	#plt.tight_layout()
	507
	508
	509	class TempsHPPlot(Plot):
	510	'''
	511	Plot for Temperatures Hybrid Experiment
	512	'''
	513
	514	CODE = 'temps_LP'
	515	plot_name = 'Temperatures'
	516	plot_type = 'scatterbuffer'
	517
	518
	519	def setup(self):
	520
	521	self.ncols = 1
	522	self.nrows = 1
	523	self.nplots = 1
	524	self.ylabel = 'Range [km]'
	525	self.xlabel = 'Temperature (K)'
	526	self.width = 3.5
	527	self.height = 6.5
	528	self.colorbar = False
	529	if not self.titles:
	530	self.titles = self.data.parameters \
	531	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	532
	533	def plot(self):
	534
	535	self.x = self.data['temps_LP'][:,-1]
	536	self.y = self.data.heights[0:self.data.NACF]
	537	self.xmin = -100
	538	self.xmax = 4500
	539	ax = self.axes[0]
	540
	541	if ax.firsttime:
	542
	543	ax.errorbar(self.x, self.y, xerr=self.data.ete, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
	544	ax.errorbar(self.data.ti, self.y, fmt='k^', xerr=self.data.eti,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
	545	plt.legend(loc='lower right')
	546	self.ystep_given = 200
	547	ax.yaxis.set_minor_locator(MultipleLocator(15))
	548	ax.grid(which='minor')
	549	#plt.tight_layout()
	550
	551
	552	else:
	553	self.clear_figures()
	554	ax.errorbar(self.x, self.y, xerr=self.data.ete, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
	555	ax.errorbar(self.data.ti, self.y, fmt='k^', xerr=self.data.eti,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
	556	plt.legend(loc='lower right')
	557	ax.yaxis.set_minor_locator(MultipleLocator(15))
	558	#plt.tight_layout()
	559
	560
	561	class FracsHPPlot(Plot):
	562	'''
	563	Plot for Composition LP
	564	'''
	565
	566	CODE = 'fracs_LP'
	567	plot_name = 'Composition'
	568	plot_type = 'scatterbuffer'
	569
	570
	571	def setup(self):
	572
	573	self.ncols = 1
	574	self.nrows = 1
	575	self.nplots = 1
	576	self.ylabel = 'Range [km]'
	577	self.xlabel = 'Frac'
	578	self.width = 3.5
	579	self.height = 6.5
	580	self.colorbar = False
	581	if not self.titles:
	582	self.titles = self.data.parameters \
	583	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	584
	585	def plot(self):
	586
	587	self.x = self.data['fracs_LP'][:,-1]
	588	self.y = self.data.heights[0:self.data.NACF]
	589
	590	self.xmin = 0
	591	self.xmax = 1
	592	ax = self.axes[0]
	593
	594	if ax.firsttime:
	595
	596	ax.errorbar(self.x, self.y[self.data.cut:], xerr=self.data.eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
	597	ax.errorbar(self.data.phe, self.y[self.data.cut:], fmt='k^', xerr=self.data.ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
	598	plt.legend(loc='lower right')
	599	self.xstep_given = 0.2
	600	self.ystep_given = 200
	601	ax.yaxis.set_minor_locator(MultipleLocator(15))
	602	ax.grid(which='minor')
	603	#plt.tight_layout()
	604
	605
	606	else:
	607	self.clear_figures()
	608	ax.errorbar(self.x, self.y[self.data.cut:], xerr=self.data.eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
	609	ax.errorbar(self.data.phe, self.y[self.data.cut:], fmt='k^', xerr=self.data.ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
	610	plt.legend(loc='lower right')
	611	ax.yaxis.set_minor_locator(MultipleLocator(15))
	612	#plt.tight_layout()
	613
	614
	615
	616	class EDensityPlot(Plot):
	617	'''
	618	Plot for electron density
	619	'''
	620
	621	CODE = 'den'
	622	plot_name = 'Electron Density'
	623	plot_type = 'scatterbuffer'
	624
	625
	626	def setup(self):
	627
	628	self.ncols = 1
	629	self.nrows = 1
	630	self.nplots = 1
	631	self.ylabel = 'Range [km]'
	632	self.xlabel = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
	633	self.width = 4
	634	self.height = 6.5
	635	self.colorbar = False
	636	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
	637	if not self.titles:
	638	self.titles = self.data.parameters \
	639	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	640
	641	def plot(self):
	642
	643
	644	self.x = self.data[self.CODE]
	645	self.y = self.data.heights
	646	self.xmin = 1000
	647	self.xmax = 10000000
	648	ax = self.axes[0]
	649
	650	if ax.firsttime:
	651	self.autoxticks=False
	652	#if self.CODE=='den':
	653	ax.errorbar(self.data.dphi, self.y[:self.data.NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
	654	#ax.errorbar(self.data.dphi, self.y[:self.data.NSHTS], xerr=self.data.sdn1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
	655
	656	ax.errorbar(self.x[:,-1], self.y[:self.data.NSHTS], fmt='k^-', xerr=self.data.sdp2,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
	657	#else:
	658	#ax.errorbar(self.data.dphi[:self.data.cut], self.y[:self.data.cut], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
	659	#ax.errorbar(self.x[:self.data.cut,-1], self.y[:self.data.cut], fmt='k^-', xerr=self.data.sdp2[:self.data.cut],elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
	660
	661	if self.CODE=='denLP':
	662	ax.errorbar(self.data.ne[self.data.cut:], self.y[self.data.cut:], xerr=self.data.ene[self.data.cut:], fmt='r^-',elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
	663
	664	plt.legend(loc='upper right')
	665	ax.set_xscale("log", nonposx='clip')
	666	grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
	667	self.ystep_given=100
	668	if self.CODE=='denLP':
	669	self.ystep_given=200
	670	ax.set_yticks(grid_y_ticks,minor=True)
	671	ax.grid(which='minor')
	672	#plt.tight_layout()
	673
	674
	675
	676	else:
	677
	678	self.clear_figures()
	679	#if self.CODE=='den':
	680	ax.errorbar(self.data.dphi, self.y[:self.data.NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
	681	#ax.errorbar(self.data.dphi, self.y[:self.data.NSHTS], xerr=self.data.sdn1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
	682
	683	ax.errorbar(self.x[:,-1], self.y[:self.data.NSHTS], fmt='k^-', xerr=self.data.sdp2,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
	684	ax.errorbar(self.x[:,-2], self.y[:self.data.NSHTS], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
	685	#else:
	686	#ax.errorbar(self.data.dphi[:self.data.cut], self.y[:self.data.cut], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
	687	#ax.errorbar(self.x[:self.data.cut,-1], self.y[:self.data.cut], fmt='k^-', xerr=self.data.sdp2[:self.data.cut],elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
	688	#ax.errorbar(self.x[:self.data.cut,-2], self.y[:self.data.cut], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
	689
	690	if self.CODE=='denLP':
	691	ax.errorbar(self.data.ne[self.data.cut:], self.y[self.data.cut:], fmt='r^-', xerr=self.data.ene[self.data.cut:],elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
	692
	693	ax.set_xscale("log", nonposx='clip')
	694	grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
	695	ax.set_yticks(grid_y_ticks,minor=True)
	696	ax.grid(which='minor')
	697	plt.legend(loc='upper right')
	698	#plt.tight_layout()
	699
	700	class FaradayAnglePlot(Plot):
	701	'''
	702	Plot for electron density
	703	'''
	704
	705	CODE = 'FaradayAngle'
	706	plot_name = 'Faraday Angle'
	707	plot_type = 'scatterbuffer'
	708
	709
	710	def setup(self):
	711
	712	self.ncols = 1
	713	self.nrows = 1
	714	self.nplots = 1
	715	self.ylabel = 'Range [km]'
	716	self.xlabel = 'Faraday Angle (º)'
	717	self.width = 4
	718	self.height = 6.5
	719	self.colorbar = False
	720	if not self.titles:
	721	self.titles = self.data.parameters \
	722	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	723
	724	def plot(self):
	725
	726
	727	self.x = self.data[self.CODE]
	728	self.y = self.data.heights
	729	self.xmin = -180
	730	self.xmax = 180
	731	ax = self.axes[0]
	732
	733	if ax.firsttime:
	734	self.autoxticks=False
	735	#if self.CODE=='den':
	736	ax.plot(self.x, self.y,marker='o',color='g',linewidth=1.0,markersize=2)
	737
	738	grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
	739	self.ystep_given=100
	740	if self.CODE=='denLP':
	741	self.ystep_given=200
	742	ax.set_yticks(grid_y_ticks,minor=True)
	743	ax.grid(which='minor')
	744	#plt.tight_layout()
	745	else:
	746
	747	self.clear_figures()
	748	#if self.CODE=='den':
	749	#print(numpy.shape(self.x))
	750	ax.plot(self.x[:,-1], self.y, marker='o',color='g',linewidth=1.0, markersize=2)
	751
	752	grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
	753	ax.set_yticks(grid_y_ticks,minor=True)
	754	ax.grid(which='minor')
	755
	756	class EDensityHPPlot(EDensityPlot):
	757
	758	'''
	759	Plot for Electron Density Hybrid Experiment
	760	'''
	761
	762	CODE = 'denLP'
	763	plot_name = 'Electron Density'
	764	plot_type = 'scatterbuffer'
	765
	766
	767	class ACFsPlot(Plot):
	768	'''
	769	Plot for ACFs Double Pulse Experiment
	770	'''
	771
	772	CODE = 'acfs'
	773	plot_name = 'ACF'
	774	plot_type = 'scatterbuffer'
	775
	776
	777	def setup(self):
	778	#self.xaxis = 'time'
	779	self.ncols = 1
	780	self.nrows = 1
	781	self.nplots = 1
	782	self.ylabel = 'Range [km]'
	783	self.xlabel = 'lags (ms)'
	784	self.width = 3.5
	785	self.height = 6
	786	self.colorbar = False
	787	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
	788	if not self.titles:
	789	self.titles = self.data.parameters \
	790	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	791
	792	def plot(self):
	793
	794	self.x = self.data.lags_to_plot
	795	self.y = self.data['acfs'][:,-1]
	796
	797
	798	self.xmin = 0.0
	799	self.xmax = 2.0
	800
	801	ax = self.axes[0]
	802
	803	if ax.firsttime:
	804
	805	for i in range(self.data.NSHTS):
	806	x_aux = numpy.isfinite(self.x[i,:])
	807	y_aux = numpy.isfinite(self.y[i,:])
	808	yerr_aux = numpy.isfinite(self.data.acfs_error_to_plot[i,:])
	809	x_igcej_aux = numpy.isfinite(self.data.x_igcej_to_plot[i,:])
	810	y_igcej_aux = numpy.isfinite(self.data.y_igcej_to_plot[i,:])
	811	x_ibad_aux = numpy.isfinite(self.data.x_ibad_to_plot[i,:])
	812	y_ibad_aux = numpy.isfinite(self.data.y_ibad_to_plot[i,:])
	813	if self.x[i,:][~numpy.isnan(self.x[i,:])].shape[0]>2:
	814	ax.errorbar(self.x[i,x_aux], self.y[i,y_aux], yerr=self.data.acfs_error_to_plot[i,x_aux],color='b',marker='o',linewidth=1.0,markersize=2)
	815	ax.plot(self.data.x_igcej_to_plot[i,x_igcej_aux],self.data.y_igcej_to_plot[i,y_igcej_aux],'x',color='red',markersize=2)
	816	ax.plot(self.data.x_ibad_to_plot[i,x_ibad_aux],self.data.y_ibad_to_plot[i,y_ibad_aux],'X',color='red',markersize=2)
	817
	818	self.xstep_given = (self.xmax-self.xmin)/(self.data.DPL-1)
	819	self.ystep_given = 50
	820	ax.yaxis.set_minor_locator(MultipleLocator(15))
	821	ax.grid(which='minor')
	822
	823
	824
	825	else:
	826	self.clear_figures()
	827
	828	for i in range(self.data.NSHTS):
	829	x_aux = numpy.isfinite(self.x[i,:])
	830	y_aux = numpy.isfinite(self.y[i,:])
	831	yerr_aux = numpy.isfinite(self.data.acfs_error_to_plot[i,:])
	832	x_igcej_aux = numpy.isfinite(self.data.x_igcej_to_plot[i,:])
	833	y_igcej_aux = numpy.isfinite(self.data.y_igcej_to_plot[i,:])
	834	x_ibad_aux = numpy.isfinite(self.data.x_ibad_to_plot[i,:])
	835	y_ibad_aux = numpy.isfinite(self.data.y_ibad_to_plot[i,:])
	836	if self.x[i,:][~numpy.isnan(self.x[i,:])].shape[0]>2:
	837	ax.errorbar(self.x[i,x_aux], self.y[i,y_aux], yerr=self.data.acfs_error_to_plot[i,x_aux],linewidth=1.0,markersize=2,color='b',marker='o')
	838	ax.plot(self.data.x_igcej_to_plot[i,x_igcej_aux],self.data.y_igcej_to_plot[i,y_igcej_aux],'x',color='red',markersize=2)
	839	ax.plot(self.data.x_ibad_to_plot[i,x_ibad_aux],self.data.y_ibad_to_plot[i,y_ibad_aux],'X',color='red',markersize=2)
	840	ax.yaxis.set_minor_locator(MultipleLocator(15))
	841
	842
	843
	844
	845	class ACFsLPPlot(Plot):
	846	'''
	847	Plot for ACFs Double Pulse Experiment
	848	'''
	849
	850	CODE = 'acfs_LP'
	851	plot_name = 'ACF'
	852	plot_type = 'scatterbuffer'
	853
	854
	855	def setup(self):
	856	#self.xaxis = 'time'
	857	self.ncols = 1
	858	self.nrows = 1
	859	self.nplots = 1
	860	self.ylabel = 'Range [km]'
	861	self.xlabel = 'lags (ms)'
	862	self.width = 3.5
	863	self.height = 7
	864	self.colorbar = False
	865	if not self.titles:
	866	self.titles = self.data.parameters \
	867	if self.data.parameters else ['{}'.format(self.CODE.upper())]
	868
	869
	870
	871	def plot(self):
	872
	873	self.x = self.data.lags_LP_to_plot
	874	self.y = self.data['acfs_LP'][:,-1]
	875
	876	self.xmin = 0.0
	877	self.xmax = 1.5
	878
	879	ax = self.axes[0]
	880
	881	if ax.firsttime:
	882
	883	for i in range(self.data.NACF):
	884	x_aux = numpy.isfinite(self.x[i,:])
	885	y_aux = numpy.isfinite(self.y[i,:])
	886	yerr_aux = numpy.isfinite(self.data.errors[i,:])
	887
	888	if self.x[i,:][~numpy.isnan(self.x[i,:])].shape[0]>2:
	889	ax.errorbar(self.x[i,x_aux], self.y[i,y_aux], yerr=self.data.errors[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
	890
	891	#self.xstep_given = (self.xmax-self.xmin)/(self.data.NLAG-1)
	892	self.xstep_given=0.3
	893	self.ystep_given = 200
	894	ax.yaxis.set_minor_locator(MultipleLocator(15))
	895	ax.grid(which='minor')
	896
	897	else:
	898	self.clear_figures()
	899
	900	for i in range(self.data.NACF):
	901	x_aux = numpy.isfinite(self.x[i,:])
	902	y_aux = numpy.isfinite(self.y[i,:])
	903	yerr_aux = numpy.isfinite(self.data.errors[i,:])
	904
	905	if self.x[i,:][~numpy.isnan(self.x[i,:])].shape[0]>2:
	906	ax.errorbar(self.x[i,x_aux], self.y[i,y_aux], yerr=self.data.errors[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
	907
	908	ax.yaxis.set_minor_locator(MultipleLocator(15))
	909
	910
	911	class CrossProductsPlot(Plot):
	912	'''
	913	Plot for cross products
	914	'''
	915
	916	CODE = 'crossprod'
	917	plot_name = 'Cross Products'
	918	plot_type = 'scatterbuffer'
	919
	920
	921	def setup(self):
	922
	923	self.ncols = 3
	924	self.nrows = 1
	925	self.nplots = 3
	926	self.ylabel = 'Range [km]'
	927
	928	self.width = 3.5*self.nplots
	929	self.height = 5.5
	930	self.colorbar = False
	931	self.titles = []
	932
	933	def plot(self):
	934
	935	self.x = self.data['crossprod'][:,-1,:,:,:,:]
	936
	937
	938
	939
	940	self.y = self.data.heights[0:self.data.NDP]
	941
	942
	943
	944	for n, ax in enumerate(self.axes):
	945
	946	self.xmin=numpy.min(numpy.concatenate((self.x[n][0,20:30,0,0],self.x[n][1,20:30,0,0],self.x[n][2,20:30,0,0],self.x[n][3,20:30,0,0])))
	947	self.xmax=numpy.max(numpy.concatenate((self.x[n][0,20:30,0,0],self.x[n][1,20:30,0,0],self.x[n][2,20:30,0,0],self.x[n][3,20:30,0,0])))
	948
	949
	950	if ax.firsttime:
	951
	952	self.autoxticks=False
	953	if n==0:
	954	label1='kax'
	955	label2='kay'
	956	label3='kbx'
	957	label4='kby'
	958	self.xlimits=[(self.xmin,self.xmax)]
	959	elif n==1:
	960	label1='kax2'
	961	label2='kay2'
	962	label3='kbx2'
	963	label4='kby2'
	964	self.xlimits.append((self.xmin,self.xmax))
	965	elif n==2:
	966	label1='kaxay'
	967	label2='kbxby'
	968	label3='kaxbx'
	969	label4='kaxby'
	970	self.xlimits.append((self.xmin,self.xmax))
	971
	972
	973	ax.plotline1 = ax.plot(self.x[n][0,:,0,0], self.y, color='r',linewidth=2.0, label=label1)
	974	ax.plotline2 = ax.plot(self.x[n][1,:,0,0], self.y, color='k',linewidth=2.0, label=label2)
	975	ax.plotline3 = ax.plot(self.x[n][2,:,0,0], self.y, color='b',linewidth=2.0, label=label3)
	976	ax.plotline4 = ax.plot(self.x[n][3,:,0,0], self.y, color='m',linewidth=2.0, label=label4)
	977	ax.legend(loc='upper right')
	978	ax.set_xlim(self.xmin, self.xmax)
	979	self.titles.append('{}'.format(self.plot_name.upper()))
	980	#plt.tight_layout()
	981
	982
	983	else:
	984
	985	if n==0:
	986	self.xlimits=[(self.xmin,self.xmax)]
	987	else:
	988	self.xlimits.append((self.xmin,self.xmax))
	989
	990	ax.set_xlim(self.xmin, self.xmax)
	991
	992
	993	ax.plotline1[0].set_data(self.x[n][0,:,0,0],self.y)
	994	ax.plotline2[0].set_data(self.x[n][1,:,0,0],self.y)
	995	ax.plotline3[0].set_data(self.x[n][2,:,0,0],self.y)
	996	ax.plotline4[0].set_data(self.x[n][3,:,0,0],self.y)
	997	self.titles.append('{}'.format(self.plot_name.upper()))
	998	#plt.tight_layout()
	999
	1000
	1001
	1002	class CrossProductsLPPlot(Plot):
	1003	'''
	1004	Plot for cross products LP
	1005	'''
	1006
	1007	CODE = 'crossprodlp'
	1008	plot_name = 'Cross Products LP'
	1009	plot_type = 'scatterbuffer'
	1010
	1011
	1012	def setup(self):
	1013
	1014	self.ncols = 2
	1015	self.nrows = 1
	1016	self.nplots = 2
	1017	self.ylabel = 'Range [km]'
	1018	self.xlabel = 'dB'
	1019	self.width = 3.5*self.nplots
	1020	self.height = 5.5
	1021	self.colorbar = False
	1022	self.titles = []
	1023	self.plotline_array=numpy.zeros((2,self.data.NLAG),dtype=object)
	1024	def plot(self):
	1025
	1026
	1027	self.x = self.data[self.CODE][:,-1,:,:]
	1028
	1029
	1030	self.y = self.data.heights[0:self.data.NRANGE]
	1031
	1032
	1033	label_array=numpy.array(['lag '+ str(x) for x in range(self.data.NLAG)])
	1034	color_array=['r','k','g','b','c','m','y','orange','steelblue','purple','peru','darksalmon','grey','limegreen','olive','midnightblue']
	1035
	1036
	1037	for n, ax in enumerate(self.axes):
	1038
	1039	self.xmin=30
	1040	self.xmax=70
	1041	#print(self.x[0,12:15,n])
	1042	#input()
	1043	#self.xmin=numpy.min(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
	1044	#self.xmax=numpy.max(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
	1045
	1046	#print("before",self.plotline_array)
	1047
	1048	if ax.firsttime:
	1049
	1050	self.autoxticks=False
	1051
	1052
	1053	for i in range(self.data.NLAG):
	1054	#print(i)
	1055	#print(numpy.shape(self.x))
	1056	self.plotline_array[n,i], = ax.plot(self.x[i,:,n], self.y, color=color_array[i],linewidth=1.0, label=label_array[i])
	1057	#ax.plotline1 = ax.plot(self.x[0,:,n], self.y, color='r',linewidth=2.0, label=label_array[0])
	1058	#ax.plotline2 = ax.plot(self.x[n][1,:,0,0], self.y, color='k',linewidth=2.0, label=label2)
	1059	#ax.plotline3 = ax.plot(self.x[n][2,:,0,0], self.y, color='b',linewidth=2.0, label=label3)
	1060	#ax.plotline4 = ax.plot(self.x[n][3,:,0,0], self.y, color='m',linewidth=2.0, label=label4)
	1061
	1062
	1063	#print(self.plotline_array)
	1064
	1065
	1066
	1067	ax.legend(loc='upper right')
	1068	ax.set_xlim(self.xmin, self.xmax)
	1069	if n==0:
	1070	self.titles.append('{} CH0'.format(self.plot_name.upper()))
	1071	if n==1:
	1072	self.titles.append('{} CH1'.format(self.plot_name.upper()))
	1073
	1074	#plt.tight_layout()
	1075
	1076	else:
	1077	#print(self.plotline_array)
	1078	for i in range(self.data.NLAG):
	1079
	1080	self.plotline_array[n,i].set_data(self.x[i,:,n],self.y)
	1081
	1082
	1083
	1084	#ax.plotline1[0].set_data(self.x[n][0,:,0,0],self.y)
	1085	#ax.plotline2[0].set_data(self.x[n][1,:,0,0],self.y)
	1086	#ax.plotline3[0].set_data(self.x[n][2,:,0,0],self.y)
	1087	#ax.plotline4[0].set_data(self.x[n][3,:,0,0],self.y)
	1088
	1089	if n==0:
	1090	self.titles.append('{} CH0'.format(self.plot_name.upper()))
	1091	if n==1:
	1092	self.titles.append('{} CH1'.format(self.plot_name.upper()))
	1093
	1094	#plt.tight_layout()
	1095
	1096
	1097	class NoiseDPPlot(NoisePlot):
	1098	'''
	1099	Plot for noise Double Pulse
	1100	'''
	1101
	1102	CODE = 'noisedp'
	1103	plot_name = 'Noise'
	1104	plot_type = 'scatterbuffer'
	1105
	1106
	1107	class XmitWaveformPlot(Plot):
	1108	'''
	1109	Plot for xmit waveform
	1110	'''
	1111
	1112	CODE = 'xmit'
	1113	plot_name = 'Xmit Waveform'
	1114	plot_type = 'scatterbuffer'
	1115
	1116
	1117	def setup(self):
	1118
	1119	self.ncols = 1
	1120	self.nrows = 1
	1121	self.nplots = 1
	1122	self.ylabel = ''
	1123	self.xlabel = 'Number of Lag'
	1124	self.width = 5.5
	1125	self.height = 3.5
	1126	self.colorbar = False
	1127	if not self.titles:
	1128	self.titles = self.data.parameters \
	1129	if self.data.parameters else ['{}'.format(self.plot_name.upper())]
	1130
	1131	def plot(self):
	1132
	1133	self.x = numpy.arange(0,self.data.NLAG,1,'float32')
	1134	self.y = self.data['xmit'][:,-1,:]
	1135
	1136	self.xmin = 0
	1137	self.xmax = self.data.NLAG-1
	1138	self.ymin = -1.0
	1139	self.ymax = 1.0
	1140	ax = self.axes[0]
	1141
	1142	if ax.firsttime:
	1143	ax.plotline0=ax.plot(self.x,self.y[0,:],color='blue')
	1144	ax.plotline1=ax.plot(self.x,self.y[1,:],color='red')
	1145	secax=ax.secondary_xaxis(location=0.5)
	1146	secax.xaxis.tick_bottom()
	1147	secax.tick_params( labelleft=False, labeltop=False,
	1148	labelright=False, labelbottom=False)
	1149
	1150	self.xstep_given = 3
	1151	self.ystep_given = .25
	1152	secax.set_xticks(numpy.linspace(self.xmin, self.xmax, 6)) #only works on matplotlib.version>3.2
	1153
	1154	else:
	1155	ax.plotline0[0].set_data(self.x,self.y[0,:])
	1156	ax.plotline1[0].set_data(self.x,self.y[1,:])

schainpy/model/io/jroIO_dat.py created 644 +683 0

This diff has been collapsed as it changes many lines, (683 lines changed) Show them Hide them
		@@ -0,0 +1,683
	1	'''
	2	Created on Jun 9, 2020
	3
	4	@author: Roberto Flores
	5	'''
	6
	7	import os
	8	import sys
	9	import time
	10
	11	import struct
	12
	13
	14	import datetime
	15
	16	import numpy
	17
	18
	19	import schainpy.admin
	20	from schainpy.model.io.jroIO_base import LOCALTIME, Reader
	21	from schainpy.model.data.jroheaderIO import BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
	22	from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
	23	from schainpy.model.data.jrodata import Voltage, Parameters
	24	from schainpy.utils import log
	25
	26
	27	class DatReader(Reader, ProcessingUnit):
	28
	29	def __init__(self):
	30
	31	ProcessingUnit.__init__(self)
	32	self.basicHeaderObj = BasicHeader(LOCALTIME)
	33	self.systemHeaderObj = SystemHeader()
	34	self.radarControllerHeaderObj = RadarControllerHeader()
	35	self.processingHeaderObj = ProcessingHeader()
	36	self.dataOut = Parameters()
	37	#print(self.basicHeaderObj.timezone)
	38	#self.counter_block=0
	39	self.format='dat'
	40	self.flagNoMoreFiles = 0
	41	self.filename = None
	42	self.intervals = set()
	43	#self.datatime = datetime.datetime(1900,1,1)
	44
	45	self.filefmt = "*%Y%m%d*****"
	46
	47	self.padding=numpy.zeros(1,'int32')
	48	self.hsize=numpy.zeros(1,'int32')
	49	self.bufsize=numpy.zeros(1,'int32')
	50	self.nr=numpy.zeros(1,'int32')
	51	self.ngates=numpy.zeros(1,'int32') ### ### ### 2
	52	self.time1=numpy.zeros(1,'uint64') # pos 3
	53	self.time2=numpy.zeros(1,'uint64') # pos 4
	54	self.lcounter=numpy.zeros(1,'int32')
	55	self.groups=numpy.zeros(1,'int32')
	56	self.system=numpy.zeros(4,'int8') # pos 7
	57	self.h0=numpy.zeros(1,'float32')
	58	self.dh=numpy.zeros(1,'float32')
	59	self.ipp=numpy.zeros(1,'float32')
	60	self.process=numpy.zeros(1,'int32')
	61	self.tx=numpy.zeros(1,'int32')
	62
	63	self.ngates1=numpy.zeros(1,'int32') ### ### ### 13
	64	self.time0=numpy.zeros(1,'uint64') # pos 14
	65	self.nlags=numpy.zeros(1,'int32')
	66	self.nlags1=numpy.zeros(1,'int32')
	67	self.txb=numpy.zeros(1,'float32') ### ### ### 17
	68	self.time3=numpy.zeros(1,'uint64') # pos 18
	69	self.time4=numpy.zeros(1,'uint64') # pos 19
	70	self.h0_=numpy.zeros(1,'float32')
	71	self.dh_=numpy.zeros(1,'float32')
	72	self.ipp_=numpy.zeros(1,'float32')
	73	self.txa_=numpy.zeros(1,'float32')
	74
	75	self.pad=numpy.zeros(100,'int32')
	76
	77	self.nbytes=numpy.zeros(1,'int32')
	78	self.limits=numpy.zeros(1,'int32')
	79	self.ngroups=numpy.zeros(1,'int32') ### ### ### 27
	80	#Make the header list
	81	#header=[hsize,bufsize,nr,ngates,time1,time2,lcounter,groups,system,h0,dh,ipp,process,tx,padding,ngates1,time0,nlags,nlags1,padding,txb,time3,time4,h0_,dh_,ipp_,txa_,pad,nbytes,limits,padding,ngroups]
	82	self.header=[self.hsize,self.bufsize,self.nr,self.ngates,self.time1,self.time2,self.lcounter,self.groups,self.system,self.h0,self.dh,self.ipp,self.process,self.tx,self.ngates1,self.padding,self.time0,self.nlags,self.nlags1,self.padding,self.txb,self.time3,self.time4,self.h0_,self.dh_,self.ipp_,self.txa_,self.pad,self.nbytes,self.limits,self.padding,self.ngroups]
	83
	84
	85
	86	def setup(self, **kwargs):
	87
	88	self.set_kwargs(**kwargs)
	89
	90
	91	if self.path is None:
	92	raise ValueError('The path is not valid')
	93
	94	self.open_file = open
	95	self.open_mode = 'rb'
	96
	97
	98
	99	if self.format is None:
	100	raise ValueError('The format is not valid')
	101	elif self.format.lower() in ('dat'):
	102	self.ext = '.dat'
	103	elif self.format.lower() in ('out'):
	104	self.ext = '.out'
	105
	106
	107	log.log("Searching files in {}".format(self.path), self.name)
	108	self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
	109	self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
	110	#print(self.path)
	111	#print(self.filenameList)
	112	#input()
	113
	114
	115	self.setNextFile()
	116
	117	def readFirstHeader(self):
	118	'''Read header and data'''
	119
	120	#self.flag_same_file=1
	121	self.counter_block=0
	122	self.parseHeader()
	123	self.parseData()
	124	self.blockIndex = 0
	125
	126	return
	127
	128	def parseHeader(self):
	129	'''
	130	'''
	131
	132	for i in range(len(self.header)):
	133	for j in range(len(self.header[i])):
	134	#print("len(header[i]) ",len(header[i]))
	135	#input()
	136	temp=self.fp.read(int(self.header[i].itemsize))
	137	if isinstance(self.header[i][0], numpy.int32):
	138	#print(struct.unpack('i', temp)[0])
	139	self.header[i][0]=struct.unpack('i', temp)[0]
	140	if isinstance(self.header[i][0], numpy.uint64):
	141	self.header[i][0]=struct.unpack('q', temp)[0]
	142	if isinstance(self.header[i][0], numpy.int8):
	143	self.header[i][0]=struct.unpack('B', temp)[0]
	144	if isinstance(self.header[i][0], numpy.float32):
	145	self.header[i][0]=struct.unpack('f', temp)[0]
	146
	147	self.fp.seek(0,0)
	148	if int(self.header[1][0])==int(81864):
	149	self.experiment='DP'
	150
	151	elif int(self.header[1][0])==int(185504):
	152	self.experiment='HP'
	153
	154
	155	self.total_blocks=os.stat(self.filename).st_size//self.header[1][0]
	156
	157
	158	def parseData(self):
	159	'''
	160	'''
	161	if self.experiment=='DP':
	162	self.header[15][0]=66
	163	self.header[18][0]=16
	164	self.header[17][0]=11
	165	self.header[2][0]=2
	166
	167
	168	self.noise=numpy.zeros(self.header[2][0],'float32') #self.header[2][0]
	169	#tmpx=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	170	self.kax=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	171	self.kay=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	172	self.kbx=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	173	self.kby=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	174	self.kax2=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	175	self.kay2=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	176	self.kbx2=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	177	self.kby2=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	178	self.kaxbx=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	179	self.kaxby=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	180	self.kaybx=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	181	self.kayby=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	182	self.kaxay=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	183	self.kbxby=numpy.zeros((self.header[15][0],self.header[17][0],2),'float32')
	184	self.output_LP_real=numpy.zeros((self.header[18][0],200,self.header[2][0]),'float32')
	185	self.output_LP_imag=numpy.zeros((self.header[18][0],200,self.header[2][0]),'float32')
	186	self.final_cross_products=[self.kax,self.kay,self.kbx,self.kby,self.kax2,self.kay2,self.kbx2,self.kby2,self.kaxbx,self.kaxby,self.kaybx,self.kayby,self.kaxay,self.kbxby]
	187	#self.final_cross_products=[tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx,tmpx]
	188
	189	#print("pos: ",self.fp.tell())
	190
	191
	192	def readNextBlock(self):
	193
	194	while True:
	195	self.flagDiscontinuousBlock = 0
	196	#print(os.stat(self.filename).st_size)
	197	#print(os.stat(self.filename).st_size//self.header[1][0])
	198	#os.stat(self.fp)
	199	if self.counter_block == self.total_blocks:
	200
	201	self.setNextFile()
	202
	203	self.readBlock()
	204	#self.counter_block+=1
	205
	206	if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
	207	(self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
	208
	209	#print(self.datatime)
	210	#print(datetime.datetime.combine(self.startDate, self.startTime))
	211	#print(datetime.datetime.combine(self.endDate, self.endTime))
	212	#print("warning")
	213	log.warning(
	214	'Reading Block No. {}/{} -> {} [Skipping]'.format(
	215	self.counter_block,
	216	self.total_blocks,
	217	self.datatime.ctime()),
	218	'DATReader')
	219	continue
	220	break
	221
	222	log.log(
	223	'Reading Block No. {}/{} -> {}'.format(
	224	self.counter_block,
	225	self.total_blocks,
	226	self.datatime.ctime()),
	227	'DATReader')
	228
	229	return 1
	230
	231	def readBlock(self):
	232	'''
	233	'''
	234
	235	self.npos=self.counter_block*self.header[1][0]
	236	#print(self.counter_block)
	237	self.fp.seek(self.npos, 0)
	238	self.counter_block+=1
	239	#print("fpos1: ",self.fp.tell())
	240
	241	self.read_header()
	242
	243	#put by hand because old files didn't save it in the header
	244	if self.experiment=='DP':
	245	self.header[15][0]=66
	246	self.header[18][0]=16
	247	self.header[17][0]=11
	248	self.header[2][0]=2
	249	#########################################
	250
	251	if self.experiment=="HP":
	252	self.long_pulse_products()
	253
	254	self.read_cross_products()
	255
	256
	257	self.read_noise()
	258
	259
	260	return
	261
	262
	263
	264	def read_header(self):
	265
	266
	267	for i in range(len(self.header)):
	268	for j in range(len(self.header[i])):
	269	#print("len(header[i]) ",len(header[i]))
	270	#input()
	271	temp=self.fp.read(int(self.header[i].itemsize))
	272	#if(b''==temp):
	273	# self.setNextFile()
	274	# self.flag_same_file=0
	275	if isinstance(self.header[i][0], numpy.int32):
	276	#print(struct.unpack('i', temp)[0])
	277	self.header[i][0]=struct.unpack('i', temp)[0]
	278	if isinstance(self.header[i][0], numpy.uint64):
	279	self.header[i][0]=struct.unpack('q', temp)[0]
	280	if isinstance(self.header[i][0], numpy.int8):
	281	self.header[i][0]=struct.unpack('B', temp)[0]
	282	if isinstance(self.header[i][0], numpy.float32):
	283	self.header[i][0]=struct.unpack('f', temp)[0]
	284	#else:
	285	# continue
	286	#self.fp.seek(self.npos_aux, 0)
	287	# break
	288
	289	#print("fpos2: ",self.fp.tell())
	290	#log.success('Parameters found: {}'.format(self.parameters),
	291	# 'DATReader')
	292	#print("Success")
	293	#self.TimeBlockSeconds_for_dp_power = self.header[4][0]#-((self.dataOut.nint-1)self.dataOut.NAVG2)
	294	#print(dataOut.TimeBlockSeconds_for_dp_power)
	295
	296	#self.datatime=datetime.datetime.fromtimestamp(self.header[4][0]).strftime("%Y-%m-%d %H:%M:%S")
	297	#print(self.header[4][0])
	298	self.datatime=datetime.datetime.fromtimestamp(self.header[4][0])
	299	#print(self.header[1][0])
	300
	301	def long_pulse_products(self):
	302	temp=self.fp.read(self.header[18][0]self.header[2][0]200*8)
	303	ii=0
	304
	305	for l in range(self.header[18][0]): #lag
	306	for r in range(self.header[2][0]): # channels
	307	for k in range(200): #RANGE## generalizar
	308	self.output_LP_real[l,k,r]=struct.unpack('f', temp[ii:ii+4])[0]
	309	ii=ii+4
	310	self.output_LP_imag[l,k,r]=struct.unpack('f', temp[ii:ii+4])[0]
	311	ii=ii+4
	312
	313	#print(self.output_LP_real[1,1,1])
	314	#print(self.output_LP_imag[1,1,1])
	315	def read_cross_products(self):
	316
	317	for ind in range(len(self.final_cross_products)): #final cross products
	318	temp=self.fp.read(self.header[17][0]2self.header[15][0]4) #4 bytes
	319	#if(b''==temp):
	320	# self.setNextFile()
	321	# self.flag_same_file=0
	322	ii=0
	323	#print("kabxys.shape ",kabxys.shape)
	324	#print(kabxys)
	325	#print("fpos3: ",self.fp.tell())
	326	for l in range(self.header[17][0]): #lag
	327	#print("fpos3: ",self.fp.tell())
	328	for fl in range(2): # unflip and flip
	329	for k in range(self.header[15][0]): #RANGE
	330	#print("fpos3: ",self.fp.tell())
	331	self.final_cross_products[ind][k,l,fl]=struct.unpack('f', temp[ii:ii+4])[0]
	332	ii=ii+4
	333	#print("fpos2: ",self.fp.tell())
	334
	335
	336
	337	def read_noise(self):
	338
	339	temp=self.fp.read(self.header[2][0]4) #4 bytes self.header[2][0]
	340	for ii in range(self.header[2][0]): #self.header[2][0]
	341	self.noise[ii]=struct.unpack('f', temp[ii4:(ii+1)4])[0]
	342
	343	#print("fpos5: ",self.fp.tell())
	344
	345
	346
	347	def set_output(self):
	348	'''
	349	Storing data from buffer to dataOut object
	350	'''
	351	#print("fpos2: ",self.fp.tell())
	352	##self.dataOut.header = self.header
	353	#this is put by hand because it isn't saved in the header
	354	if self.experiment=='DP':
	355	self.dataOut.NRANGE=0
	356	self.dataOut.NSCAN=132
	357	self.dataOut.heightList=self.header[10][0]*(numpy.arange(self.header[15][0]))
	358	elif self.experiment=='HP':
	359	self.dataOut.output_LP=self.output_LP_real+1.j*self.output_LP_imag
	360	self.dataOut.NRANGE=200
	361	self.dataOut.NSCAN=128
	362	self.dataOut.heightList=self.header[10][0]*(numpy.arange(90)) #NEEEDS TO BE GENERALIZED
	363	#########################################
	364	#print(self.dataOut.output_LP[1,1,1])
	365	self.dataOut.MAXNRANGENDT=self.header[3][0]
	366	self.dataOut.NDP=self.header[15][0]
	367	self.dataOut.DPL=self.header[17][0]
	368	self.dataOut.DH=self.header[10][0]
	369	self.dataOut.NAVG=self.header[7][0]
	370	self.dataOut.H0=self.header[9][0]
	371	self.dataOut.NR=self.header[2][0]
	372	self.dataOut.NLAG=self.header[18][0]
	373	#self.dataOut.tmpx=self.tmpx
	374	#self.dataOut.timeZone = 5
	375	#self.dataOut.final_cross_products=self.final_cross_products
	376	self.dataOut.kax=self.kax
	377	#print(self.dataOut.kax[1,1,1])
	378	self.dataOut.kay=self.kay
	379	self.dataOut.kbx=self.kbx
	380	self.dataOut.kby=self.kby
	381	self.dataOut.kax2=self.kax2
	382	self.dataOut.kay2=self.kay2
	383	self.dataOut.kbx2=self.kbx2
	384	self.dataOut.kby2=self.kby2
	385	self.dataOut.kaxbx=self.kaxbx
	386	self.dataOut.kaxby=self.kaxby
	387	self.dataOut.kaybx=self.kaybx
	388	self.dataOut.kayby=self.kayby
	389	self.dataOut.kaxay=self.kaxay
	390	self.dataOut.kbxby=self.kbxby
	391	self.dataOut.noise_final=self.noise
	392	#print("NOISE",self.noise)
	393
	394
	395	self.dataOut.useLocalTime=True
	396
	397	#self.dataOut.experiment=self.experiment
	398	#print(self.datatime)
	399	#print(self.dataOut.datatime)
	400
	401
	402	#self.dataOut.utctime = (self.datatime - datetime.datetime(1970, 1, 1)).total_seconds()
	403	#self.dataOut.utctimeInit = self.dataOut.utctime
	404
	405
	406
	407	self.dataOut.lt=self.datatime.hour
	408
	409
	410	#print(RadarControllerHeader().ippSeconds)
	411	#print(RadarControllerHeader().ipp)
	412	#self.dataOut.utctime=time.gmtime(self.header[4][0])- datetime.datetime(1970, 1, 1)
	413	#self.dataOut.utctime=self.dataOut.utctime.total_seconds()
	414	#time1 = self.header[4][0] # header.time1
	415	#print("time1: ",time1)
	416	#print(self.header[4][0])
	417	#date = time.ctime(time1)
	418	#print("DADSADA",time.strptime(date))
	419	#print("date_before: ",date)
	420	#bd_time=time.gmtime(time1)
	421	#print(time.mktime(bd_time))
	422	#self.dataOut.utctime=time.mktime(bd_time)
	423	self.dataOut.utctime = self.header[4][0]
	424	#self.dataOut.datatime=a
	425	#print(datetime.datetime.utcfromtimestamp(self.dataOut.utctime))
	426	#self.dataOut.TimeBlockDate=self.datatime.ctime()
	427	self.dataOut.TimeBlockSeconds=time.mktime(time.strptime(self.dataOut.datatime.ctime()))
	428
	429	#self.dataOut.heightList = self.ranges
	430	#self.dataOut.utctime = (self.datatime - datetime.datetime(1970, 1, 1)).total_seconds()
	431	#self.dataOut.utctimeInit = self.dataOut.utctime
	432	#self.dataOut.paramInterval = min(self.intervals)
	433	#self.dataOut.useLocalTime = False
	434	self.dataOut.flagNoData = False
	435	self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
	436	#print(self.dataOut.channelIndexList)
	437	self.dataOut.channelList=list(range(0,self.header[2][0]))
	438	#print(self.dataOut.channelList)
	439	#print(self.datatime)
	440	#print(self.dataOut.final_cross_products[0])
	441
	442
	443	#self.dataOut.heightList=self.header[10][0]*(numpy.arange(self.header[15][0]))
	444
	445	#print(numpy.shape(self.dataOut.heightList))
	446
	447
	448	def getData(self):
	449	'''
	450	Storing data from databuffer to dataOut object
	451	'''
	452
	453	if not self.readNextBlock():
	454	self.dataOut.flagNoData = True
	455	return 0
	456
	457	self.set_output()
	458
	459	return 1
	460
	461	def run(self, **kwargs):
	462
	463	if not(self.isConfig):
	464	self.setup(**kwargs)
	465	self.isConfig = True
	466	#print("fpos1: ",self.fp.tell())
	467	self.getData()
	468
	469	return
	470
	471	@MPDecorator
	472	class DatWriter(Operation):
	473
	474
	475	def __init__(self):
	476
	477	Operation.__init__(self)
	478	#self.dataOut = Voltage()
	479	self.counter = 0
	480	self.path = None
	481	self.fp = None
	482	return
	483	#self.ext= '.dat'
	484
	485	def run(self, dataOut, path, format='dat', experiment=None, **kwargs):
	486	print(dataOut.flagNoData)
	487	print(dataOut.datatime.ctime())
	488	print(dataOut.TimeBlockDate)
	489	input()
	490	#if dataOut.flag_save:
	491	self.experiment=experiment
	492	self.path=path
	493	if self.experiment=='DP':
	494	dataOut.header[1][0]=81864
	495	elif self.experiment=='HP':
	496	dataOut.header[1][0]=185504#173216
	497	#dataOut.header[1][0]=bufsize
	498	self.dataOut = dataOut
	499	#print(self.dataOut.nint)
	500	#self.bufsize=bufsize
	501	if format == 'dat':
	502	self.ext = '.dat'
	503	if format == 'out':
	504	self.ext = '.out'
	505	self.putData()
	506
	507	return
	508
	509
	510
	511	def setFile(self):
	512	'''
	513	Create new out file object
	514	'''
	515
	516	#self.dataOut.TimeBlockSeconds=time.mktime(time.strptime(self.dataOut.TimeBlockDate))
	517	date = datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds)
	518
	519	#print("date",date)
	520
	521	filename = '{}{}{}'.format('jro',
	522	date.strftime('%Y%m%d_%H%M%S'),
	523	self.ext)
	524	#print(filename)
	525	#print(self.path)
	526
	527	self.fullname = os.path.join(self.path, filename)
	528
	529	if os.path.isfile(self.fullname) :
	530	log.warning(
	531	'Destination file {} already exists, previous file deleted.'.format(
	532	self.fullname),
	533	'DatWriter')
	534	os.remove(self.fullname)
	535
	536	try:
	537	log.success(
	538	'Creating file: {}'.format(self.fullname),
	539	'DatWriter')
	540	if not os.path.exists(self.path):
	541	os.makedirs(self.path)
	542	#self.fp = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
	543	self.fp = open(self.fullname,'wb')
	544
	545	except ValueError as e:
	546	log.error(
	547	'Impossible to create *.out file',
	548	'DatWriter')
	549	return
	550
	551	return 1
	552
	553	def writeBlock(self):
	554
	555	#self.dataOut.paramInterval=2
	556	#startTime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
	557	#print(startTime)
	558	#endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
	559
	560	self.dataOut.header[0].astype('int32').tofile(self.fp)
	561	self.dataOut.header[1].astype('int32').tofile(self.fp)
	562	self.dataOut.header[2].astype('int32').tofile(self.fp)
	563	self.dataOut.header[3].astype('int32').tofile(self.fp)
	564	self.dataOut.header[4].astype('uint64').tofile(self.fp)
	565	self.dataOut.header[5].astype('uint64').tofile(self.fp)
	566	self.dataOut.header[6].astype('int32').tofile(self.fp)
	567	self.dataOut.header[7].astype('int32').tofile(self.fp)
	568	#print(dataOut.header[7])
	569	self.dataOut.header[8].astype('int8').tofile(self.fp)
	570	self.dataOut.header[9].astype('float32').tofile(self.fp)
	571	self.dataOut.header[10].astype('float32').tofile(self.fp)
	572	self.dataOut.header[11].astype('float32').tofile(self.fp)
	573	self.dataOut.header[12].astype('int32').tofile(self.fp)
	574	self.dataOut.header[13].astype('int32').tofile(self.fp)
	575	self.dataOut.header[14].astype('int32').tofile(self.fp)
	576	self.dataOut.header[15].astype('int32').tofile(self.fp)
	577	self.dataOut.header[16].astype('uint64').tofile(self.fp)
	578	self.dataOut.header[17].astype('int32').tofile(self.fp)
	579	self.dataOut.header[18].astype('int32').tofile(self.fp)
	580	self.dataOut.header[19].astype('int32').tofile(self.fp)
	581	self.dataOut.header[20].astype('float32').tofile(self.fp)
	582	self.dataOut.header[21].astype('uint64').tofile(self.fp)
	583	self.dataOut.header[22].astype('uint64').tofile(self.fp)
	584	self.dataOut.header[23].astype('float32').tofile(self.fp)
	585	self.dataOut.header[24].astype('float32').tofile(self.fp)
	586	self.dataOut.header[25].astype('float32').tofile(self.fp)
	587	self.dataOut.header[26].astype('float32').tofile(self.fp)
	588	self.dataOut.header[27].astype('int32').tofile(self.fp)
	589	self.dataOut.header[28].astype('int32').tofile(self.fp)
	590	self.dataOut.header[29].astype('int32').tofile(self.fp)
	591	self.dataOut.header[30].astype('int32').tofile(self.fp)
	592	self.dataOut.header[31].astype('int32').tofile(self.fp)
	593	#print("tell before 1 ",self.fp.tell())
	594	#input()
	595
	596	if self.experiment=="HP":
	597	#print("INSIDE")
	598	#tmp=numpy.zeros(1,dtype='complex64')
	599	#print("tmp ",tmp)
	600	#input()
	601	#print(dataOut.NLAG)
	602	#print(dataOut.NR)
	603	#print(dataOut.NRANGE)
	604	for l in range(self.dataOut.NLAG): #lag
	605	for r in range(self.dataOut.NR): # unflip and flip
	606	for k in range(self.dataOut.NRANGE): #RANGE
	607	self.dataOut.output_LP.real[l,k,r].astype('float32').tofile(self.fp)
	608	self.dataOut.output_LP.imag[l,k,r].astype('float32').tofile(self.fp)
	609
	610
	611	#print("tell before 2 ",self.outputfile.tell())
	612
	613
	614
	615
	616
	617	#print(self.dataOut.output_LP[1,1,1])
	618
	619	#print(self.dataOut.kax)
	620	final_cross_products=[self.dataOut.kax,self.dataOut.kay,self.dataOut.kbx,self.dataOut.kby,
	621	self.dataOut.kax2,self.dataOut.kay2,self.dataOut.kbx2,self.dataOut.kby2,
	622	self.dataOut.kaxbx,self.dataOut.kaxby,self.dataOut.kaybx,self.dataOut.kayby,
	623	self.dataOut.kaxay,self.dataOut.kbxby]
	624
	625	#print(self.dataOut.kax)
	626	#print("tell before crossp saving ",self.outputfile.tell())
	627	for kabxys in final_cross_products:
	628
	629	for l in range(self.dataOut.DPL): #lag
	630	for fl in range(2): # unflip and flip
	631	for k in range(self.dataOut.NDT): #RANGE
	632	kabxys[k,l,fl].astype('float32').tofile(self.fp)
	633
	634
	635	#print("tell before noise saving ",self.outputfile.tell())
	636
	637
	638	for nch in range(self.dataOut.NR):
	639	self.dataOut.noise_final[nch].astype('float32').tofile(self.fp)
	640
	641	#print("tell before noise saving ",self.fp.tell())
	642	#input()
	643
	644
	645
	646
	647	log.log(
	648	'Writing {} blocks'.format(
	649	self.counter+1),
	650	'DatWriter')
	651
	652
	653
	654
	655
	656
	657	def putData(self):
	658	#print("flagNoData",self.dataOut.flagNoData)
	659	#print("flagDiscontinuousBlock",self.dataOut.flagDiscontinuousBlock)
	660	#print(self.dataOut.flagNoData)
	661
	662	if self.dataOut.flagNoData:
	663	return 0
	664
	665	if self.dataOut.flagDiscontinuousBlock:
	666
	667	self.counter = 0
	668
	669	if self.counter == 0:
	670	self.setFile()
	671	#if self.experiment=="HP":
	672	#if self.dataOut.debris_activated==0:
	673	#self.writeBlock()
	674	#self.counter += 1
	675	#else:
	676	self.writeBlock()
	677	self.counter += 1
	678
	679	def close(self):
	680
	681	if self.counter > 0:
	682	self.fp.close()
	683	log.success('Closing file {}'.format(self.fullname), 'DatWriter')

schainpy/model/proc/jroproc_voltage_lags.py created 644 0 0

	1		NO CONTENT: new file 100644
The requested commit or file is too big and content was truncated. Show full diff

schainpy/model/data/jrodata.py +407 -8

		@@ -76,6 +76,8 def hildebrand_sekhon(data, navg):
76	76	"""
77	77
78	78	sortdata = numpy.sort(data, axis=None)
	79	#print(numpy.shape(data))
	80	#exit()
79	81	'''
80	82	lenOfData = len(sortdata)
81	83	nums_min = lenOfData*0.2
		@@ -273,13 +275,13 class JROData(GenericData):
273	275	'''
274	276	'''
275	277	return self.radarControllerHeaderObj.ippSeconds
276
	278
277	279	@ippSeconds.setter
278	280	def ippSeconds(self, ippSeconds):
279	281	'''
280	282	'''
281	283	self.radarControllerHeaderObj.ippSeconds = ippSeconds
282
	284
283	285	@property
284	286	def code(self):
285	287	'''
		@@ -370,7 +372,7 class Voltage(JROData):
370	372	self.flagShiftFFT = False
371	373	self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372	374	self.profileIndex = 0
373		self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
	375	self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
374	376	'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375	377
376	378	def getNoisebyHildebrand(self, channel=None):
		@@ -428,6 +430,103 class Voltage(JROData):
428	430	noise = property(getNoise, "I'm the 'nHeights' property.")
429	431
430	432
	433	class CrossProds(JROData):
	434
	435	# data es un numpy array de 2 dmensiones (canales, alturas)
	436	data = None
	437
	438	def __init__(self):
	439	'''
	440	Constructor
	441	'''
	442
	443	self.useLocalTime = True
	444	'''
	445	self.radarControllerHeaderObj = RadarControllerHeader()
	446	self.systemHeaderObj = SystemHeader()
	447	self.type = "Voltage"
	448	self.data = None
	449	# self.dtype = None
	450	# self.nChannels = 0
	451	# self.nHeights = 0
	452	self.nProfiles = None
	453	self.heightList = None
	454	self.channelList = None
	455	# self.channelIndexList = None
	456	self.flagNoData = True
	457	self.flagDiscontinuousBlock = False
	458	self.utctime = None
	459	self.timeZone = None
	460	self.dstFlag = None
	461	self.errorCount = None
	462	self.nCohInt = None
	463	self.blocksize = None
	464	self.flagDecodeData = False # asumo q la data no esta decodificada
	465	self.flagDeflipData = False # asumo q la data no esta sin flip
	466	self.flagShiftFFT = False
	467	self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
	468	self.profileIndex = 0
	469
	470
	471	def getNoisebyHildebrand(self, channel=None):
	472
	473
	474	if channel != None:
	475	data = self.data[channel]
	476	nChannels = 1
	477	else:
	478	data = self.data
	479	nChannels = self.nChannels
	480
	481	noise = numpy.zeros(nChannels)
	482	power = data * numpy.conjugate(data)
	483
	484	for thisChannel in range(nChannels):
	485	if nChannels == 1:
	486	daux = power[:].real
	487	else:
	488	daux = power[thisChannel, :].real
	489	noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
	490
	491	return noise
	492
	493	def getNoise(self, type=1, channel=None):
	494
	495	if type == 1:
	496	noise = self.getNoisebyHildebrand(channel)
	497
	498	return noise
	499
	500	def getPower(self, channel=None):
	501
	502	if channel != None:
	503	data = self.data[channel]
	504	else:
	505	data = self.data
	506
	507	power = data * numpy.conjugate(data)
	508	powerdB = 10 * numpy.log10(power.real)
	509	powerdB = numpy.squeeze(powerdB)
	510
	511	return powerdB
	512
	513	def getTimeInterval(self):
	514
	515	timeInterval = self.ippSeconds * self.nCohInt
	516
	517	return timeInterval
	518
	519	noise = property(getNoise, "I'm the 'nHeights' property.")
	520	timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
	521	'''
	522	def getTimeInterval(self):
	523
	524	timeInterval = self.ippSeconds * self.nCohInt
	525
	526	return timeInterval
	527
	528
	529
431	530	class Spectra(JROData):
432	531
433	532	def __init__(self):
		@@ -458,7 +557,7 class Spectra(JROData):
458	557	self.ippFactor = 1
459	558	self.beacon_heiIndexList = []
460	559	self.noise_estimation = None
461		self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
	560	self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
462	561	'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
463	562
464	563	def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
		@@ -608,7 +707,7 class Spectra(JROData):
608	707	print("This property should not be initialized")
609	708
610	709	return
611
	710
612	711	noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
613	712
614	713
		@@ -705,7 +804,7 class Fits(JROData):
705	804	return self.ipp_sec
706	805
707	806	noise = property(getNoise, "I'm the 'nHeights' property.")
708
	807
709	808
710	809	class Correlation(JROData):
711	810
		@@ -886,6 +985,7 class Parameters(Spectra):
886	985	else:
887	986	return self.paramInterval
888	987
	988
889	989	def setValue(self, value):
890	990
891	991	print("This property should not be initialized")
		@@ -928,6 +1028,10 class PlotterData(object):
928	1028	self.plottypes = ['cspc', 'spc', 'noise', 'rti']
929	1029	elif code == 'rti':
930	1030	self.plottypes = ['noise', 'rti']
	1031	elif code == 'crossprod':
	1032	self.plottypes = ['crossprod', 'kay']
	1033	elif code == 'spectrogram':
	1034	self.plottypes = ['spc', 'spectrogram']
931	1035	else:
932	1036	self.plottypes = [code]
933	1037
		@@ -976,9 +1080,11 class PlotterData(object):
976	1080	plot = 'snr'
977	1081	elif 'spc_moments' == plot:
978	1082	plot = 'moments'
	1083	elif 'spc_oblique' == plot:
	1084	plot = 'oblique'
979	1085	self.data[plot] = {}
980	1086
981		if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
	1087	if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data or 'oblique' in self.data:
982	1088	self.data['noise'] = {}
983	1089	self.data['rti'] = {}
984	1090	if 'noise' not in self.plottypes:
		@@ -1020,16 +1126,33 class PlotterData(object):
1020	1126	self.__heights.append(dataOut.heightList)
1021	1127	self.__all_heights.update(dataOut.heightList)
1022	1128
	1129
	1130
1023	1131	for plot in self.plottypes:
1024		if plot in ('spc', 'spc_moments', 'spc_cut'):
	1132	if plot in ('spc', 'spc_moments', 'spc_cut', 'spc_oblique'):
	1133
	1134
	1135	self.shift1 = dataOut.Oblique_params[0][1]
	1136	self.shift2 = dataOut.Oblique_params[0][4]
	1137	self.shift1_error = dataOut.Oblique_param_errors[0][1]
	1138	self.shift2_error = dataOut.Oblique_param_errors[0][4]
	1139
1025	1140	z = dataOut.data_spc/dataOut.normFactor
	1141	#print(dataOut.normFactor)
	1142	#print(z[0,3,15])
	1143	#print("here")
	1144	#print(dataOut.data_spc[0,0,0])
	1145	#exit()
1026	1146	buffer = 10*numpy.log10(z)
1027	1147	if plot == 'cspc':
1028	1148	buffer = (dataOut.data_spc, dataOut.data_cspc)
	1149	self.nFactor=dataOut.normFactor
1029	1150	if plot == 'noise':
1030	1151	buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1031	1152	if plot in ('rti', 'spcprofile'):
1032	1153	buffer = dataOut.getPower()
	1154	#print(buffer[0,0])
	1155	#exit()
1033	1156	if plot == 'snr_db':
1034	1157	buffer = dataOut.data_SNR
1035	1158	if plot == 'snr':
		@@ -1048,6 +1171,277 class PlotterData(object):
1048	1171	buffer = dataOut.data_output
1049	1172	if plot == 'param':
1050	1173	buffer = dataOut.data_param
	1174	if plot == 'spectrogram':
	1175	maxHei = 1350 #11
	1176	#maxHei = 2500
	1177	maxHei = 0
	1178	#maxHei = 990 #12
	1179	###maxHei = 990
	1180	indb = numpy.where(dataOut.heightList <= maxHei)
	1181	hei = indb[0][-1]
	1182	#hei = 19
	1183	print(hei)
	1184	#hei = 0
	1185	factor = dataOut.nIncohInt
	1186	#print(factor)
	1187
	1188	#exit(1)
	1189	z = dataOut.data_spc[:,:,hei] / factor
	1190
	1191	#for j in range(z.shape[1]):
	1192	#z[:,j] = z[:,j]/hildebrand_sekhon(z[], self.nCohInt)
	1193
	1194	##z = z/hildebrand_sekhon(z, factor)
	1195	noise = numpy.zeros(dataOut.nChannels)
	1196	for i in range(dataOut.nChannels):
	1197	#daux = numpy.sort(pair0[i,:,:],axis= None)
	1198	noise[i]=hildebrand_sekhon( z[i,:] ,dataOut.nIncohInt)
	1199	#for j in range(z.shape[1]):
	1200	#z[:,j] = z[:,j]/noise
	1201
	1202	#print(z.shape[1])
	1203	norm_factor = numpy.copy(z[:,int(z.shape[1]/2)])#/z[:,int(z.shape[1]/2)])*8000
	1204	#print(norm_factor)
	1205	#print(z[0,315:325])
	1206	#norm_factor = norm_factor.reshape((z.shape[0],z.shape[1]))
	1207	#print(norm_factor)
	1208	#exit(1)
	1209	#print(z.shape[1])
	1210
	1211	#for j in range(z.shape[1]):
	1212	#z[:,j] = z[:,j]/norm_factor
	1213
	1214	#print(z[0,315:325])
	1215	#exit(1)
	1216
	1217	#z = numpy.mean(dataOut.data_spc[:,:,:],axis=2) / factor
	1218	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
	1219	#avg = numpy.average(z, axis=1)
	1220	#print((dataOut.data_spc.shape))
	1221	#exit(1)
	1222	self.hei = hei
	1223	self.heightList = dataOut.heightList
	1224	self.DH = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
	1225	self.nProfiles = dataOut.nProfiles
	1226	#print(dataOut.heightList)
	1227
	1228
	1229	buffer = 10 * numpy.log10(z)
	1230
	1231
	1232	###buffer = z
	1233	import matplotlib.pyplot as plt
	1234	fig, axes = plt.subplots(figsize=(14, 10))
	1235	x = numpy.linspace(0,20,numpy.shape(buffer)[1])
	1236	x = numpy.fft.fftfreq(numpy.shape(buffer)[1],0.00001)
	1237	x = numpy.fft.fftshift(x)
	1238
	1239	plt.plot(x,buffer[0,:])
	1240	axes = plt.gca()
	1241	axes.set_xlim([-10000,10000])
	1242
	1243	#axes.set_xlim([0,30000])
	1244	#axes.set_ylim([-100,0.0025*1e10])
	1245	plt.show()
	1246	import time
	1247	#time.sleep(20)
	1248	#exit(1)
	1249
	1250
	1251
	1252	#if dataOut.profileIndex
	1253
	1254	if plot == 'xmit':
	1255	y_1=numpy.arctan2(dataOut.output_LP[:,0,2].imag,dataOut.output_LP[:,0,2].real)* 180 / (numpy.pi*10)
	1256	y_2=numpy.abs(dataOut.output_LP[:,0,2])
	1257	norm=numpy.max(y_2)
	1258	norm=max(norm,0.1)
	1259	y_2=y_2/norm
	1260
	1261	buffer = numpy.vstack((y_1,y_2))
	1262	self.NLAG = dataOut.NLAG
	1263
	1264	if plot == 'crossprod':
	1265	buffer = dataOut.crossprods
	1266	self.NDP = dataOut.NDP
	1267
	1268	if plot == 'crossprodlp':
	1269	buffer = 10*numpy.log10(numpy.abs(dataOut.output_LP))
	1270	self.NRANGE = dataOut.NRANGE
	1271	self.NLAG = dataOut.NLAG
	1272
	1273
	1274	if plot == 'noisedp':
	1275	buffer = 10*numpy.log10(dataOut.noise_final)
	1276	#print(buffer)
	1277
	1278	if plot == 'FaradayAngle':
	1279	buffer = numpy.degrees(dataOut.phi)
	1280	#print(buffer)
	1281
	1282	if plot == 'RTIDP':
	1283	buffer = dataOut.data_for_RTI_DP
	1284	self.NDP = dataOut.NDP
	1285
	1286	if plot == 'RTILP':
	1287	buffer = dataOut.data_for_RTI_LP
	1288	self.NRANGE = dataOut.NRANGE
	1289
	1290
	1291	if plot == 'denrti':
	1292	buffer = dataOut.DensityFinal
	1293
	1294
	1295	if plot == 'denrtiLP':
	1296
	1297	#buffer = numpy.reshape(numpy.concatenate((dataOut.ph2[:dataOut.cut],dataOut.ne[dataOut.cut:dataOut.NACF])),(1,-1))
	1298	buffer = dataOut.DensityFinal
	1299	#self.flagDataAsBlock = dataOut.flagDataAsBlock
	1300	#self.NDP = dataOut.NDP
	1301	if plot == 'den':
	1302	buffer = dataOut.ph2[:dataOut.NSHTS]
	1303	self.dphi=dataOut.dphi[:dataOut.NSHTS]
	1304	self.sdp2=dataOut.sdp2[:dataOut.NSHTS]
	1305	self.sdn1=dataOut.sdn1[:dataOut.NSHTS]#/self.dphi
	1306	self.NSHTS=dataOut.NSHTS
	1307	'''
	1308	flag1=False
	1309	flag0=True
	1310	for i in range(12,dataOut.NSHTS):
	1311	print("H: ",i*15)
	1312	print(abs((dataOut.sdn1[i]/(dataOut.dphi[i]*2))100))
	1313	if flag0:
	1314	if abs((dataOut.sdn1[i]/dataOut.dphi[i]))<0.0005*abs(dataOut.dphi[i]):
	1315	print("***************************** FIRST: ",(i)15,"****************************")
	1316	flag1=True
	1317	flag0=False
	1318	#pass
	1319	#print("**************************************GOOD**************************************")
	1320	#else:
	1321	#print("***************************************",(i-1)15,"****************************************")
	1322	#break
	1323	if flag1:
	1324	if abs((dataOut.sdn1[i]/dataOut.dphi[i]))>0.0005*abs(dataOut.dphi[i]):
	1325	print("***************************** LAST: ",(i-1)15,"****************************")
	1326	break
	1327	#print("H: ",i*15)
	1328	#print(dataOut.sdn1[i])
	1329	'''
	1330	if plot == 'denLP':
	1331	buffer = dataOut.ph2[:dataOut.NSHTS]
	1332	self.dphi=dataOut.dphi[:dataOut.NSHTS]
	1333	self.sdp2=dataOut.sdp2[:dataOut.NSHTS]
	1334	self.ne=dataOut.ne[:dataOut.NACF]
	1335	self.ene=dataOut.ene[:dataOut.NACF]dataOut.ne[:dataOut.NACF]0.434
	1336	#self.ene=10**dataOut.ene[:dataOut.NACF]
	1337	self.NSHTS=dataOut.NSHTS
	1338	self.cut=dataOut.cut
	1339
	1340	if plot == 'ETemp':
	1341	#buffer = dataOut.ElecTempClean
	1342	buffer = dataOut.ElecTempFinal
	1343	if plot == 'ITemp':
	1344	#buffer = dataOut.IonTempClean
	1345	buffer = dataOut.IonTempFinal
	1346	if plot == 'ETempLP':
	1347	#buffer = dataOut.IonTempClean
	1348	#buffer = numpy.reshape(numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:])),(1,-1))
	1349	buffer = dataOut.ElecTempFinal
	1350	#print(buffer)
	1351	if plot == 'ITempLP':
	1352	#buffer = dataOut.IonTempClean
	1353	#buffer = numpy.reshape(numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:])),(1,-1))
	1354	buffer = dataOut.IonTempFinal
	1355
	1356	if plot == 'HFracLP':
	1357	#buffer = dataOut.IonTempClean
	1358	#buffer = numpy.reshape(numpy.concatenate((dataOut.phy2[:dataOut.cut],dataOut.ph[dataOut.cut:])),(1,-1))
	1359	buffer = dataOut.PhyFinal
	1360	if plot == 'HeFracLP':
	1361	#buffer = dataOut.IonTempClean
	1362	#nan_array=numpy.empty((dataOut.cut))
	1363	#nan_array[:]=numpy.nan
	1364	#buffer = numpy.reshape(numpy.concatenate((nan_array,dataOut.phe[dataOut.cut:])),(1,-1))
	1365	buffer = dataOut.PheFinal
	1366
	1367
	1368
	1369
	1370
	1371	if plot =='acfs':
	1372	buffer = dataOut.acfs_to_plot
	1373	self.acfs_error_to_plot=dataOut.acfs_error_to_plot
	1374	self.lags_to_plot=dataOut.lags_to_plot
	1375	self.x_igcej_to_plot=dataOut.x_igcej_to_plot
	1376	self.x_ibad_to_plot=dataOut.x_ibad_to_plot
	1377	self.y_igcej_to_plot=dataOut.y_igcej_to_plot
	1378	self.y_ibad_to_plot=dataOut.y_ibad_to_plot
	1379	self.NSHTS = dataOut.NSHTS
	1380	self.DPL = dataOut.DPL
	1381	if plot =='acfs_LP':
	1382
	1383	aux=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
	1384	self.errors=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
	1385	self.lags_LP_to_plot=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
	1386	'''
	1387	for i in range(dataOut.NACF):
	1388	for j in range(dataOut.IBITS):
	1389	aux[i,j]=dataOut.fit_array_real[i,j]/dataOut.fit_array_real[i,0]
	1390	aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
	1391	'''
	1392	for i in range(dataOut.NACF):
	1393	for j in range(dataOut.IBITS):
	1394	if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0:
	1395	aux[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0]
	1396	aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
	1397	self.lags_LP_to_plot[i,j]=dataOut.lags_LP[j]
	1398	self.errors[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0]*dataOut.DH
	1399	else:
	1400	aux[i,j]=numpy.nan
	1401	self.lags_LP_to_plot[i,j]=numpy.nan
	1402	self.errors[i,j]=numpy.nan
	1403
	1404
	1405
	1406	buffer = aux
	1407
	1408	#self.lags_LP_to_plot=dataOut.lags_LP
	1409
	1410	self.NACF = dataOut.NACF
	1411	self.NLAG = dataOut.NLAG
	1412
	1413	if plot == 'tempsDP':
	1414
	1415	buffer = dataOut.te2
	1416	self.ete2 = dataOut.ete2
	1417	self.ti2 = dataOut.ti2
	1418	self.eti2 = dataOut.eti2
	1419
	1420	self.NSHTS = dataOut.NSHTS
	1421
	1422	if plot == 'temps_LP':
	1423
	1424	buffer = numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:]))
	1425	self.ete = numpy.concatenate((dataOut.ete2[:dataOut.cut],dataOut.ete[dataOut.cut:]))
	1426	self.ti = numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:]))
	1427	self.eti = numpy.concatenate((dataOut.eti2[:dataOut.cut],dataOut.eti[dataOut.cut:]))
	1428
	1429	self.NACF = dataOut.NACF
	1430
	1431
	1432	if plot == 'fracs_LP':
	1433
	1434	aux_nan=numpy.zeros(dataOut.cut,'float32')
	1435	aux_nan[:]=numpy.nan
	1436	buffer = numpy.concatenate((aux_nan,dataOut.ph[dataOut.cut:]))
	1437	self.eph = numpy.concatenate((aux_nan,dataOut.eph[dataOut.cut:]))
	1438	self.phe = dataOut.phe[dataOut.cut:]
	1439	self.ephe = dataOut.ephe[dataOut.cut:]
	1440
	1441	self.NACF = dataOut.NACF
	1442	self.cut = dataOut.cut
	1443
	1444
1051	1445	if plot == 'scope':
1052	1446	buffer = dataOut.data
1053	1447	self.flagDataAsBlock = dataOut.flagDataAsBlock
		@@ -1076,6 +1470,10 class PlotterData(object):
1076	1470	elif plot == 'spc_moments':
1077	1471	self.data['spc'][tm] = buffer
1078	1472	self.data['moments'][tm] = dataOut.moments
	1473	elif plot == 'spc_oblique':
	1474	self.data['spc'][tm] = buffer
	1475	self.data['shift1'][tm] = dataOut.Oblique_params[0]
	1476	self.data['shift2'][tm] = dataOut.Oblique_params[3]
1079	1477	else:
1080	1478	if self.buffering:
1081	1479	self.data[plot][tm] = buffer
		@@ -1141,6 +1539,7 class PlotterData(object):
1141	1539	meta['interval'] = float(self.interval)
1142	1540	meta['localtime'] = self.localtime
1143	1541	meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
	1542
1144	1543	if 'spc' in self.data or 'cspc' in self.data:
1145	1544	meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1146	1545	else:

schainpy/model/data/jroheaderIO.py +5 -2

                  timeZone = None
                  dstFlag = None
                  errorCount = None
-                 datatime = None
+                 F = None
                  structure = BASIC_STRUCTURE
                  __LOCALTIME = None
                      self.expType = int(header['nExpType'][0])
                      self.nTx = int(header['nNTx'][0])
                      self.ipp = float(header['fIpp'][0])
+                     #print(self.ipp)
                      self.txA = float(header['fTxA'][0])
                      self.txB = float(header['fTxB'][0])
                      self.nWindows = int(header['nNumWindows'][0])
                  def get_ippSeconds(self):
                      '''
                      '''
                      ippSeconds = 2.0 * 1000 * self.ipp / SPEED_OF_LIGHT
                      return ippSeconds
                      self.nWindows = int(header['nNumWindows'][0])
                      self.processFlags = header['nProcessFlags']
                      self.nCohInt = int(header['nCoherentIntegrations'][0])
                      self.nIncohInt = int(header['nIncoherentIntegrations'][0])
                      self.totalSpectra = int(header['nTotalSpectra'][0])
              def get_dtype_width(index):
-                 return DTYPE_WIDTH[index]
  No newline at end of file
+                 return DTYPE_WIDTH[index]

schainpy/model/graphics/__init__.py +1 0

              from .jroplot_heispectra import *
              from .jroplot_correlation import *
              from .jroplot_parameters import *
+             from .jroplot_voltage_lags import *

schainpy/model/graphics/jroplot_base.py +20 -16

                      self.zmin = kwargs.get('zmin', None)
                      self.zmax = kwargs.get('zmax', None)
                      self.zlimits = kwargs.get('zlimits', None)
+                     self.xlimits = kwargs.get('xlimits', None)
+                     self.xstep_given = kwargs.get('xstep_given', None)
+                     self.ystep_given = kwargs.get('ystep_given', None)
+                     self.autoxticks = kwargs.get('autoxticks', True)
                      self.xmin = kwargs.get('xmin', None)
                      self.xmax = kwargs.get('xmax', None)
                      self.xrange = kwargs.get('xrange', 12)
                      self.__throttle_plot = apply_throttle(self.throttle)
                      self.data = PlotterData(
                          self.CODE, self.throttle, self.exp_code, self.localtime, self.buffering, snr=self.showSNR)
                      if self.server:
                          if not self.server.startswith('tcp://'):
                              self.server = 'tcp://{}'.format(self.server)
                      self.setup()
-                     self.time_label = 'LT' if self.localtime else 'UTC'
+                     self.time_label = 'LT' if self.localtime else 'UTC'
                      if self.width is None:
                          self.width = 8
                      '''
                      Set min and max values, labels, ticks and titles
                      '''
                      for n, ax in enumerate(self.axes):
                          if ax.firsttime:
                              if self.xaxis != 'time':
                      self.plot()
                      self.format()
                      for n, fig in enumerate(self.figures):
                          if self.nrows == 0 or self.nplots == 0:
                              log.warning('No data', self.name)
                              fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
                              fig.canvas.manager.set_window_title(self.CODE)
                              continue
                          fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
                                                                               self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
                          fig.canvas.draw()
                          if self.save:
                              self.save_figure(n)
                      if self.server:
                          self.send_to_server()
                      figname = os.path.join(
                          self.save,
                          self.save_code,
-                         '{}_{}.png'.format(
+                         '{}_{}.png'.format(
                              self.save_code,
                              self.getDateTime(self.data.max_time).strftime(
                                  '%Y%m%d_%H%M%S'
                          return
                      self.sender_time = self.data.tm
                      attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
                      for attr in attrs:
                          value = getattr(self, attr)
                      except:
                          tm = self.sender_queue.get()
                          self.sender_queue.put(self.data.tm)
                      while True:
                          if self.sender_queue.empty():
                              break
                      self.ncols: number of cols
                      self.nplots: number of plots (channels or pairs)
                      self.ylabel: label for Y axes
-                     self.titles: list of axes title
+                     self.titles: list of axes title
                      '''
                      raise NotImplementedError
                      Must be defined in the child class
                      '''
                      raise NotImplementedError
                  def run(self, dataOut, **kwargs):
                      '''
                      Main plotting routine
                      '''
+                     print("time_inside_plot: ",dataOut.datatime)
+                     print(dataOut.flagNoData)
                      if self.isConfig is False:
                          self.__setup(**kwargs)
                              self.poll.register(self.socket, zmq.POLLIN)
                      tm = getattr(dataOut, self.attr_time)
-                     if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
+                     if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
                          self.save_time = tm
                          self.__plot()
                          self.tmin += self.xrange*60*60
                              dt = self.getDateTime(tm)
                              if self.xmin is None:
                                  self.tmin = tm
-                                 self.xmin = dt.hour
+                                 self.xmin = dt.hour
                              minutes = (self.xmin-int(self.xmin)) * 60
                              seconds = (minutes - int(minutes)) * 60
                              self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
                          self.__plot()
                      if self.data and not self.data.flagNoData and self.pause:
                          figpause(10)

schainpy/model/graphics/jroplot_parameters.py +11 -3

                  colormap = 'jet'
                  plot_type = 'pcolor'
+             class SpectralFitObliquePlot(SpectraPlot):
+                 '''
+                 Plot for Spectral Oblique
+                 '''
+                 CODE = 'spc_moments'
+                 colormap = 'jet'
+                 plot_type = 'pcolor'
              class SnrPlot(RTIPlot):
                  '''
                      self.nrows = self.data.shape(self.CODE)[0]
                      self.nplots = self.nrows
                      self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
                      if not self.xlabel:
                          self.xlabel = 'Time'
                      if self.showSNR:
                          self.nrows += 1
                          self.nplots += 1
                      self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
                      self.titles = ['{} {}'.format(
                          self.data.parameters[x], title) for x in self.channels]

schainpy/model/graphics/jroplot_spectra.py +551 -21

This diff has been collapsed as it changes many lines, (572 lines changed) Show them Hide them
			@@ -22,6 +22,7 class SpectraPlot(Plot):
	22	22	plot_type = 'pcolor'
	23	23
	24	24	def setup(self):
		25
	25	26	self.nplots = len(self.data.channels)
	26	27	self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
	27	28	self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
			@@ -31,10 +32,13 class SpectraPlot(Plot):
	31	32	self.width = 4 * self.ncols
	32	33	else:
	33	34	self.width = 3.5 * self.ncols
	34		self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
		35	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
	35	36	self.ylabel = 'Range [km]'
	36	37
	37	38	def plot(self):
		39
		40	#print(self.xaxis)
		41	#exit(1)
	38	42	if self.xaxis == "frequency":
	39	43	x = self.data.xrange[0]
	40	44	self.xlabel = "Frequency (kHz)"
			@@ -51,19 +55,25 class SpectraPlot(Plot):
	51	55
	52	56	self.titles = []
	53	57
		58
	54	59	y = self.data.heights
	55	60	self.y = y
	56	61	z = self.data['spc']
	57	62
		63	self.CODE2 = 'spc_oblique'
		64
		65
	58	66	for n, ax in enumerate(self.axes):
	59	67	noise = self.data['noise'][n][-1]
	60	68	if self.CODE == 'spc_moments':
	61	69	mean = self.data['moments'][n, :, 1, :][-1]
		70
	62	71	if ax.firsttime:
	63	72	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
	64	73	self.xmin = self.xmin if self.xmin else -self.xmax
	65	74	self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
	66	75	self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
		76	#print(numpy.shape(x))
	67	77	ax.plt = ax.pcolormesh(x, y, z[n].T,
	68	78	vmin=self.zmin,
	69	79	vmax=self.zmax,
			@@ -77,15 +87,122 class SpectraPlot(Plot):
	77	87	color="k", linestyle="dashed", lw=1)[0]
	78	88	if self.CODE == 'spc_moments':
	79	89	ax.plt_mean = ax.plot(mean, y, color='k')[0]
		90
	80	91	else:
		92
	81	93	ax.plt.set_array(z[n].T.ravel())
	82	94	if self.showprofile:
	83	95	ax.plt_profile.set_data(self.data['rti'][n][-1], y)
	84	96	ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
	85	97	if self.CODE == 'spc_moments':
	86	98	ax.plt_mean.set_data(mean, y)
		99
	87	100	self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
	88	101
		102	class SpectraObliquePlot(Plot):
		103	'''
		104	Plot for Spectra data
		105	'''
		106
		107	CODE = 'spc'
		108	colormap = 'jet'
		109	plot_type = 'pcolor'
		110
		111	def setup(self):
		112	self.xaxis = "oblique"
		113	self.nplots = len(self.data.channels)
		114	self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
		115	self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
		116	self.height = 2.6 * self.nrows
		117	self.cb_label = 'dB'
		118	if self.showprofile:
		119	self.width = 4 * self.ncols
		120	else:
		121	self.width = 3.5 * self.ncols
		122	self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
		123	self.ylabel = 'Range [km]'
		124
		125	def plot(self):
		126
		127	#print(self.xaxis)
		128	#exit(1)
		129	if self.xaxis == "frequency":
		130	x = self.data.xrange[0]
		131	self.xlabel = "Frequency (kHz)"
		132	elif self.xaxis == "time":
		133	x = self.data.xrange[1]
		134	self.xlabel = "Time (ms)"
		135	else:
		136	x = self.data.xrange[2]
		137	self.xlabel = "Velocity (m/s)"
		138
		139	if self.CODE == 'spc_moments':
		140	x = self.data.xrange[2]
		141	self.xlabel = "Velocity (m/s)"
		142
		143	self.titles = []
		144	#self.xlabel = "Velocidad (m/s)"
		145	#self.ylabel = 'Rango (km)'
		146
		147
		148	y = self.data.heights
		149	self.y = y
		150	z = self.data['spc']
		151
		152	self.CODE2 = 'spc_oblique'
		153
		154
		155	for n, ax in enumerate(self.axes):
		156	noise = self.data['noise'][n][-1]
		157	if self.CODE == 'spc_moments':
		158	mean = self.data['moments'][n, :, 1, :][-1]
		159	if self.CODE2 == 'spc_oblique':
		160	shift1 = self.data.shift1
		161	shift2 = self.data.shift2
		162	if ax.firsttime:
		163	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		164	self.xmin = self.xmin if self.xmin else -self.xmax
		165	self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
		166	self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
		167	#print(numpy.shape(x))
		168	ax.plt = ax.pcolormesh(x, y, z[n].T,
		169	vmin=self.zmin,
		170	vmax=self.zmax,
		171	cmap=plt.get_cmap(self.colormap)
		172	)
		173
		174	if self.showprofile:
		175	ax.plt_profile = self.pf_axes[n].plot(
		176	self.data['rti'][n][-1], y)[0]
		177	ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
		178	color="k", linestyle="dashed", lw=1)[0]
		179	if self.CODE == 'spc_moments':
		180	ax.plt_mean = ax.plot(mean, y, color='k')[0]
		181
		182	if self.CODE2 == 'spc_oblique':
		183	#ax.plt_shift1 = ax.plot(shift1, y, color='k', marker='x', linestyle='None', markersize=0.5)[0]
		184	#ax.plt_shift2 = ax.plot(shift2, y, color='m', marker='x', linestyle='None', markersize=0.5)[0]
		185	self.ploterr1 = ax.errorbar(shift1, y, xerr=self.data.shift1_error,fmt='k^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
		186	self.ploterr2 = ax.errorbar(shift2, y, xerr=self.data.shift2_error,fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
		187
		188	else:
		189	self.ploterr1.remove()
		190	self.ploterr2.remove()
		191	ax.plt.set_array(z[n].T.ravel())
		192	if self.showprofile:
		193	ax.plt_profile.set_data(self.data['rti'][n][-1], y)
		194	ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
		195	if self.CODE == 'spc_moments':
		196	ax.plt_mean.set_data(mean, y)
		197	if self.CODE2 == 'spc_oblique':
		198	#ax.plt_shift1.set_data(shift1, y)
		199	#ax.plt_shift2.set_data(shift2, y)
		200	#ax.clf()
		201	self.ploterr1 = ax.errorbar(shift1, y, xerr=self.data.shift1_error,fmt='k^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
		202	self.ploterr2 = ax.errorbar(shift2, y, xerr=self.data.shift2_error,fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
		203
		204	self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
		205	#self.titles.append('{}'.format('Velocidad Doppler'))
	89	206
	90	207	class CrossSpectraPlot(Plot):
	91	208
			@@ -103,7 +220,7 class CrossSpectraPlot(Plot):
	103	220	self.nrows = len(self.data.pairs)
	104	221	self.nplots = self.nrows * 4
	105	222	self.width = 3.1 * self.ncols
	106		self.height = ~~2.6~~ * self.nrows
		223	self.height = 5 * self.nrows
	107	224	self.ylabel = 'Range [km]'
	108	225	self.showprofile = False
	109	226	self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
			@@ -119,30 +236,44 class CrossSpectraPlot(Plot):
	119	236	else:
	120	237	x = self.data.xrange[2]
	121	238	self.xlabel = "Velocity (m/s)"
	122
		239
	123	240	self.titles = []
	124	241
		242
	125	243	y = self.data.heights
	126	244	self.y = y
	127	245	nspc = self.data['spc']
		246	#print(numpy.shape(self.data['spc']))
	128	247	spc = self.data['cspc'][0]
		248	#print(numpy.shape(spc))
		249	#exit()
	129	250	cspc = self.data['cspc'][1]
		251	#print(numpy.shape(cspc))
		252	#exit()
	130	253
	131	254	for n in range(self.nrows):
	132	255	noise = self.data['noise'][:,-1]
	133	256	pair = self.data.pairs[n]
		257	#print(pair)
		258	#exit()
	134	259	ax = self.axes[4 * n]
	135		if ax.firsttime:
	136		self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
	137		self.xmin = self.xmin if self.xmin else -self.xmax
	138		self.zmin = self.zmin if self.zmin else numpy.nanmin(~~nspc~~)
	139		self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
		260	if ax.firsttime:
		261	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		262	#self.xmin = self.xmin if self.xmin else -self.xmax
		263	self.xmin = self.xmin if self.xmin else numpy.nanmin(x)
		264	self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
		265	self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
		266	#print(numpy.nanmin(x))
		267	#print(self.xmax)
		268	#print(self.xmin)
		269	#exit()
		270	#self.xmin=-.1
	140	271	ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
	141	272	vmin=self.zmin,
	142	273	vmax=self.zmax,
	143	274	cmap=plt.get_cmap(self.colormap)
	144		)
	145		else:
		275	)
		276	else:
	146	277	ax.plt.set_array(nspc[pair[0]].T.ravel())
	147	278	self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
	148	279
			@@ -153,10 +284,10 class CrossSpectraPlot(Plot):
	153	284	vmax=self.zmax,
	154	285	cmap=plt.get_cmap(self.colormap)
	155	286	)
	156		else:
		287	else:
	157	288	ax.plt.set_array(nspc[pair[1]].T.ravel())
	158	289	self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
	159
		290
	160	291	out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
	161	292	coh = numpy.abs(out)
	162	293	phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
			@@ -172,19 +303,345 class CrossSpectraPlot(Plot):
	172	303	ax.plt.set_array(coh.T.ravel())
	173	304	self.titles.append(
	174	305	'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
	175
		306
	176	307	ax = self.axes[4 * n + 3]
	177	308	if ax.firsttime:
	178	309	ax.plt = ax.pcolormesh(x, y, phase.T,
	179	310	vmin=-180,
	180	311	vmax=180,
	181		cmap=plt.get_cmap(self.colormap_phase)
		312	cmap=plt.get_cmap(self.colormap_phase)
	182	313	)
	183	314	else:
	184	315	ax.plt.set_array(phase.T.ravel())
	185	316	self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
	186	317
	187	318
		319	class CrossSpectra4Plot(Plot):
		320
		321	CODE = 'cspc'
		322	colormap = 'jet'
		323	plot_type = 'pcolor'
		324	zmin_coh = None
		325	zmax_coh = None
		326	zmin_phase = None
		327	zmax_phase = None
		328
		329	def setup(self):
		330
		331	self.ncols = 4
		332	self.nrows = len(self.data.pairs)
		333	self.nplots = self.nrows * 4
		334	self.width = 3.1 * self.ncols
		335	self.height = 5 * self.nrows
		336	self.ylabel = 'Range [km]'
		337	self.showprofile = False
		338	self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
		339
		340	def plot(self):
		341
		342	if self.xaxis == "frequency":
		343	x = self.data.xrange[0]
		344	self.xlabel = "Frequency (kHz)"
		345	elif self.xaxis == "time":
		346	x = self.data.xrange[1]
		347	self.xlabel = "Time (ms)"
		348	else:
		349	x = self.data.xrange[2]
		350	self.xlabel = "Velocity (m/s)"
		351
		352	self.titles = []
		353
		354
		355	y = self.data.heights
		356	self.y = y
		357	nspc = self.data['spc']
		358	#print(numpy.shape(self.data['spc']))
		359	spc = self.data['cspc'][0]
		360	#print(numpy.shape(nspc))
		361	#exit()
		362	#nspc[1,:,:] = numpy.flip(nspc[1,:,:],axis=0)
		363	#print(numpy.shape(spc))
		364	#exit()
		365	cspc = self.data['cspc'][1]
		366
		367	#xflip=numpy.flip(x)
		368	#print(numpy.shape(cspc))
		369	#exit()
		370
		371	for n in range(self.nrows):
		372	noise = self.data['noise'][:,-1]
		373	pair = self.data.pairs[n]
		374	#print(pair)
		375	#exit()
		376	ax = self.axes[4 * n]
		377	if ax.firsttime:
		378	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		379	self.xmin = self.xmin if self.xmin else -self.xmax
		380	self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
		381	self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
		382	ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
		383	vmin=self.zmin,
		384	vmax=self.zmax,
		385	cmap=plt.get_cmap(self.colormap)
		386	)
		387	else:
		388	#print(numpy.shape(nspc[pair[0]].T))
		389	#exit()
		390	ax.plt.set_array(nspc[pair[0]].T.ravel())
		391	self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
		392
		393	ax = self.axes[4 * n + 1]
		394
		395	if ax.firsttime:
		396	ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
		397	vmin=self.zmin,
		398	vmax=self.zmax,
		399	cmap=plt.get_cmap(self.colormap)
		400	)
		401	else:
		402
		403	ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
		404	self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
		405
		406	out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
		407	coh = numpy.abs(out)
		408	phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
		409
		410	ax = self.axes[4 * n + 2]
		411	if ax.firsttime:
		412	ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
		413	vmin=0,
		414	vmax=1,
		415	cmap=plt.get_cmap(self.colormap_coh)
		416	)
		417	else:
		418	ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
		419	self.titles.append(
		420	'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
		421
		422	ax = self.axes[4 * n + 3]
		423	if ax.firsttime:
		424	ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
		425	vmin=-180,
		426	vmax=180,
		427	cmap=plt.get_cmap(self.colormap_phase)
		428	)
		429	else:
		430	ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
		431	self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
		432
		433
		434	class CrossSpectra2Plot(Plot):
		435
		436	CODE = 'cspc'
		437	colormap = 'jet'
		438	plot_type = 'pcolor'
		439	zmin_coh = None
		440	zmax_coh = None
		441	zmin_phase = None
		442	zmax_phase = None
		443
		444	def setup(self):
		445
		446	self.ncols = 1
		447	self.nrows = len(self.data.pairs)
		448	self.nplots = self.nrows * 1
		449	self.width = 3.1 * self.ncols
		450	self.height = 5 * self.nrows
		451	self.ylabel = 'Range [km]'
		452	self.showprofile = False
		453	self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
		454
		455	def plot(self):
		456
		457	if self.xaxis == "frequency":
		458	x = self.data.xrange[0]
		459	self.xlabel = "Frequency (kHz)"
		460	elif self.xaxis == "time":
		461	x = self.data.xrange[1]
		462	self.xlabel = "Time (ms)"
		463	else:
		464	x = self.data.xrange[2]
		465	self.xlabel = "Velocity (m/s)"
		466
		467	self.titles = []
		468
		469
		470	y = self.data.heights
		471	self.y = y
		472	#nspc = self.data['spc']
		473	#print(numpy.shape(self.data['spc']))
		474	#spc = self.data['cspc'][0]
		475	#print(numpy.shape(spc))
		476	#exit()
		477	cspc = self.data['cspc'][1]
		478	#print(numpy.shape(cspc))
		479	#exit()
		480
		481	for n in range(self.nrows):
		482	noise = self.data['noise'][:,-1]
		483	pair = self.data.pairs[n]
		484	#print(pair) #exit()
		485
		486
		487
		488	out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
		489
		490	#print(out[:,53])
		491	#exit()
		492	cross = numpy.abs(out)
		493	z = cross/self.data.nFactor
		494	#print("here")
		495	#print(dataOut.data_spc[0,0,0])
		496	#exit()
		497
		498	cross = 10*numpy.log10(z)
		499	#print(numpy.shape(cross))
		500	#print(cross[0,:])
		501	#print(self.data.nFactor)
		502	#exit()
		503	#phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
		504
		505	ax = self.axes[1 * n]
		506	if ax.firsttime:
		507	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		508	self.xmin = self.xmin if self.xmin else -self.xmax
		509	self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
		510	self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
		511	ax.plt = ax.pcolormesh(x, y, cross.T,
		512	vmin=self.zmin,
		513	vmax=self.zmax,
		514	cmap=plt.get_cmap(self.colormap)
		515	)
		516	else:
		517	ax.plt.set_array(cross.T.ravel())
		518	self.titles.append(
		519	'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
		520
		521
		522	class CrossSpectra3Plot(Plot):
		523
		524	CODE = 'cspc'
		525	colormap = 'jet'
		526	plot_type = 'pcolor'
		527	zmin_coh = None
		528	zmax_coh = None
		529	zmin_phase = None
		530	zmax_phase = None
		531
		532	def setup(self):
		533
		534	self.ncols = 3
		535	self.nrows = len(self.data.pairs)
		536	self.nplots = self.nrows * 3
		537	self.width = 3.1 * self.ncols
		538	self.height = 5 * self.nrows
		539	self.ylabel = 'Range [km]'
		540	self.showprofile = False
		541	self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
		542
		543	def plot(self):
		544
		545	if self.xaxis == "frequency":
		546	x = self.data.xrange[0]
		547	self.xlabel = "Frequency (kHz)"
		548	elif self.xaxis == "time":
		549	x = self.data.xrange[1]
		550	self.xlabel = "Time (ms)"
		551	else:
		552	x = self.data.xrange[2]
		553	self.xlabel = "Velocity (m/s)"
		554
		555	self.titles = []
		556
		557
		558	y = self.data.heights
		559	self.y = y
		560	#nspc = self.data['spc']
		561	#print(numpy.shape(self.data['spc']))
		562	#spc = self.data['cspc'][0]
		563	#print(numpy.shape(spc))
		564	#exit()
		565	cspc = self.data['cspc'][1]
		566	#print(numpy.shape(cspc))
		567	#exit()
		568
		569	for n in range(self.nrows):
		570	noise = self.data['noise'][:,-1]
		571	pair = self.data.pairs[n]
		572	#print(pair) #exit()
		573
		574
		575
		576	out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
		577
		578	#print(out[:,53])
		579	#exit()
		580	cross = numpy.abs(out)
		581	z = cross/self.data.nFactor
		582	cross = 10*numpy.log10(z)
		583
		584	out_r= out.real/self.data.nFactor
		585	#out_r = 10*numpy.log10(out_r)
		586
		587	out_i= out.imag/self.data.nFactor
		588	#out_i = 10*numpy.log10(out_i)
		589	#print(numpy.shape(cross))
		590	#print(cross[0,:])
		591	#print(self.data.nFactor)
		592	#exit()
		593	#phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
		594
		595	ax = self.axes[3 * n]
		596	if ax.firsttime:
		597	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		598	self.xmin = self.xmin if self.xmin else -self.xmax
		599	self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
		600	self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
		601	ax.plt = ax.pcolormesh(x, y, cross.T,
		602	vmin=self.zmin,
		603	vmax=self.zmax,
		604	cmap=plt.get_cmap(self.colormap)
		605	)
		606	else:
		607	ax.plt.set_array(cross.T.ravel())
		608	self.titles.append(
		609	'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
		610
		611	ax = self.axes[3 * n + 1]
		612	if ax.firsttime:
		613	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		614	self.xmin = self.xmin if self.xmin else -self.xmax
		615	self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
		616	self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
		617	ax.plt = ax.pcolormesh(x, y, out_r.T,
		618	vmin=-1.e6,
		619	vmax=0,
		620	cmap=plt.get_cmap(self.colormap)
		621	)
		622	else:
		623	ax.plt.set_array(out_r.T.ravel())
		624	self.titles.append(
		625	'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
		626
		627	ax = self.axes[3 * n + 2]
		628
		629
		630	if ax.firsttime:
		631	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
		632	self.xmin = self.xmin if self.xmin else -self.xmax
		633	self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
		634	self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
		635	ax.plt = ax.pcolormesh(x, y, out_i.T,
		636	vmin=-1.e6,
		637	vmax=1.e6,
		638	cmap=plt.get_cmap(self.colormap)
		639	)
		640	else:
		641	ax.plt.set_array(out_i.T.ravel())
		642	self.titles.append(
		643	'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
		644
	188	645	class RTIPlot(Plot):
	189	646	'''
	190	647	Plot for RTI data
			@@ -202,7 +659,7 class RTIPlot(Plot):
	202	659	self.ylabel = 'Range [km]'
	203	660	self.xlabel = 'Time'
	204	661	self.cb_label = 'dB'
	205		self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
		662	self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
	206	663	self.titles = ['{} Channel {}'.format(
	207	664	self.CODE.upper(), x) for x in range(self.nrows)]
	208	665
			@@ -210,6 +667,78 class RTIPlot(Plot):
	210	667	self.x = self.data.times
	211	668	self.y = self.data.heights
	212	669	self.z = self.data[self.CODE]
		670
		671	self.z = numpy.ma.masked_invalid(self.z)
		672
		673	if self.decimation is None:
		674	x, y, z = self.fill_gaps(self.x, self.y, self.z)
		675	else:
		676	x, y, z = self.fill_gaps(*self.decimate())
		677
		678	for n, ax in enumerate(self.axes):
		679	self.zmin = self.zmin if self.zmin else numpy.min(self.z)
		680	self.zmax = self.zmax if self.zmax else numpy.max(self.z)
		681	if ax.firsttime:
		682	ax.plt = ax.pcolormesh(x, y, z[n].T,
		683	vmin=self.zmin,
		684	vmax=self.zmax,
		685	cmap=plt.get_cmap(self.colormap)
		686	)
		687	if self.showprofile:
		688	ax.plot_profile = self.pf_axes[n].plot(
		689	self.data['rti'][n][-1], self.y)[0]
		690	ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
		691	color="k", linestyle="dashed", lw=1)[0]
		692	else:
		693	ax.collections.remove(ax.collections[0])
		694	ax.plt = ax.pcolormesh(x, y, z[n].T,
		695	vmin=self.zmin,
		696	vmax=self.zmax,
		697	cmap=plt.get_cmap(self.colormap)
		698	)
		699	if self.showprofile:
		700	ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
		701	ax.plot_noise.set_data(numpy.repeat(
		702	self.data['noise'][n][-1], len(self.y)), self.y)
		703
		704
		705	class SpectrogramPlot(Plot):
		706	'''
		707	Plot for Spectrogram data
		708	'''
		709
		710	CODE = 'spectrogram'
		711	colormap = 'binary'
		712	plot_type = 'pcolorbuffer'
		713
		714	def setup(self):
		715	self.xaxis = 'time'
		716	self.ncols = 1
		717	self.nrows = len(self.data.channels)
		718	self.nplots = len(self.data.channels)
		719	#print(self.dataOut.heightList)
		720	#self.ylabel = 'Range [km]'
		721	self.xlabel = 'Time'
		722	self.cb_label = 'dB'
		723	self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
		724	self.titles = ['{} Channel {} \n H = {} km ({} - {})'.format(
		725	self.CODE.upper(), x, self.data.heightList[self.data.hei], self.data.heightList[self.data.hei],self.data.heightList[self.data.hei]+(self.data.DH*self.data.nProfiles)) for x in range(self.nrows)]
		726
		727	def plot(self):
		728	self.x = self.data.times
		729	#self.y = self.data.heights
		730	self.z = self.data[self.CODE]
		731	self.y = self.data.xrange[0]
		732	#import time
		733	#print(time.ctime(self.x))
		734
		735	'''
		736	print(numpy.shape(self.x))
		737	print(numpy.shape(self.y))
		738	print(numpy.shape(self.z))
		739	'''
		740	self.ylabel = "Frequency (kHz)"
		741
	213	742	self.z = numpy.ma.masked_invalid(self.z)
	214	743
	215	744	if self.decimation is None:
			@@ -280,7 +809,7 class PhasePlot(CoherencePlot):
	280	809
	281	810	class NoisePlot(Plot):
	282	811	'''
	283		Plot for noise
		812	Plot for noise
	284	813	'''
	285	814
	286	815	CODE = 'noise'
			@@ -296,6 +825,7 class NoisePlot(Plot):
	296	825	self.xlabel = 'Time'
	297	826	self.titles = ['Noise']
	298	827	self.colorbar = False
		828	self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.17, 'right':0.95})
	299	829
	300	830	def plot(self):
	301	831
			@@ -315,7 +845,7 class NoisePlot(Plot):
	315	845	self.axes[0].lines[ch].set_data(x, y)
	316	846
	317	847	self.ymin = numpy.nanmin(Y) - 5
	318		self.ymax = numpy.nanmax(Y) + 5
		848	self.ymax = numpy.nanmax(Y) + 10
	319	849
	320	850
	321	851	class PowerProfilePlot(Plot):
			@@ -342,10 +872,10 class PowerProfilePlot(Plot):
	342	872	self.y = y
	343	873
	344	874	x = self.data['spcprofile']
	345
		875
	346	876	if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
	347	877	if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
	348
		878
	349	879	if self.axes[0].firsttime:
	350	880	for ch in self.data.channels:
	351	881	self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
			@@ -498,7 +1028,7 class BeaconPhase(Plot):
	498	1028	server=None, folder=None, username=None, password=None,
	499	1029	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
	500	1030
	501		if dataOut.flagNoData:
		1031	if dataOut.flagNoData:
	502	1032	return dataOut
	503	1033
	504	1034	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
			@@ -638,4 +1168,4 class BeaconPhase(Plot):
	638	1168	thisDatetime=thisDatetime,
	639	1169	update_figfile=update_figfile)
	640	1170
	641		return dataOut No newline at end of file
		1171	return dataOut

schainpy/model/io/__init__.py +6 -1

              from .julIO_param import *
              from .pxIO_param import *
-             from .jroIO_simulator import *
  No newline at end of file
+             from .jroIO_simulator import *
+             ############DP############
+             from .jroIO_dat import *
+             ############DP############

schainpy/model/io/jroIO_base.py +123 -83

                  basicHeaderObj = BasicHeader(LOCALTIME)
                  try:
                      fp = open(filename, 'rb')
                  except IOError:
                      print("The file %s can't be opened" % (filename))
                  firstBasicHeaderObj = BasicHeader(LOCALTIME)
                  systemHeaderObj = SystemHeader()
                  radarControllerHeaderObj = RadarControllerHeader()
                  processingHeaderObj = ProcessingHeader()
              def isRadarFile(file):
-                 try:
+                 try:
                      year = int(file[1:5])
                      doy = int(file[5:8])
                      set = int(file[8:11])
              def getDateFromRadarFile(file):
-                 try:
+                 try:
                      year = int(file[1:5])
                      doy = int(file[5:8])
-                     set = int(file[8:11])
+                     set = int(file[8:11])
                  except:
                      return None
                  return thisDate
              def parse_format(s, fmt):
                  for i in range(fmt.count('%')):
                      x = fmt.index('%')
                      d = DT_DIRECTIVES[fmt[x:x+2]]
                  def run(self):
-                     raise NotImplementedError
+                     raise NotImplementedError
                  def getAllowedArgs(self):
                      if hasattr(self, '__attrs__'):
                      for key, value in kwargs.items():
                          setattr(self, key, value)
                  def find_folders(self, path, startDate, endDate, folderfmt, last=False):
-                     folders = [x for f in path.split(',')
+                     folders = [x for f in path.split(',')
                                 for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))]
                      folders.sort()
                      if last:
                          folders = [folders[-1]]
-                     for folder in folders:
-                         try:
-                             dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
+                     for folder in folders:
+                         try:
+                             dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
                              if dt >= startDate and dt <= endDate:
                                  yield os.path.join(path, folder)
                              else:
                              log.log('Skiping folder {}'.format(folder), self.name)
                              continue
                      return
-                 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
+                 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
                                 expLabel='', last=False):
-                     for path in folders:
+                     for path in folders:
                          files = glob.glob1(path, '*{}'.format(ext))
                          files.sort()
                          if last:
-                             if files:
+                             if files:
                                  fo = files[-1]
-                                 try:
+                                 try:
                                      dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
-                                     yield os.path.join(path, expLabel, fo)
-                                 except Exception as e:
+                                     yield os.path.join(path, expLabel, fo)
+                                 except Exception as e:
                                      pass
                                  return
                              else:
                                  return
                          for fo in files:
-                             try:
-                                 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
+                             try:
+                                 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
+                                 #print(dt)
+                                 #print(startDate)
+                                 #print(endDate)
                                  if dt >= startDate and dt <= endDate:
                                      yield os.path.join(path, expLabel, fo)
                                  else:
                                      log.log('Skiping file {}'.format(fo), self.name)
                              except Exception as e:
                                  log.log('Skiping file {}'.format(fo), self.name)
-                                 continue
+                                 continue
                  def searchFilesOffLine(self, path, startDate, endDate,
-                                        expLabel, ext, walk,
+                                        expLabel, ext, walk,
                                         filefmt, folderfmt):
                      """Search files in offline mode for the given arguments
                              path, startDate, endDate, folderfmt)
                      else:
                          folders = path.split(',')
                      return self.find_files(
-                         folders, ext, filefmt, startDate, endDate, expLabel)
+                         folders, ext, filefmt, startDate, endDate, expLabel)
                  def searchFilesOnLine(self, path, startDate, endDate,
-                                       expLabel, ext, walk,
+                                       expLabel, ext, walk,
                                        filefmt, folderfmt):
                      """Search for the last file of the last folder
                      Return:
                          generator with the full path of last filename
                      """
                      if walk:
                          folders = self.find_folders(
                              path, startDate, endDate, folderfmt, last=True)
                      else:
                          folders = path.split(',')
                      return self.find_files(
                          folders, ext, filefmt, startDate, endDate, expLabel, last=True)
                  def setNextFile(self):
                      """Set the next file to be readed open it and parse de file header"""
+                     #print("fp: ",self.fp)
                      while True:
+                         #print(self.fp)
                          if self.fp != None:
-                             self.fp.close()
+                             self.fp.close()
+                             #print("setNextFile")
+                             #print("BEFORE OPENING",self.filename)
                          if self.online:
                              newFile = self.setNextFileOnline()
                          else:
                              newFile = self.setNextFileOffline()
+                             #print("newFile: ",newFile)
                          if not(newFile):
                              if self.online:
                                  raise schainpy.admin.SchainError('Time to wait for new files reach')
                              else:
                                  if self.fileIndex == -1:
+                                     #print("OKK")
                                      raise schainpy.admin.SchainWarning('No files found in the given path')
                                  else:
                                      raise schainpy.admin.SchainWarning('No more files to read')
                          if self.verifyFile(self.filename):
                              break
+                     ##print("BEFORE OPENING",self.filename)
                      log.log('Opening file: %s' % self.filename, self.name)
                      self.readFirstHeader()
                          self.filename
                          self.fp
                          self.filesize
                      Return:
                          boolean
                      """
                      nextFile = True
                      nextDay = False
-                     for nFiles in range(self.nFiles+1):
+                     for nFiles in range(self.nFiles+1):
                          for nTries in range(self.nTries):
                              fullfilename, filename = self.checkForRealPath(nextFile, nextDay)
                              if fullfilename is not None:
                                  self.name)
                              time.sleep(self.delay)
                              nextFile = False
-                             continue
+                             continue
                          if fullfilename is not None:
                              break
-                         self.nTries = 1
-                         nextFile = True
+                         #self.nTries = 1
+                         nextFile = True
                          if nFiles == (self.nFiles - 1):
                              log.log('Trying with next day...', self.name)
                              nextDay = True
-                             self.nTries = 3
+                             self.nTries = 3
                      if fullfilename:
                          self.fileSize = os.path.getsize(fullfilename)
                          self.flagIsNewFile = 1
                          if self.fp != None:
                              self.fp.close()
+                         #print(fullfilename)
                          self.fp = self.open_file(fullfilename, self.open_mode)
                          self.flagNoMoreFiles = 0
                          self.fileIndex += 1
                          return 1
-                     else:
+                     else:
                          return 0
                  def setNextFileOffline(self):
                      """Open the next file to be readed in offline mode"""
                      try:
                          filename = next(self.filenameList)
                          self.fileIndex +=1
                      except StopIteration:
                          self.flagNoMoreFiles = 1
-                         return 0
+                         return 0
+                     #print(self.fileIndex)
+                     #print(filename)
                      self.filename = filename
                      self.fileSize = os.path.getsize(filename)
                      self.fp = self.open_file(filename, self.open_mode)
                      self.flagIsNewFile = 1
                      return 1
                  @staticmethod
                  def isDateTimeInRange(dt, startDate, endDate, startTime, endTime):
                      """Check if the given datetime is in range"""
                      if startDate <= dt.date() <= endDate:
                          if startTime <= dt.time() <= endTime:
                              return True
                      return False
                  def verifyFile(self, filename):
                      """Check for a valid file
                      Arguments:
                          filename -- full path filename
                      Return:
                          boolean
                      """
                      """Check if the next file to be readed exists"""
                      raise NotImplementedError
                  def readFirstHeader(self):
                      """Parse the file header"""
                      pass
                  def waitDataBlock(self, pointer_location, blocksize=None):
                      Return:
                          str -- fullpath of the file
                      """
                      if nextFile:
                          self.set += 1
                      if nextDay:
                          prefixFileList = ['d', 'D']
                      elif self.ext.lower() == ".pdata":  # spectra
                          prefixFileList = ['p', 'P']
+                     ##############DP##############
+                     elif self.ext.lower() == ".dat":  # dat
+                         prefixFileList = ['z', 'Z']
+                     ##############DP##############
                      # barrido por las combinaciones posibles
                      for prefixDir in prefixDirList:
                          thispath = self.path
                              if os.path.exists(fullfilename):
                                  return fullfilename, filename
-                     return None, filename
+                     return None, filename
                  def __waitNewBlock(self):
                      """
                      Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
                              return 0
                          print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
+                         #print(self.filename)
                          time.sleep(self.delay)
                      return 0
                  def __setNewBlock(self):
                      if self.fp == None:
-                         return 0
-                     if self.flagIsNewFile:
+                         return 0
+                     if self.flagIsNewFile:
                          self.lastUTTime = self.basicHeaderObj.utc
                          return 1
                      currentSize = self.fileSize - self.fp.tell()
                      neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
                      if (currentSize >= neededSize):
                          self.basicHeaderObj.read(self.fp)
                          self.lastUTTime = self.basicHeaderObj.utc
                          return 1
                      if self.__waitNewBlock():
                          self.lastUTTime = self.basicHeaderObj.utc
                          return 1
                          print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
                                                                     self.processingHeaderObj.dataBlocksPerFile,
                                                                     self.dataOut.datatime.ctime()))
+                         #################DP#################
+                         self.dataOut.TimeBlockDate=self.dataOut.datatime.ctime()
+                         self.dataOut.TimeBlockSeconds=time.mktime(time.strptime(self.dataOut.datatime.ctime()))
+                         #################DP#################
                      return 1
                  def readFirstHeader(self):
                      except IOError:
                          log.error("File {} can't be opened".format(filename), self.name)
                          return False
                      if self.online and self.waitDataBlock(0):
                          pass
                      basicHeaderObj = BasicHeader(LOCALTIME)
                      systemHeaderObj = SystemHeader()
                      radarControllerHeaderObj = RadarControllerHeader()
                          dt2 = basicHeaderObj.datatime
                          if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
                                  self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
-                             flag = False
+                             flag = False
                      fp.close()
                      return flag
                      return dateList
                  def setup(self, **kwargs):
                      self.set_kwargs(**kwargs)
                      if not self.ext.startswith('.'):
                          self.ext = '.{}'.format(self.ext)
                      if self.server is not None:
                          if 'tcp://' in self.server:
                              address = server
                              for nTries in range(self.nTries):
                                  fullpath = self.searchFilesOnLine(self.path, self.startDate,
-                                     self.endDate, self.expLabel, self.ext, self.walk,
+                                     self.endDate, self.expLabel, self.ext, self.walk,
                                      self.filefmt, self.folderfmt)
                                  try:
                                      fullpath = next(fullpath)
                                  except:
                                      fullpath = None
                                  if fullpath:
                                      break
                                  log.warning(
                                      'Waiting {} sec for a valid file in {}: try {} ...'.format(
-                                         self.delay, self.path, nTries + 1),
+                                         self.delay, self.path, nTries + 1),
                                      self.name)
                                  time.sleep(self.delay)
                              if not(fullpath):
                                  raise schainpy.admin.SchainError(
-                                     'There isn\'t any valid file in {}'.format(self.path))
+                                     'There isn\'t any valid file in {}'.format(self.path))
                              pathname, filename = os.path.split(fullpath)
                              self.year = int(filename[1:5])
                              self.doy = int(filename[5:8])
-                             self.set = int(filename[8:11]) - 1
+                             self.set = int(filename[8:11]) - 1
                          else:
                              log.log("Searching files in {}".format(self.path), self.name)
-                             self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
+                             self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
                                  self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
                          self.setNextFile()
                      return
                      self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
                      self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
                  def getFirstHeader(self):
                      raise NotImplementedError
                      """
                      Arguments:
-                         path        :
-                         startDate   :
+                         path        :
+                         startDate   :
                          endDate     :
                          startTime   :
                          endTime     :
                      dtype_width = get_dtype_width(dtype_index)
                      return dtype_width
                  def getProcessFlags(self):
                      processFlags = 0
                      self.basicHeaderObj.size = self.basicHeaderSize  # bytes
                      self.basicHeaderObj.version = self.versionFile
-                     self.basicHeaderObj.dataBlock = self.nTotalBlocks
+                     self.basicHeaderObj.dataBlock = self.nTotalBlocks
                      utc = numpy.floor(self.dataOut.utctime)
-                     milisecond = (self.dataOut.utctime - utc) * 1000.0
+                     milisecond = (self.dataOut.utctime - utc) * 1000.0
                      self.basicHeaderObj.utc = utc
                      self.basicHeaderObj.miliSecond = milisecond
                      self.basicHeaderObj.timeZone = self.dataOut.timeZone
                      if self.dataOut.datatime.date() > self.fileDate:
                          setFile = 0
                          self.nTotalBlocks = 0
                      filen = '{}{:04d}{:03d}{:03d}{}'.format(
-                             self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
+                             self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
                      filename = os.path.join(path, subfolder, filen)
                      self.ext = ext.lower()
                      self.path = path
                      if set is None:
                          self.setFile = -1
                      else:
-                         self.setFile = set - 1
+                         self.setFile = set - 1
                      self.blocksPerFile = blocksPerFile
                      self.profilesPerBlock = profilesPerBlock

schainpy/model/io/jroIO_madrigal.py +104 -70

                  'sciRemarks': '',
                  'instRemarks': ''
                  }
              DEF_HEADER = {
                  'kindatDesc': '',
                  'analyst': 'Jicamarca User',
                          for k, v in list(iterable.items())}
                  elif isinstance(iterable, (list, tuple)):
                      return [str(v) if isinstance(v, basestring) else v for v in iterable]
                  return iterable
                      ProcessingUnit.__init__(self)
-                     self.dataOut = Parameters()
+                     self.dataOut = Parameters()
                      self.counter_records = 0
                      self.nrecords = None
                      self.flagNoMoreFiles = 0
-                     self.filename = None
+                     self.filename = None
                      self.intervals = set()
                      self.datatime = datetime.datetime(1900,1,1)
                      self.format = None
                      self.filefmt = "***%Y%m%d*******"
                  def setup(self, **kwargs):
                      self.set_kwargs(**kwargs)
                      self.oneDDict = load_json(self.oneDDict)
                      self.twoDDict = load_json(self.twoDDict)
                          for nTries in range(self.nTries):
                              fullpath = self.searchFilesOnLine(self.path, self.startDate,
-                                 self.endDate, self.expLabel, self.ext, self.walk,
+                                 self.endDate, self.expLabel, self.ext, self.walk,
                                  self.filefmt, self.folderfmt)
                              try:
                                  fullpath = next(fullpath)
                              except:
                                  fullpath = None
                              if fullpath:
                                  break
                              log.warning(
                                  'Waiting {} sec for a valid file in {}: try {} ...'.format(
-                                     self.delay, self.path, nTries + 1),
+                                     self.delay, self.path, nTries + 1),
                                  self.name)
                              time.sleep(self.delay)
                          if not(fullpath):
                              raise schainpy.admin.SchainError(
-                                 'There isn\'t any valid file in {}'.format(self.path))
+                                 'There isn\'t any valid file in {}'.format(self.path))
                      else:
                          log.log("Searching files in {}".format(self.path), self.name)
-                         self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
+                         self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
                              self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
                      self.setNextFile()
                  def readFirstHeader(self):
                      self.parseHeader()
                      self.parseData()
                      self.blockIndex = 0
-                     return
+                     return
                  def parseHeader(self):
                      '''
                      if s_parameters:
                          log.success('Spatial parameters found: {}'.format(s_parameters),
                                      'MADReader')
                      for param in list(self.oneDDict.keys()):
                          if param.lower() not in self.parameters:
                              log.warning(
                                      param),
                                  'MADReader')
                              self.oneDDict.pop(param, None)
                      for param, value in list(self.twoDDict.items()):
                          if param.lower() not in self.parameters:
                              log.warning(
                      while True:
                          self.flagDiscontinuousBlock = 0
                          if self.counter_records == self.nrecords:
-                             self.setNextFile()
+                             self.setNextFile()
                          self.readBlock()
                          if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
                             (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
                              log.warning(
                                  if self.counter_records == self.nrecords:
                                      break
                                  continue
-                             self.intervals.add((datatime-self.datatime).seconds)
+                             self.intervals.add((datatime-self.datatime).seconds)
                              break
                      elif self.ext == '.hdf5':
                          datatime = datetime.datetime.utcfromtimestamp(
                          if datatime.date()>self.datatime.date():
                              self.flagDiscontinuousBlock = 1
                          self.datatime = datatime
-                         self.counter_records += 1
+                         self.counter_records += 1
                      self.buffer = numpy.array(dum)
                      return
                  def set_output(self):
                      '''
                      Storing data from buffer to dataOut object
-                     '''
+                     '''
                      parameters = [None for __ in self.parameters]
-                     for param, attr in list(self.oneDDict.items()):
+                     for param, attr in list(self.oneDDict.items()):
                          x = self.parameters.index(param.lower())
                          setattr(self.dataOut, attr, self.buffer[0][x])
                      for param, value in list(self.twoDDict.items()):
-                         dummy = numpy.zeros(self.ranges.shape) + numpy.nan
+                         dummy = numpy.zeros(self.ranges.shape) + numpy.nan
                          if self.ext == '.txt':
                              x = self.parameters.index(param.lower())
-                             y = self.parameters.index(self.independentParam.lower())
+                             y = self.parameters.index(self.independentParam.lower())
                              ranges = self.buffer[:,y]
                              #if self.ranges.size == ranges.size:
                              #    continue
                              ranges = self.buffer[self.independentParam.lower()]
                              index = numpy.where(numpy.in1d(self.ranges, ranges))[0]
                              dummy[index] = self.buffer[param.lower()]
                          if isinstance(value, str):
-                             if value not in self.independentParam:
+                             if value not in self.independentParam:
                                  setattr(self.dataOut, value, dummy.reshape(1,-1))
-                         elif isinstance(value, list):
+                         elif isinstance(value, list):
                              self.output[value[0]][value[1]] = dummy
                              parameters[value[1]] = param
                      for key, value in list(self.output.items()):
                          setattr(self.dataOut, key, numpy.array(value))
                      self.dataOut.parameters = [s for s in parameters if s]
                      self.dataOut.heightList = self.ranges
                      self.dataOut.utctime = (self.datatime - datetime.datetime(1970, 1, 1)).total_seconds()
-                     self.dataOut.utctimeInit = self.dataOut.utctime
+                     self.dataOut.utctimeInit = self.dataOut.utctime
                      self.dataOut.paramInterval = min(self.intervals)
-                     self.dataOut.useLocalTime = False
-                     self.dataOut.flagNoData = False
+                     self.dataOut.useLocalTime = False
+                     self.dataOut.flagNoData = False
                      self.dataOut.nrecords = self.nrecords
                      self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
              @MPDecorator
              class MADWriter(Operation):
                  '''Writing module for Madrigal files
              type: external
              Inputs:
                  __attrs__ = ['path', 'oneDDict', 'ind2DList', 'twoDDict','metadata', 'format', 'blocks']
                  missing = -32767
                  def __init__(self):
                      Operation.__init__(self)
                  def run(self, dataOut, path, oneDDict, ind2DList='[]', twoDDict='{}',
                          metadata='{}', format='cedar', **kwargs):
+                     #if dataOut.AUX==1: #Modified
                      if not self.isConfig:
                          self.setup(path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs)
                          self.isConfig = True
-                     self.dataOut = dataOut
-                     self.putData()
+                     self.dataOut = dataOut
+                     self.putData()
                      return 1
                  def setup(self, path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs):
                      '''
-                     Configure Operation
+                     Configure Operation
                      '''
                      self.path = path
                      self.blocks = kwargs.get('blocks', None)
                      self.counter = 0
                      self.oneDDict = load_json(oneDDict)
                      self.twoDDict = load_json(twoDDict)
                      self.ind2DList = load_json(ind2DList)
-                     meta = load_json(metadata)
+                     meta = load_json(metadata)
                      self.kinst = meta.get('kinst')
                      self.kindat = meta.get('kindat')
                      self.catalog = meta.get('catalog', DEF_CATALOG)
                      elif format == 'hdf5':
                          self.ext = '.hdf5'
                          self.extra_args = {'ind2DList': self.ind2DList}
-                     self.keys = [k.lower() for k in self.twoDDict]
+                     self.keys = [k.lower() for k in self.twoDDict]
                      if 'range' in self.keys:
                          self.keys.remove('range')
                      if 'gdalt' in self.keys:
                      self.mnemonic = MNEMONICS[self.kinst]   #TODO get mnemonic from madrigal
                      date = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
+                     #if self.dataOut.input_dat_type:
+                         #date=datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds_for_dp_power)
+                     #print("date",date)
                      filename = '{}{}{}'.format(self.mnemonic,
                                                 date.strftime('%Y%m%d_%H%M%S'),
                                                 self.ext)
                      self.fullname = os.path.join(self.path, filename)
-                     if os.path.isfile(self.fullname) :
+                     if os.path.isfile(self.fullname) :
                          log.warning(
                              'Destination file {} already exists, previous file deleted.'.format(
                                  self.fullname),
                              'MADWriter')
                          os.remove(self.fullname)
                      try:
                          log.success(
                              'Creating file: {}'.format(self.fullname),
                          if not os.path.exists(self.path):
                              os.makedirs(self.path)
                          self.fp = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
                      except ValueError as e:
                          log.error(
                              'Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile"',
                      attributes.
                      Allowed parameters in: parcodes.tab
                      '''
+                     #self.dataOut.paramInterval=2
                      startTime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
                      endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
+                     #if self.dataOut.input_dat_type:
+                     #if self.dataOut.experiment=="DP":
+                         #startTime=datetime.datetime.fromtimestamp(self.dataOut.TimeBlockSeconds_for_dp_power)
+                         #endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval)
+                     #print("2: ",startTime)
+                     #print(endTime)
                      heights = self.dataOut.heightList
+                     #print(self.blocks)
+                     #print(startTime)
+                     #print(endTime)
+                     #print(heights)
+                     #input()
                      if self.ext == '.dat':
                          for key, value in list(self.twoDDict.items()):
                              if isinstance(value, str):
                              out[key] = tmp.flatten()[:len(heights)]
                          elif isinstance(value, (tuple, list)):
                              attr, x = value
-                             data = getattr(self.dataOut, attr)
+                             data = getattr(self.dataOut, attr)
+                             #print(x)
+                             #print(len(heights))
+                             #print(data[int(x)][:len(heights)])
+                             #print(numpy.shape(out))
+                             #print(numpy.shape(data))
                              out[key] = data[int(x)][:len(heights)]
                      a = numpy.array([out[k] for k in self.keys])
+                     #print(a)
                      nrows = numpy.array([numpy.isnan(a[:, x]).all() for x in range(len(heights))])
                      index = numpy.where(nrows == False)[0]
+                     #print(startTime.minute)
                      rec = madrigal.cedar.MadrigalDataRecord(
                          self.kinst,
                          self.kindat,
                          len(index),
                          **self.extra_args
                      )
-                     # Setting 1d values
+                     #print("rec",rec)
+                     # Setting 1d values
                      for key in self.oneDDict:
                          rec.set1D(key, getattr(self.dataOut, self.oneDDict[key]))
                      # Setting 2d values
                      nrec = 0
-                     for n in index:
+                     for n in index:
                          for key in out:
                              rec.set2D(key, nrec, out[key][n])
-                         nrec += 1
+                         nrec += 1
                      self.fp.append(rec)
-                     if self.ext == '.hdf5' and self.counter % 500 == 0 and self.counter > 0:
+                     if self.ext == '.hdf5' and self.counter %2 == 0 and self.counter > 0:
+                         #print("here")
                          self.fp.dump()
                      if self.counter % 20 == 0 and self.counter > 0:
+                         #self.fp.write()
                          log.log(
                              'Writing {} records'.format(
                                  self.counter),
                  def setHeader(self):
                      '''
                      Create an add catalog and header to cedar file
-                     '''
+                     '''
                      log.success('Closing file {}'.format(self.fullname), 'MADWriter')
                      if self.ext == '.dat':
                      else:
                          self.fp.dump()
                          self.fp.close()
-                     header = madrigal.cedar.CatalogHeaderCreator(self.fullname)
+                     header = madrigal.cedar.CatalogHeaderCreator(self.fullname)
                      header.createCatalog(**self.catalog)
                      header.createHeader(**self.header)
                      header.write()
                  def putData(self):
                      if self.dataOut.flagNoData:
-                         return 0
+                         return 0
                      if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks:
                          if self.counter > 0:
                              self.setHeader()
                      if self.counter == 0:
                          self.setFile()
                      self.writeBlock()
-                     self.counter += 1
+                     self.counter += 1
                  def close(self):
-                     if self.counter > 0:
-                         self.setHeader()
  No newline at end of file
+                     if self.counter > 0:
+                         self.setHeader()

schainpy/model/io/jroIO_param.py +41 -40

                  This unit reads HDF5 files created with `HDFWriter` operation contains
                  by default two groups Data and Metadata all variables would be saved as `dataOut`
-                 attributes.
+                 attributes.
                  It is possible to read any HDF5 file by given the structure in the `description`
                  parameter, also you can add extra values to metadata with the parameter `extras`.
                      Dictionary with the description of the HDF5 file
                  extras : dict, optional
                      Dictionary with extra metadata to be be added to `dataOut`
                  Examples
                  --------
                  desc = {
                      'Data': {
                          'data_output': ['u', 'v', 'w'],
                  extras = {
                      'timeZone': 300
                  }
                  reader = project.addReadUnit(
                      name='HDFReader',
                      path='/path/to/files',
                      self.set_kwargs(**kwargs)
                      if not self.ext.startswith('.'):
-                         self.ext = '.{}'.format(self.ext)
+                         self.ext = '.{}'.format(self.ext)
                      if self.online:
                          log.log("Searching files in online mode...", self.name)
                          for nTries in range(self.nTries):
                              fullpath = self.searchFilesOnLine(self.path, self.startDate,
-                                 self.endDate, self.expLabel, self.ext, self.walk,
+                                 self.endDate, self.expLabel, self.ext, self.walk,
                                  self.filefmt, self.folderfmt)
                              try:
                                  fullpath = next(fullpath)
                              except:
                                  fullpath = None
                              if fullpath:
                                  break
                              log.warning(
                                  'Waiting {} sec for a valid file in {}: try {} ...'.format(
-                                     self.delay, self.path, nTries + 1),
+                                     self.delay, self.path, nTries + 1),
                                  self.name)
                              time.sleep(self.delay)
                          if not(fullpath):
                              raise schainpy.admin.SchainError(
-                                 'There isn\'t any valid file in {}'.format(self.path))
+                                 'There isn\'t any valid file in {}'.format(self.path))
                          pathname, filename = os.path.split(fullpath)
                          self.year = int(filename[1:5])
                          self.doy = int(filename[5:8])
-                         self.set = int(filename[8:11]) - 1
+                         self.set = int(filename[8:11]) - 1
                      else:
                          log.log("Searching files in {}".format(self.path), self.name)
-                         self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
+                         self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
                              self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
                      self.setNextFile()
                      return
                  def readFirstHeader(self):
                      '''Read metadata and data'''
-                     self.__readMetadata()
+                     self.__readMetadata()
                      self.__readData()
                      self.__setBlockList()
                      if 'type' in self.meta:
                          self.dataOut = eval(self.meta['type'])()
                      for attr in self.meta:
                          setattr(self.dataOut, attr, self.meta[attr])
                      self.blockIndex = 0
                      return
                  def __setBlockList(self):
                      else:
                          grp = self.fp['Metadata']
                          for name in grp:
-                             meta[name] = grp[name].value
+                             meta[name] = grp[name].value
                      if self.extras:
                          for key, value in self.extras.items():
                  def __readData(self):
                      data = {}
                      if self.description:
                          for key, value in self.description['Data'].items():
                              if isinstance(value, str):
                                  array = numpy.array(array)
                              else:
                                  log.warning('Unknown type: {}'.format(name))
                              if name in self.description:
                                  key = self.description[name]
                              else:
                      self.data = data
                      return
                  def getData(self):
                      for attr in self.data:
                  The HDF5 file contains by default two groups Data and Metadata where
                  you can save any `dataOut` attribute specified by `dataList` and `metadataList`
                  parameters, data attributes are normaly time dependent where the metadata
-                 are not.
-                 It is possible to customize the structure of the HDF5 file with the
+                 are not.
+                 It is possible to customize the structure of the HDF5 file with the
                  optional description parameter see the examples.
                  Parameters:
                      If True the name of the files corresponds to the timestamp of the data
                  description : dict, optional
                      Dictionary with the desired description of the HDF5 file
                  Examples
                  --------
                  desc = {
                      'data_output': {'winds': ['z', 'w', 'v']},
                      'utctime': 'timestamps',
                          'heightList': 'heights'
                      }
                  }
                  writer = proc_unit.addOperation(name='HDFWriter')
                  writer.addParameter(name='path', value='/path/to/file')
                  writer.addParameter(name='blocksPerFile', value='32')
                  lastTime = None
                  def __init__(self):
                      Operation.__init__(self)
                      return
                              dsDict['shape'] = dataAux.shape
                              dsDict['dsNumber'] = dataAux.shape[0]
                              dsDict['dtype'] = dataAux.dtype
                          dsList.append(dsDict)
                      self.dsList = dsList
                          self.lastTime = currentTime
                          self.currentDay = dataDay
                          return False
                      timeDiff = currentTime - self.lastTime
                      #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
                  def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
                          dataList=[], setType=None, description={}):
+                     print("hdf",dataOut.flagNoData)
+                     print(dataOut.datatime.ctime())
                      self.dataOut = dataOut
                      if not(self.isConfig):
-                         self.setup(path=path, blocksPerFile=blocksPerFile,
+                         self.setup(path=path, blocksPerFile=blocksPerFile,
                                     metadataList=metadataList, dataList=dataList,
                                     setType=setType, description=description)
                      self.putData()
                      return
                  def setNextFile(self):
                      ext = self.ext
                      path = self.path
                      setFile = self.setFile
                              return 'pair{:02d}'.format(x)
                          else:
                              return 'channel{:02d}'.format(x)
                  def writeMetadata(self, fp):
                      if self.description:
                      return
                  def writeData(self, fp):
                      if self.description:
                          if 'Data' in self.description:
                              grp = fp.create_group('Data')
                      dtsets = []
                      data = []
                      for dsInfo in self.dsList:
                          if dsInfo['nDim'] == 0:
                              ds = grp.create_dataset(
-                                 self.getLabel(dsInfo['variable']),
+                                 self.getLabel(dsInfo['variable']),
                                  (self.blocksPerFile, ),
-                                 chunks=True,
+                                 chunks=True,
                                  dtype=numpy.float64)
                              dtsets.append(ds)
                              data.append((dsInfo['variable'], -1))
                                  sgrp = grp
                              for i in range(dsInfo['dsNumber']):
                                  ds = sgrp.create_dataset(
-                                     self.getLabel(dsInfo['variable'], i),
+                                     self.getLabel(dsInfo['variable'], i),
                                      (self.blocksPerFile, ) + dsInfo['shape'][1:],
                                      chunks=True,
                                      dtype=dsInfo['dtype'])
                      fp.flush()
                      log.log('Creating file: {}'.format(fp.filename), self.name)
                      self.ds = dtsets
                      self.data = data
                      self.firsttime = True

schainpy/model/io/jroIO_voltage.py +26 -10

                      """
                      ProcessingUnit.__init__(self)
                      self.ext = ".r"
                      self.optchar = "D"
                      self.basicHeaderObj = BasicHeader(LOCALTIME)
                      self.systemHeaderObj = SystemHeader()
                      self.radarControllerHeaderObj = RadarControllerHeader()
                      self.processingHeaderObj = ProcessingHeader()
                      self.lastUTTime = 0
-                     self.profileIndex = 2**32 - 1
+                     self.profileIndex = 2**32 - 1
                      self.dataOut = Voltage()
                      self.selBlocksize = None
                      self.selBlocktime = None
+                     ##print("1--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                  def createObjByDefault(self):
+                     ##print("2--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      dataObj = Voltage()
                      return dataObj
                  def __hasNotDataInBuffer(self):
+                     ##print("3--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      if self.profileIndex >= self.processingHeaderObj.profilesPerBlock * self.nTxs:
                          return 1
                          Return:
                              None
                      """
+                     ##print("4--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      pts2read = self.processingHeaderObj.profilesPerBlock * \
                          self.processingHeaderObj.nHeights * self.systemHeaderObj.nChannels
                      self.blocksize = pts2read
                  def readBlock(self):
                      """
                          readBlock lee el bloque de datos desde la posicion actual del puntero del archivo
                          (self.fp) y actualiza todos los parametros relacionados al bloque de datos
                              self.flagIsNewBlock
                              self.nTotalBlocks
-                         Exceptions:
+                         Exceptions:
                              Si un bloque leido no es un bloque valido
                      """
+                     ##print("5--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      # if self.server is not None:
                      #     self.zBlock = self.receiver.recv()
                      #     self.zHeader = self.zBlock[:24]
                      return 1
                  def getFirstHeader(self):
+                     ##print("6--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      self.getBasicHeader()
                      self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
+                     #self.dataOut.ippSeconds_general=self.radarControllerHeaderObj.ippSeconds
+                     #print(self.nTxs)
                      if self.nTxs > 1:
+                         #print(self.radarControllerHeaderObj.ippSeconds)
                          self.dataOut.radarControllerHeaderObj.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
+                         #print(self.radarControllerHeaderObj.ippSeconds)
                      # Time interval and code are propierties of dataOut. Its value depends of radarControllerHeaderObj.
                      #         self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt
                      self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
                  def reshapeData(self):
+                     ##print("7--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      if self.nTxs < 0:
                          return
                  def readFirstHeaderFromServer(self):
+                     ##print("8--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      self.getFirstHeader()
                      self.firstHeaderSize = self.basicHeaderObj.size
                      self.getBlockDimension()
                  def getFromServer(self):
+                     ##print("9--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      self.flagDiscontinuousBlock = 0
                      self.profileIndex = 0
                      self.flagIsNewBlock = 1
                              self.flagDiscontinuousBlock
                              self.flagIsNewBlock
                      """
+                     ##print("10--OKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK")
                      if self.flagNoMoreFiles:
                          self.dataOut.flagNoData = True
                          return 0
                          self.dataOut.data = self.datablock[:, self.profileIndex, :]
                          self.dataOut.profileIndex = self.profileIndex
                          self.profileIndex += 1
                      else:
                          self.dataOut.flagDataAsBlock = True
                          self.dataOut.nProfiles = self.dataOut.data.shape[1]
+             #######################DP#######################
+                     self.dataOut.CurrentBlock=self.nReadBlocks
+                     self.dataOut.LastBlock=self.processingHeaderObj.dataBlocksPerFile
+             #######################DP#######################
                      self.dataOut.flagNoData = False
-                     self.getBasicHeader()
+                     #self.getBasicHeader()
                      self.dataOut.realtime = self.online
                      self.processingHeaderObj.processFlags = self.getProcessFlags()
                      self.setBasicHeader()
-  No newline at end of file

schainpy/model/proc/__init__.py +6 0

              from .jroproc_spectra_acf import *
              from .bltrproc_parameters import *
              from .pxproc_parameters import *
+             ###########DP###########
+             from .jroproc_voltage_lags import *
+             ###########DP###########
+             from .jroproc_spectra_lags_faraday import *

schainpy/model/proc/jroproc_base.py +1 -1

                          self.op_type = 'external'
                          self.name = BaseClass.__name__
                          self.__doc__ = BaseClass.__doc__
                          if 'plot' in self.name.lower() and not self.name.endswith('_'):
                              self.name = '{}{}'.format(self.CODE.upper(), 'Plot')

schainpy/model/proc/jroproc_parameters.py 755 ➡ 644 +535 -169

This diff has been collapsed as it changes many lines, (704 lines changed) Show them Hide them
			@@ -302,6 +302,12 class SpectralFilters(Operation):
	302	302	dataOut.spcparam_range[0]=FrecRange
	303	303	return dataOut
	304	304
		305
		306	from scipy.optimize import fmin
		307	import itertools
		308	from scipy.optimize import curve_fit
		309
		310
	305	311	class GaussianFit(Operation):
	306	312
	307	313	'''
			@@ -321,135 +327,198 class GaussianFit(Operation):
	321	327	self.i=0
	322	328
	323	329
	324		def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
		330	# def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
		331	def run(self, dataOut, SNRdBlimit=-9, method='generalized'):
	325	332	"""This routine will find a couple of generalized Gaussians to a power spectrum
		333	methods: generalized, squared
	326	334	input: spc
	327	335	output:
	328		Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1~~,noise~~
		336	noise, amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1
	329	337	"""
	330
		338	print ('Entering ',method,' double Gaussian fit')
	331	339	self.spc = dataOut.data_pre[0].copy()
	332	340	self.Num_Hei = self.spc.shape[2]
	333	341	self.Num_Bin = self.spc.shape[1]
	334	342	self.Num_Chn = self.spc.shape[0]
	335		Vrange = dataOut.abscissaList
	336
	337		GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
	338		SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
	339		SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
	340		SPC_ch1[:] = numpy.NaN
	341		SPC_ch2[:] = numpy.NaN
	342
	343	343
	344	344	start_time = time.time()
	345	345
	346		noise_ = dataOut.spc_noise[0].copy()
	347
	348
	349	346	pool = Pool(processes=self.Num_Chn)
	350		args = [(~~Vrange~~, Ch, pnoise, ~~noise_~~, ~~num_intg~~, SNRlimit) for Ch in range(self.Num_Chn)]
		347	args = [(dataOut.spc_range[2], ich, dataOut.spc_noise[ich], dataOut.nIncohInt, SNRdBlimit) for ich in range(self.Num_Chn)]
	351	348	objs = [self for __ in range(self.Num_Chn)]
	352	349	attrs = list(zip(objs, args))
	353		~~gauSPC~~ = pool.map(target, attrs)
	354		dataOut.SPCparam = numpy.asarray(SPCparam)
	355
	356		''' Parameters:
	357		1. Amplitude
	358		2. Shift
	359		3. Width
	360		4. Power
	361		'''
		350	DGauFitParam = pool.map(target, attrs)
		351	# Parameters:
		352	# 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
		353	dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
		354
		355	# Double Gaussian Curves
		356	gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
		357	gau0[:] = numpy.NaN
		358	gau1 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
		359	gau1[:] = numpy.NaN
		360	x_mtr = numpy.transpose(numpy.tile(dataOut.getVelRange(1)[:-1], (self.Num_Hei,1)))
		361	for iCh in range(self.Num_Chn):
		362	N0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,0]] * self.Num_Bin))
		363	N1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,1]] * self.Num_Bin))
		364	A0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,0]] * self.Num_Bin))
		365	A1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,1]] * self.Num_Bin))
		366	v0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,0]] * self.Num_Bin))
		367	v1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,1]] * self.Num_Bin))
		368	s0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,0]] * self.Num_Bin))
		369	s1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,1]] * self.Num_Bin))
		370	if method == 'generalized':
		371	p0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,0]] * self.Num_Bin))
		372	p1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,1]] * self.Num_Bin))
		373	elif method == 'squared':
		374	p0 = 2.
		375	p1 = 2.
		376	gau0[iCh] = A0numpy.exp(-0.5numpy.abs((x_mtr-v0)/s0)**p0)+N0
		377	gau1[iCh] = A1numpy.exp(-0.5numpy.abs((x_mtr-v1)/s1)**p1)+N1
		378	dataOut.GaussFit0 = gau0
		379	dataOut.GaussFit1 = gau1
		380	print(numpy.shape(gau0))
		381	hei = 26
		382	print(dataOut.heightList[hei])
		383	#import matplotlib.pyplot as plt
		384	plt.plot(self.spc[0,:,hei])
		385	plt.plot(dataOut.GaussFit0[0,:,hei])
		386	plt.plot(dataOut.GaussFit1[0,:,hei])
		387	plt.plot(dataOut.GaussFit0[0,:,hei]+dataOut.GaussFit1[0,:,hei])
		388
		389	plt.show()
		390	time.sleep(60)
		391	#print(gau0)
		392
		393	print('Leaving ',method ,' double Gaussian fit')
		394	return dataOut
	362	395
	363	396	def FitGau(self, X):
	364
	365		Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
	366
	367		SPCparam = []
	368		SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
	369		SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
	370		SPC_ch1[:] = 0#numpy.NaN
	371		SPC_ch2[:] = 0#numpy.NaN
	372
	373
	374
		397	# print('Entering FitGau')
		398	# Assigning the variables
		399	Vrange, ch, wnoise, num_intg, SNRlimit = X
		400	# Noise Limits
		401	noisebl = wnoise * 0.9
		402	noisebh = wnoise * 1.1
		403	# Radar Velocity
		404	Va = max(Vrange)
		405	deltav = Vrange[1] - Vrange[0]
		406	x = numpy.arange(self.Num_Bin)
		407
		408	# print ('stop 0')
		409
		410	# 5 parameters, 2 Gaussians
		411	DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
		412	DGauFitParam[:] = numpy.NaN
		413
		414	# SPCparam = []
		415	# SPC_ch1 = numpy.zeros([self.Num_Bin,self.Num_Hei])
		416	# SPC_ch2 = numpy.zeros([self.Num_Bin,self.Num_Hei])
		417	# SPC_ch1[:] = 0 #numpy.NaN
		418	# SPC_ch2[:] = 0 #numpy.NaN
		419	# print ('stop 1')
	375	420	for ht in range(self.Num_Hei):
	376
	377
		421	# print (ht)
		422	# print ('stop 2')
		423	# Spectra at each range
	378	424	spc = numpy.asarray(self.spc)[ch,:,ht]
		425	snr = ( spc.mean() - wnoise ) / wnoise
		426	snrdB = 10.*numpy.log10(snr)
	379	427
		428	#print ('stop 3')
		429	if snrdB < SNRlimit :
		430	# snr = numpy.NaN
		431	# SPC_ch1[:,ht] = 0#numpy.NaN
		432	# SPC_ch1[:,ht] = 0#numpy.NaN
		433	# SPCparam = (SPC_ch1,SPC_ch2)
		434	# print ('SNR less than SNRth')
		435	continue
		436	# wnoise = hildebrand_sekhon(spc,num_intg)
		437	# print ('stop 2.01')
	380	438	#############################################
	381	439	# normalizing spc and noise
	382	440	# This part differs from gg1
	383		spc_norm_max = max(spc)
		441	# spc_norm_max = max(spc) #commented by D. Scipión 19.03.2021
	384	442	#spc = spc / spc_norm_max
	385		pnoise = pnoise #/ spc_norm_max
		443	# pnoise = pnoise #/ spc_norm_max #commented by D. Scipión 19.03.2021
	386	444	#############################################
	387	445
		446	# print ('stop 2.1')
	388	447	fatspectra=1.0
		448	# noise per channel.... we might want to use the noise at each range
	389	449
	390		wnoise = noise_ #/ spc_norm_max
		450	# wnoise = noise_ #/ spc_norm_max #commented by D. Scipión 19.03.2021
	391	451	#wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
	392	452	#if wnoise>1.1*pnoise: # to be tested later
	393	453	# wnoise=pnoise
	394		noisebl=wnoise*0.9;
	395		noisebh=wnoise*1.1
	396		spc=spc-wnoise
		454	# noisebl = wnoise*0.9
		455	# noisebh = wnoise*1.1
		456	spc = spc - wnoise # signal
	397	457
	398		minx=numpy.argmin(spc)
		458	# print ('stop 2.2')
		459	minx = numpy.argmin(spc)
	399	460	#spcs=spc.copy()
	400		spcs=numpy.roll(spc,-minx)
	401		cum=numpy.cumsum(spcs)
	402		tot_noise=wnoise * self.Num_Bin #64;
	403
	404		snr = sum(spcs)/tot_noise
	405		snrdB=10.*numpy.log10(snr)
	406
	407		if snrdB < SNRlimit :
	408		snr = numpy.NaN
	409		~~SPC_ch1[:,~~ht] = 0#numpy.NaN
	410		SPC_ch1[:,ht] = 0#numpy.NaN
	411		SPCparam = (SPC_ch1,SPC_ch2)
	412		continue
		461	spcs = numpy.roll(spc,-minx)
		462	cum = numpy.cumsum(spcs)
		463	# tot_noise = wnoise * self.Num_Bin #64;
		464
		465	# print ('stop 2.3')
		466	# snr = sum(spcs) / tot_noise
		467	# snrdB = 10.*numpy.log10(snr)
		468	#print ('stop 3')
		469	# if snrdB < SNRlimit :
		470	# snr = numpy.NaN
		471	# SPC_ch1[:,ht] = 0#numpy.NaN
		472	# SPC_ch1[:,ht] = 0#numpy.NaN
		473	# SPCparam = (SPC_ch1,SPC_ch2)
		474	# print ('SNR less than SNRth')
		475	# continue
	413	476
	414	477
	415	478	#if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
	416	479	# return [None,]4,[None,]4,None,snrdB,None,None,[None,]5,[None,]9,None
	417
	418		cummax=max(cum);
	419		epsi=0.08*fatspectra # cumsum to narrow down the energy region
	420		cumlo=cummax*epsi;
	421		cumhi=cummax*(1-epsi)
	422		powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
	423
	424
		480	# print ('stop 4')
		481	cummax = max(cum)
		482	epsi = 0.08 * fatspectra # cumsum to narrow down the energy region
		483	cumlo = cummax * epsi
		484	cumhi = cummax * (1-epsi)
		485	powerindex = numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
		486
		487	# print ('stop 5')
	425	488	if len(powerindex) < 1:# case for powerindex 0
		489	# print ('powerindex < 1')
		490	continue
		491	powerlo = powerindex[0]
		492	powerhi = powerindex[-1]
		493	powerwidth = powerhi-powerlo
		494	if powerwidth <= 1:
		495	# print('powerwidth <= 1')
	426	496	continue
	427		powerlo=powerindex[0]
	428		powerhi=powerindex[-1]
	429		powerwidth=powerhi-powerlo
	430	497
	431		firstpeak=powerlo+powerwidth/10.# first gaussian energy location
	432		~~secondpeak~~=~~powerhi~~-powerwidth/10.#~~second~~ gaussian energy location
	433		midpeak=(firstpeak+secondpeak)/2.
	434		firstamp=spcs[int(firstpeak)]
	435		~~secondamp~~=spcs[int(~~second~~peak)]
	436		~~midamp~~=spcs[int(midpeak)]
		498	# print ('stop 6')
		499	firstpeak = powerlo + powerwidth/10.# first gaussian energy location
		500	secondpeak = powerhi - powerwidth/10. #second gaussian energy location
		501	midpeak = (firstpeak + secondpeak)/2.
		502	firstamp = spcs[int(firstpeak)]
		503	secondamp = spcs[int(secondpeak)]
		504	midamp = spcs[int(midpeak)]
	437	505
	438		x=numpy.arange( self.Num_Bin )
	439		y_data=spc+wnoise
		506	y_data = spc + wnoise
	440	507
	441	508	''' single Gaussian '''
	442		shift0=numpy.mod(midpeak+minx, self.Num_Bin )
	443		width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
	444		power0=2.
	445		amplitude0=midamp
	446		state0=[shift0,width0,amplitude0,power0,wnoise]
	447		bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
	448		lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
	449
	450		chiSq1=lsq1[1];
	451
	452
		509	shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
		510	width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
		511	power0 = 2.
		512	amplitude0 = midamp
		513	state0 = [shift0,width0,amplitude0,power0,wnoise]
		514	bnds = ((0,self.Num_Bin-1),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
		515	lsq1 = fmin_l_bfgs_b(self.misfit1, state0, args=(y_data,x,num_intg), bounds=bnds, approx_grad=True)
		516	# print ('stop 7.1')
		517	# print (bnds)
		518
		519	chiSq1=lsq1[1]
		520
		521	# print ('stop 8')
	453	522	if fatspectra<1.0 and powerwidth<4:
	454	523	choice=0
	455	524	Amplitude0=lsq1[0][2]
			@@ -464,133 +533,378 class GaussianFit(Operation):
	464	533	#return (numpy.array([shift0,width0,Amplitude0,p0]),
	465	534	# numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
	466	535
	467		''' two gaussians '''
		536	# print ('stop 9')
		537	''' two Gaussians '''
	468	538	#shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
	469		shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
	470		shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
	471		width0=powerwidth/6.;
	472		width1=width0
	473		power0=2.;
	474		power1=power0
	475		amplitude0=firstamp;
	476		amplitude1=secondamp
	477		state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
		539	shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
		540	shift1 = numpy.mod(secondpeak+minx, self.Num_Bin )
		541	width0 = powerwidth/6.
		542	width1 = width0
		543	power0 = 2.
		544	power1 = power0
		545	amplitude0 = firstamp
		546	amplitude1 = secondamp
		547	state0 = [shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
	478	548	#bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
	479		bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
		549	bnds=((0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
	480	550	#bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
	481	551
		552	# print ('stop 10')
	482	553	lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
	483	554
		555	# print ('stop 11')
		556	chiSq2 = lsq2[1]
	484	557
	485		chiSq2=lsq2[1];
	486
	487
		558	# print ('stop 12')
	488	559
	489		oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
		560	oneG = (chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
	490	561
		562	# print ('stop 13')
	491	563	if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
	492	564	if oneG:
	493		choice=0
		565	choice = 0
	494	566	else:
	495		w1=lsq2[0][1]; w2=lsq2[0][5]
	496		a1=lsq2[0][2]; a2=lsq2[0][6]
	497		p1=lsq2[0][3]; p2=lsq2[0][7]
	498		s1=(2*(1+1./p1))scipy.special.gamma(1./p1)/p1;
	499		s2=(2*(1+1./p2))scipy.special.gamma(1./p2)/p2;
	500		gp1=a1w1s1; gp2=a2w2s2 # power content of each ggaussian with proper p scaling
		567	w1 = lsq2[0][1]; w2 = lsq2[0][5]
		568	a1 = lsq2[0][2]; a2 = lsq2[0][6]
		569	p1 = lsq2[0][3]; p2 = lsq2[0][7]
		570	s1 = (2*(1+1./p1))scipy.special.gamma(1./p1)/p1
		571	s2 = (2*(1+1./p2))scipy.special.gamma(1./p2)/p2
		572	gp1 = a1w1s1; gp2 = a2w2s2 # power content of each ggaussian with proper p scaling
	501	573
	502	574	if gp1>gp2:
	503	575	if a1>0.7*a2:
	504		choice=1
		576	choice = 1
	505	577	else:
	506		choice=2
		578	choice = 2
	507	579	elif gp2>gp1:
	508	580	if a2>0.7*a1:
	509		choice=2
		581	choice = 2
	510	582	else:
	511		choice=1
		583	choice = 1
	512	584	else:
	513		choice=numpy.argmax([a1,a2])+1
		585	choice = numpy.argmax([a1,a2])+1
	514	586	#else:
	515	587	#choice=argmin([std2a,std2b])+1
	516	588
	517	589	else: # with low SNR go to the most energetic peak
	518		choice=numpy.argmax([lsq1[0][2]lsq1[0][1],lsq2[0][2]lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
	519
	520
	521		shift0=lsq2[0][0];
	522		vel0=Vrange[0] + shift0*(~~Vrange~~[1]-~~Vrange~~[0])
	523		shift1=lsq2[0][4];
	524		vel1=Vrange[0] + shift1*(~~Vrange~~[1]-~~Vrange~~[0])
	525
	526		max_vel = 1.0
	527
		590	choice = numpy.argmax([lsq1[0][2]lsq1[0][1],lsq2[0][2]lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
		591
		592	# print ('stop 14')
		593	shift0 = lsq2[0][0]
		594	vel0 = Vrange[0] + shift0 * deltav
		595	shift1 = lsq2[0][4]
		596	# vel1=Vrange[0] + shift1 * deltav
		597
		598	# max_vel = 1.0
		599	# Va = max(Vrange)
		600	# deltav = Vrange[1]-Vrange[0]
		601	# print ('stop 15')
	528	602	#first peak will be 0, second peak will be 1
	529		if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
	530		shift0=lsq2[0][0]
	531		~~width0~~=lsq2[0][1]
	532		~~Amplitude0~~=lsq2[0][2]
	533		p0=lsq2[0][3]
	534
	535		shift1=lsq2[0][4]
	536		~~width1~~=lsq2[0][5]
	537		~~Amplitude1~~=lsq2[0][6]
	538		p1=lsq2[0][7]
	539		~~noise~~=lsq2[0][8]
		603	# if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.Scipión 19.03.2021
		604	if vel0 > -Va and vel0 < Va : #first peak is in the correct range
		605	shift0 = lsq2[0][0]
		606	width0 = lsq2[0][1]
		607	Amplitude0 = lsq2[0][2]
		608	p0 = lsq2[0][3]
		609
		610	shift1 = lsq2[0][4]
		611	width1 = lsq2[0][5]
		612	Amplitude1 = lsq2[0][6]
		613	p1 = lsq2[0][7]
		614	noise = lsq2[0][8]
	540	615	else:
	541		shift1=lsq2[0][0]
	542		width1=lsq2[0][1]
	543		Amplitude1=lsq2[0][2]
	544		p1=lsq2[0][3]
		616	shift1 = lsq2[0][0]
		617	width1 = lsq2[0][1]
		618	Amplitude1 = lsq2[0][2]
		619	p1 = lsq2[0][3]
	545	620
	546		shift0=lsq2[0][4]
	547		width0=lsq2[0][5]
	548		Amplitude0=lsq2[0][6]
	549		p0=lsq2[0][7]
	550		noise=lsq2[0][8]
		621	shift0 = lsq2[0][4]
		622	width0 = lsq2[0][5]
		623	Amplitude0 = lsq2[0][6]
		624	p0 = lsq2[0][7]
		625	noise = lsq2[0][8]
	551	626
	552	627	if Amplitude0<0.05: # in case the peak is noise
	553		shift0,width0,Amplitude0,p0 = [0,0,0,0]~~#4*[numpy.NaN]~~
		628	shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
	554	629	if Amplitude1<0.05:
	555		shift1,width1,Amplitude1,p1 = [0,0,0,0]~~#4*[numpy.NaN]~~
		630	shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
		631
		632	# print ('stop 16 ')
		633	# SPC_ch1[:,ht] = noise + Amplitude0numpy.exp(-0.5(abs(x-shift0)/width0)**p0)
		634	# SPC_ch2[:,ht] = noise + Amplitude1numpy.exp(-0.5(abs(x-shift1)/width1)**p1)
		635	# SPCparam = (SPC_ch1,SPC_ch2)
		636
		637	DGauFitParam[0,ht,0] = noise
		638	DGauFitParam[0,ht,1] = noise
		639	DGauFitParam[1,ht,0] = Amplitude0
		640	DGauFitParam[1,ht,1] = Amplitude1
		641	DGauFitParam[2,ht,0] = Vrange[0] + shift0 * deltav
		642	DGauFitParam[2,ht,1] = Vrange[0] + shift1 * deltav
		643	DGauFitParam[3,ht,0] = width0 * deltav
		644	DGauFitParam[3,ht,1] = width1 * deltav
		645	DGauFitParam[4,ht,0] = p0
		646	DGauFitParam[4,ht,1] = p1
		647
		648	# print (DGauFitParam.shape)
		649	# print ('Leaving FitGau')
		650	return DGauFitParam
		651	# return SPCparam
		652	# return GauSPC
	556	653
		654	def y_model1(self,x,state):
		655	shift0, width0, amplitude0, power0, noise = state
		656	model0 = amplitude0numpy.exp(-0.5abs((x - shift0)/width0)**power0)
		657	model0u = amplitude0numpy.exp(-0.5abs((x - shift0 - self.Num_Bin)/width0)**power0)
		658	model0d = amplitude0numpy.exp(-0.5abs((x - shift0 + self.Num_Bin)/width0)**power0)
		659	return model0 + model0u + model0d + noise
	557	660
	558		SPC_ch1[:,ht] = noise + Amplitude0numpy.exp(-0.5(abs(x-shift0))/width0)**p0
	559		SPC_ch2[:,ht] = noise + Amplitude1numpy.exp(-0.5(abs(x-shift1))/width1)**p1
	560		SPCparam = (SPC_ch1,SPC_ch2)
		661	def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
		662	shift0, width0, amplitude0, power0, shift1, width1, amplitude1, power1, noise = state
		663	model0 = amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)
		664	model0u = amplitude0numpy.exp(-0.5abs((x - shift0 - self.Num_Bin)/width0)**power0)
		665	model0d = amplitude0numpy.exp(-0.5abs((x - shift0 + self.Num_Bin)/width0)**power0)
	561	666
		667	model1 = amplitude1numpy.exp(-0.5abs((x - shift1)/width1)**power1)
		668	model1u = amplitude1numpy.exp(-0.5abs((x - shift1 - self.Num_Bin)/width1)**power1)
		669	model1d = amplitude1numpy.exp(-0.5abs((x - shift1 + self.Num_Bin)/width1)**power1)
		670	return model0 + model0u + model0d + model1 + model1u + model1d + noise
	562	671
	563		return GauSPC
		672	def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
	564	673
	565		def y_model1(self,x,state):
	566		shift0,width0,amplitude0,power0,noise=state
	567		model0=amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)
		674	return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))*2)#/(64-5.) # /(64-5.) can be commented
	568	675
	569		model0u=amplitude0numpy.exp(-0.5abs((x-shift0- self.Num_Bin )/width0)**power0)
		676	def misfit2(self,state,y_data,x,num_intg):
		677	return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))*2)#/(64-9.)
	570	678
	571		model0d=amplitude0numpy.exp(-0.5abs((x-shift0+ self.Num_Bin )/width0)**power0)
	572		return model0+model0u+model0d+noise
	573	679
	574		def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
	575		shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
	576		model0=amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)
		680	class Oblique_Gauss_Fit(Operation):
		681
		682	def __init__(self):
		683	Operation.__init__(self)
	577	684
	578		model0u=amplitude0numpy.exp(-0.5abs((x-shift0- self.Num_Bin )/width0)**power0)
	579	685
	580		model0d=amplitude0numpy.exp(-0.5abs((x-shift0+ self.Num_Bin )/width0)**power0)
	581		model1=amplitude1numpy.exp(-0.5abs((x-shift1)/width1)**power1)
	582	686
	583		model1u=amplitude1numpy.exp(-0.5abs((x-shift1- self.Num_Bin )/width1)**power1)
		687	def Gauss_fit(self,spc,x,nGauss):
	584	688
	585		model1d=amplitude1numpy.exp(-0.5abs((x-shift1+ self.Num_Bin )/width1)**power1)
	586		return model0+model0u+model0d+model1+model1u+model1d+noise
	587	689
	588		def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
		690	def gaussian(x, a, b, c, d):
		691	val = a * numpy.exp(-(x - b)*2 / (2c**2)) + d
		692	return val
	589	693
	590		return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))*2)#/(64-5.) # /(64-5.) can be commented
		694	if nGauss == 'first':
		695	spc_1_aux = numpy.copy(spc[:numpy.argmax(spc)+1])
		696	spc_2_aux = numpy.flip(spc_1_aux)
		697	spc_3_aux = numpy.concatenate((spc_1_aux,spc_2_aux[1:]))
		698
		699	len_dif = len(x)-len(spc_3_aux)
		700
		701	spc_zeros = numpy.ones(len_dif)*spc_1_aux[0]
		702
		703	spc_new = numpy.concatenate((spc_3_aux,spc_zeros))
		704
		705	y = spc_new
		706
		707	elif nGauss == 'second':
		708	y = spc
		709
		710
		711	# estimate starting values from the data
		712	a = y.max()
		713	b = x[numpy.argmax(y)]
		714	if nGauss == 'first':
		715	c = 1.#b#b#numpy.std(spc)
		716	elif nGauss == 'second':
		717	c = b
		718	else:
		719	print("ERROR")
		720
		721	d = numpy.mean(y[-100:])
		722
		723	# define a least squares function to optimize
		724	def minfunc(params):
		725	return sum((y-gaussian(x,params[0],params[1],params[2],params[3]))**2)
		726
		727	# fit
		728	popt = fmin(minfunc,[a,b,c,d],disp=False)
		729	#popt,fopt,niter,funcalls = fmin(minfunc,[a,b,c,d])
		730
		731
		732	return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
		733
		734
		735	def Gauss_fit_2(self,spc,x,nGauss):
		736
		737
		738	def gaussian(x, a, b, c, d):
		739	val = a * numpy.exp(-(x - b)*2 / (2c**2)) + d
		740	return val
		741
		742	if nGauss == 'first':
		743	spc_1_aux = numpy.copy(spc[:numpy.argmax(spc)+1])
		744	spc_2_aux = numpy.flip(spc_1_aux)
		745	spc_3_aux = numpy.concatenate((spc_1_aux,spc_2_aux[1:]))
		746
		747	len_dif = len(x)-len(spc_3_aux)
		748
		749	spc_zeros = numpy.ones(len_dif)*spc_1_aux[0]
		750
		751	spc_new = numpy.concatenate((spc_3_aux,spc_zeros))
		752
		753	y = spc_new
		754
		755	elif nGauss == 'second':
		756	y = spc
		757
		758
		759	# estimate starting values from the data
		760	a = y.max()
		761	b = x[numpy.argmax(y)]
		762	if nGauss == 'first':
		763	c = 1.#b#b#numpy.std(spc)
		764	elif nGauss == 'second':
		765	c = b
		766	else:
		767	print("ERROR")
		768
		769	d = numpy.mean(y[-100:])
		770
		771	# define a least squares function to optimize
		772	popt,pcov = curve_fit(gaussian,x,y,p0=[a,b,c,d])
		773	#popt,fopt,niter,funcalls = fmin(minfunc,[a,b,c,d])
		774
		775
		776	#return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
		777	return gaussian(x, popt[0], popt[1], popt[2], popt[3]),popt[0], popt[1], popt[2], popt[3]
		778
		779	def Double_Gauss_fit(self,spc,x,A1,B1,C1,A2,B2,C2,D):
		780
		781	def double_gaussian(x, a1, b1, c1, a2, b2, c2, d):
		782	val = a1 * numpy.exp(-(x - b1)*2 / (2c1*2)) + a2 numpy.exp(-(x - b2)*2 / (2c2**2)) + d
		783	return val
		784
		785
		786	y = spc
		787
		788	# estimate starting values from the data
		789	a1 = A1
		790	b1 = B1
		791	c1 = C1#numpy.std(spc)
		792
		793	a2 = A2#y.max()
		794	b2 = B2#x[numpy.argmax(y)]
		795	c2 = C2#numpy.std(spc)
		796	d = D
		797
		798	# define a least squares function to optimize
		799	def minfunc(params):
		800	return sum((y-double_gaussian(x,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))**2)
		801
		802	# fit
		803	popt = fmin(minfunc,[a1,b1,c1,a2,b2,c2,d],disp=False)
		804
		805	return double_gaussian(x, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
		806
		807	def Double_Gauss_fit_2(self,spc,x,A1,B1,C1,A2,B2,C2,D):
		808
		809	def double_gaussian(x, a1, b1, c1, a2, b2, c2, d):
		810	val = a1 * numpy.exp(-(x - b1)*2 / (2c1*2)) + a2 numpy.exp(-(x - b2)*2 / (2c2**2)) + d
		811	return val
		812
		813
		814	y = spc
		815
		816	# estimate starting values from the data
		817	a1 = A1
		818	b1 = B1
		819	c1 = C1#numpy.std(spc)
		820
		821	a2 = A2#y.max()
		822	b2 = B2#x[numpy.argmax(y)]
		823	c2 = C2#numpy.std(spc)
		824	d = D
		825
		826	# fit
		827
		828	popt,pcov = curve_fit(double_gaussian,x,y,p0=[a1,b1,c1,a2,b2,c2,d])
		829
		830	error = numpy.sqrt(numpy.diag(pcov))
		831
		832	return popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], error[0], error[1], error[2], error[3], error[4], error[5], error[6]
		833
		834
		835
		836
		837	def run(self, dataOut):
		838
		839	pwcode = 1
		840
		841	if dataOut.flagDecodeData:
		842	pwcode = numpy.sum(dataOut.code[0]**2)
		843	#normFactor = min(self.nFFTPoints,self.nProfiles)self.nIncohIntself.nCohIntpwcodeself.windowOfFilter
		844	normFactor = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * pwcode * dataOut.windowOfFilter
		845	factor = normFactor
		846	z = dataOut.data_spc / factor
		847	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
		848	dataOut.power = numpy.average(z, axis=1)
		849	dataOut.powerdB = 10 * numpy.log10(dataOut.power)
		850
		851
		852	x = dataOut.getVelRange(0)
		853	#print(aux)
		854	#print(numpy.shape(aux))
		855	#exit(1)
		856
		857	#print(numpy.shape(dataOut.data_spc))
		858
		859	dataOut.Oblique_params = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
		860	dataOut.Oblique_param_errors = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
		861
		862	dataOut.VelRange = x
		863
		864
		865	l1=range(22,36)
		866	l2=range(58,99)
		867
		868	for hei in itertools.chain(l1, l2):
		869	#print("INIT")
		870	#print(hei)
		871
		872	try:
		873	spc = dataOut.data_spc[0,:,hei]
		874
		875	spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
		876
		877	spc_diff = spc - spc_fit
		878	spc_diff[spc_diff < 0] = 0
		879
		880	spc_fit_diff, A2, B2, C2, D2 = self.Gauss_fit_2(spc_diff,x,'second')
		881
		882	D = (D1+D2)
		883
		884	dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_param_errors[0,0,hei],dataOut.Oblique_param_errors[0,1,hei],dataOut.Oblique_param_errors[0,2,hei],dataOut.Oblique_param_errors[0,3,hei],dataOut.Oblique_param_errors[0,4,hei],dataOut.Oblique_param_errors[0,5,hei],dataOut.Oblique_param_errors[0,6,hei] = self.Double_Gauss_fit_2(spc,x,A1,B1,C1,A2,B2,C2,D)
		885	#spc_double_fit,dataOut.Oblique_params = self.Double_Gauss_fit(spc,x,A1,B1,C1,A2,B2,C2,D)
		886	#print(dataOut.Oblique_params)
		887	except:
		888	###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
		889	pass
		890	#print("DONE")
		891	'''
		892	print(dataOut.Oblique_params[1])
		893	print(dataOut.Oblique_params[4])
		894	import matplotlib.pyplot as plt
		895	plt.plot(x,spc_double_fit)
		896	plt.show()
		897	import time
		898	time.sleep(5)
		899	plt.close()
		900	'''
		901
		902
		903
		904
		905
		906	return dataOut
	591	907
	592		def misfit2(self,state,y_data,x,num_intg):
	593		return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))*2)#/(64-9.)
	594	908
	595	909
	596	910
			@@ -3998,3 +4312,55 class SMOperations():
	3998	4312	# error[indInvalid1] = 13
	3999	4313	#
	4000	4314	# return heights, error
		4315
		4316
		4317
		4318	class IGRFModel(Operation):
		4319	"""Operation to calculate Geomagnetic parameters.
		4320
		4321	Parameters:
		4322	-----------
		4323	None
		4324
		4325	Example
		4326	--------
		4327
		4328	op = proc_unit.addOperation(name='IGRFModel', optype='other')
		4329
		4330	"""
		4331
		4332	def __init__(self, **kwargs):
		4333
		4334	Operation.__init__(self, **kwargs)
		4335
		4336	self.aux=1
		4337
		4338	def run(self,dataOut):
		4339
		4340	try:
		4341	from schainpy.model.proc import mkfact_short_2020
		4342	except:
		4343	log.warning('You should install "mkfact_short_2020" module to process IGRF Model')
		4344
		4345	if self.aux==1:
		4346
		4347	#dataOut.TimeBlockSeconds_First_Time=time.mktime(time.strptime(dataOut.TimeBlockDate))
		4348	#### we do not use dataOut.datatime.ctime() because it's the time of the second (next) block
		4349	dataOut.TimeBlockSeconds_First_Time=dataOut.TimeBlockSeconds
		4350	dataOut.bd_time=time.gmtime(dataOut.TimeBlockSeconds_First_Time)
		4351	dataOut.year=dataOut.bd_time.tm_year+(dataOut.bd_time.tm_yday-1)/364.0
		4352	dataOut.ut=dataOut.bd_time.tm_hour+dataOut.bd_time.tm_min/60.0+dataOut.bd_time.tm_sec/3600.0
		4353
		4354	self.aux=0
		4355
		4356	dataOut.h=numpy.arange(0.0,15.0*dataOut.MAXNRANGENDT,15.0,dtype='float32')
		4357	dataOut.bfm=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
		4358	dataOut.bfm=numpy.array(dataOut.bfm,order='F')
		4359	dataOut.thb=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
		4360	dataOut.thb=numpy.array(dataOut.thb,order='F')
		4361	dataOut.bki=numpy.zeros(dataOut.MAXNRANGENDT,dtype='float32')
		4362	dataOut.bki=numpy.array(dataOut.bki,order='F')
		4363
		4364	mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
		4365
		4366	return dataOut

schainpy/model/proc/jroproc_spectra.py +1 -1

                          dataOut.utctime = avgdatatime
                          dataOut.flagNoData = False
-                     return dataOut
  No newline at end of file
+                     return dataOut

schainpy/model/proc/jroproc_spectra_lags.py +1 -1

                      self.dataOut.noise_estimation = noise.copy()
-                     return 1
  No newline at end of file
+                     return 1

schainpy/model/proc/jroproc_voltage.py 0 0

	1		NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff

schainpy/model/proc/pxproc_parameters.py 0 0

	1		NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff

schainpy/scripts/PPD.py 0 0

	1		NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages