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			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	sortdata = numpy.sort(data, axis=None)	78	sortdata = numpy.sort(data, axis=None)
		79	#print(numpy.shape(data))
		80	#exit()
79	'''	81	'''
80	lenOfData = len(sortdata)	82	lenOfData = len(sortdata)
81	nums_min = lenOfData*0.2	83	nums_min = lenOfData*0.2
@@ -273,13 +275,13 class JROData(GenericData):
273	'''	275	'''
274	'''	276	'''
275	return self.radarControllerHeaderObj.ippSeconds	277	return self.radarControllerHeaderObj.ippSeconds
276		278
277	@ippSeconds.setter	279	@ippSeconds.setter
278	def ippSeconds(self, ippSeconds):	280	def ippSeconds(self, ippSeconds):
279	'''	281	'''
280	'''	282	'''
281	self.radarControllerHeaderObj.ippSeconds = ippSeconds	283	self.radarControllerHeaderObj.ippSeconds = ippSeconds
282		284
283	@property	285	@property
284	def code(self):	286	def code(self):
285	'''	287	'''
@@ -370,7 +372,7 class Voltage(JROData):
370	self.flagShiftFFT = False	372	self.flagShiftFFT = False
371	self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil	373	self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372	self.profileIndex = 0	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	'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']	376	'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375		377
376	def getNoisebyHildebrand(self, channel=None):	378	def getNoisebyHildebrand(self, channel=None):
@@ -428,6 +430,103 class Voltage(JROData):
428	noise = property(getNoise, "I'm the 'nHeights' property.")	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	class Spectra(JROData):	530	class Spectra(JROData):
432		531
433	def __init__(self):	532	def __init__(self):
@@ -458,7 +557,7 class Spectra(JROData):
458	self.ippFactor = 1	557	self.ippFactor = 1
459	self.beacon_heiIndexList = []	558	self.beacon_heiIndexList = []
460	self.noise_estimation = None	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	'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']	561	'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
463		562
464	def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):	563	def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
@@ -608,7 +707,7 class Spectra(JROData):
608	print("This property should not be initialized")	707	print("This property should not be initialized")
609		708
610	return	709	return
611		710
612	noise = property(getNoise, setValue, "I'm the 'nHeights' property.")	711	noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
613		712
614		713
@@ -705,7 +804,7 class Fits(JROData):
705	return self.ipp_sec	804	return self.ipp_sec
706		805
707	noise = property(getNoise, "I'm the 'nHeights' property.")	806	noise = property(getNoise, "I'm the 'nHeights' property.")
708		807
709		808
710	class Correlation(JROData):	809	class Correlation(JROData):
711		810
@@ -886,6 +985,7 class Parameters(Spectra):
886	else:	985	else:
887	return self.paramInterval	986	return self.paramInterval
888		987
		988
889	def setValue(self, value):	989	def setValue(self, value):
890		990
891	print("This property should not be initialized")	991	print("This property should not be initialized")
@@ -928,6 +1028,10 class PlotterData(object):
928	self.plottypes = ['cspc', 'spc', 'noise', 'rti']	1028	self.plottypes = ['cspc', 'spc', 'noise', 'rti']
929	elif code == 'rti':	1029	elif code == 'rti':
930	self.plottypes = ['noise', 'rti']	1030	self.plottypes = ['noise', 'rti']
		1031	elif code == 'crossprod':
		1032	self.plottypes = ['crossprod', 'kay']
		1033	elif code == 'spectrogram':
		1034	self.plottypes = ['spc', 'spectrogram']
931	else:	1035	else:
932	self.plottypes = [code]	1036	self.plottypes = [code]
933		1037
@@ -976,9 +1080,11 class PlotterData(object):
976	plot = 'snr'	1080	plot = 'snr'
977	elif 'spc_moments' == plot:	1081	elif 'spc_moments' == plot:
978	plot = 'moments'	1082	plot = 'moments'
		1083	elif 'spc_oblique' == plot:
		1084	plot = 'oblique'
979	self.data[plot] = {}	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	self.data['noise'] = {}	1088	self.data['noise'] = {}
983	self.data['rti'] = {}	1089	self.data['rti'] = {}
984	if 'noise' not in self.plottypes:	1090	if 'noise' not in self.plottypes:
@@ -1020,16 +1126,33 class PlotterData(object):
1020	self.__heights.append(dataOut.heightList)	1126	self.__heights.append(dataOut.heightList)
1021	self.__all_heights.update(dataOut.heightList)	1127	self.__all_heights.update(dataOut.heightList)
1022		1128
		1129
		1130
1023	for plot in self.plottypes:	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	z = dataOut.data_spc/dataOut.normFactor	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	buffer = 10*numpy.log10(z)	1146	buffer = 10*numpy.log10(z)
1027	if plot == 'cspc':	1147	if plot == 'cspc':
1028	buffer = (dataOut.data_spc, dataOut.data_cspc)	1148	buffer = (dataOut.data_spc, dataOut.data_cspc)
		1149	self.nFactor=dataOut.normFactor
1029	if plot == 'noise':	1150	if plot == 'noise':
1030	buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)	1151	buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1031	if plot in ('rti', 'spcprofile'):	1152	if plot in ('rti', 'spcprofile'):
1032	buffer = dataOut.getPower()	1153	buffer = dataOut.getPower()
		1154	#print(buffer[0,0])
		1155	#exit()
1033	if plot == 'snr_db':	1156	if plot == 'snr_db':
1034	buffer = dataOut.data_SNR	1157	buffer = dataOut.data_SNR
1035	if plot == 'snr':	1158	if plot == 'snr':
@@ -1048,6 +1171,277 class PlotterData(object):
1048	buffer = dataOut.data_output	1171	buffer = dataOut.data_output
1049	if plot == 'param':	1172	if plot == 'param':
1050	buffer = dataOut.data_param	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	if plot == 'scope':	1445	if plot == 'scope':
1052	buffer = dataOut.data	1446	buffer = dataOut.data
1053	self.flagDataAsBlock = dataOut.flagDataAsBlock	1447	self.flagDataAsBlock = dataOut.flagDataAsBlock
@@ -1076,6 +1470,10 class PlotterData(object):
1076	elif plot == 'spc_moments':	1470	elif plot == 'spc_moments':
1077	self.data['spc'][tm] = buffer	1471	self.data['spc'][tm] = buffer
1078	self.data['moments'][tm] = dataOut.moments	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	else:	1477	else:
1080	if self.buffering:	1478	if self.buffering:
1081	self.data[plot][tm] = buffer	1479	self.data[plot][tm] = buffer
@@ -1141,6 +1539,7 class PlotterData(object):
1141	meta['interval'] = float(self.interval)	1539	meta['interval'] = float(self.interval)
1142	meta['localtime'] = self.localtime	1540	meta['localtime'] = self.localtime
1143	meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())	1541	meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
		1542
1144	if 'spc' in self.data or 'cspc' in self.data:	1543	if 'spc' in self.data or 'cspc' in self.data:
1145	meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())	1544	meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1146	else:	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'
                     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(
                             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
                     '''
                     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:
                         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
                             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	plot_type = 'pcolor'	22	plot_type = 'pcolor'
	23		23
	24	def setup(self):	24	def setup(self):
			25
	25	self.nplots = len(self.data.channels)	26	self.nplots = len(self.data.channels)
	26	self.ncols = int(numpy.sqrt(self.nplots) + 0.9)	27	self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
	27	self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)	28	self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
	@@ -31,10 +32,13 class SpectraPlot(Plot):
	31	self.width = 4 * self.ncols	32	self.width = 4 * self.ncols
	32	else:	33	else:
	33	self.width = 3.5 * self.ncols	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	self.ylabel = 'Range [km]'	36	self.ylabel = 'Range [km]'
	36		37
	37	def plot(self):	38	def plot(self):
			39
			40	#print(self.xaxis)
			41	#exit(1)
	38	if self.xaxis == "frequency":	42	if self.xaxis == "frequency":
	39	x = self.data.xrange[0]	43	x = self.data.xrange[0]
	40	self.xlabel = "Frequency (kHz)"	44	self.xlabel = "Frequency (kHz)"
	@@ -51,19 +55,25 class SpectraPlot(Plot):
	51		55
	52	self.titles = []	56	self.titles = []
	53		57
			58
	54	y = self.data.heights	59	y = self.data.heights
	55	self.y = y	60	self.y = y
	56	z = self.data['spc']	61	z = self.data['spc']
	57		62
			63	self.CODE2 = 'spc_oblique'
			64
			65
	58	for n, ax in enumerate(self.axes):	66	for n, ax in enumerate(self.axes):
	59	noise = self.data['noise'][n][-1]	67	noise = self.data['noise'][n][-1]
	60	if self.CODE == 'spc_moments':	68	if self.CODE == 'spc_moments':
	61	mean = self.data['moments'][n, :, 1, :][-1]	69	mean = self.data['moments'][n, :, 1, :][-1]
			70
	62	if ax.firsttime:	71	if ax.firsttime:
	63	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)	72	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
	64	self.xmin = self.xmin if self.xmin else -self.xmax	73	self.xmin = self.xmin if self.xmin else -self.xmax
	65	self.zmin = self.zmin if self.zmin else numpy.nanmin(z)	74	self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
	66	self.zmax = self.zmax if self.zmax else numpy.nanmax(z)	75	self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
			76	#print(numpy.shape(x))
	67	ax.plt = ax.pcolormesh(x, y, z[n].T,	77	ax.plt = ax.pcolormesh(x, y, z[n].T,
	68	vmin=self.zmin,	78	vmin=self.zmin,
	69	vmax=self.zmax,	79	vmax=self.zmax,
	@@ -77,15 +87,122 class SpectraPlot(Plot):
	77	color="k", linestyle="dashed", lw=1)[0]	87	color="k", linestyle="dashed", lw=1)[0]
	78	if self.CODE == 'spc_moments':	88	if self.CODE == 'spc_moments':
	79	ax.plt_mean = ax.plot(mean, y, color='k')[0]	89	ax.plt_mean = ax.plot(mean, y, color='k')[0]
			90
	80	else:	91	else:
			92
	81	ax.plt.set_array(z[n].T.ravel())	93	ax.plt.set_array(z[n].T.ravel())
	82	if self.showprofile:	94	if self.showprofile:
	83	ax.plt_profile.set_data(self.data['rti'][n][-1], y)	95	ax.plt_profile.set_data(self.data['rti'][n][-1], y)
	84	ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)	96	ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
	85	if self.CODE == 'spc_moments':	97	if self.CODE == 'spc_moments':
	86	ax.plt_mean.set_data(mean, y)	98	ax.plt_mean.set_data(mean, y)
			99
	87	self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))	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	class CrossSpectraPlot(Plot):	207	class CrossSpectraPlot(Plot):
	91		208
	@@ -103,7 +220,7 class CrossSpectraPlot(Plot):
	103	self.nrows = len(self.data.pairs)	220	self.nrows = len(self.data.pairs)
	104	self.nplots = self.nrows * 4	221	self.nplots = self.nrows * 4
	105	self.width = 3.1 * self.ncols	222	self.width = 3.1 * self.ncols
	106	self.height = ~~2.6~~ * self.nrows	223	self.height = 5 * self.nrows
	107	self.ylabel = 'Range [km]'	224	self.ylabel = 'Range [km]'
	108	self.showprofile = False	225	self.showprofile = False
	109	self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})	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	else:	236	else:
	120	x = self.data.xrange[2]	237	x = self.data.xrange[2]
	121	self.xlabel = "Velocity (m/s)"	238	self.xlabel = "Velocity (m/s)"
	122		239
	123	self.titles = []	240	self.titles = []
	124		241
			242
	125	y = self.data.heights	243	y = self.data.heights
	126	self.y = y	244	self.y = y
	127	nspc = self.data['spc']	245	nspc = self.data['spc']
			246	#print(numpy.shape(self.data['spc']))
	128	spc = self.data['cspc'][0]	247	spc = self.data['cspc'][0]
			248	#print(numpy.shape(spc))
			249	#exit()
	129	cspc = self.data['cspc'][1]	250	cspc = self.data['cspc'][1]
			251	#print(numpy.shape(cspc))
			252	#exit()
	130		253
	131	for n in range(self.nrows):	254	for n in range(self.nrows):
	132	noise = self.data['noise'][:,-1]	255	noise = self.data['noise'][:,-1]
	133	pair = self.data.pairs[n]	256	pair = self.data.pairs[n]
			257	#print(pair)
			258	#exit()
	134	ax = self.axes[4 * n]	259	ax = self.axes[4 * n]
	135	if ax.firsttime:	260	if ax.firsttime:
	136	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)	261	self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
	137	self.xmin = self.xmin if self.xmin else -self.xmax	262	#self.xmin = self.xmin if self.xmin else -self.xmax
	138	self.zmin = self.zmin if self.zmin else numpy.nanmin(~~nspc~~)	263	self.xmin = self.xmin if self.xmin else numpy.nanmin(x)
	139	self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)	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	ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,	271	ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
	141	vmin=self.zmin,	272	vmin=self.zmin,
	142	vmax=self.zmax,	273	vmax=self.zmax,
	143	cmap=plt.get_cmap(self.colormap)	274	cmap=plt.get_cmap(self.colormap)
	144	)	275	)
	145	else:	276	else:
	146	ax.plt.set_array(nspc[pair[0]].T.ravel())	277	ax.plt.set_array(nspc[pair[0]].T.ravel())
	147	self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))	278	self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
	148		279
	@@ -153,10 +284,10 class CrossSpectraPlot(Plot):
	153	vmax=self.zmax,	284	vmax=self.zmax,
	154	cmap=plt.get_cmap(self.colormap)	285	cmap=plt.get_cmap(self.colormap)
	155	)	286	)
	156	else:	287	else:
	157	ax.plt.set_array(nspc[pair[1]].T.ravel())	288	ax.plt.set_array(nspc[pair[1]].T.ravel())
	158	self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))	289	self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
	159		290
	160	out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])	291	out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
	161	coh = numpy.abs(out)	292	coh = numpy.abs(out)
	162	phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi	293	phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
	@@ -172,19 +303,345 class CrossSpectraPlot(Plot):
	172	ax.plt.set_array(coh.T.ravel())	303	ax.plt.set_array(coh.T.ravel())
	173	self.titles.append(	304	self.titles.append(
	174	'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))	305	'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
	175		306
	176	ax = self.axes[4 * n + 3]	307	ax = self.axes[4 * n + 3]
	177	if ax.firsttime:	308	if ax.firsttime:
	178	ax.plt = ax.pcolormesh(x, y, phase.T,	309	ax.plt = ax.pcolormesh(x, y, phase.T,
	179	vmin=-180,	310	vmin=-180,
	180	vmax=180,	311	vmax=180,
	181	cmap=plt.get_cmap(self.colormap_phase)	312	cmap=plt.get_cmap(self.colormap_phase)
	182	)	313	)
	183	else:	314	else:
	184	ax.plt.set_array(phase.T.ravel())	315	ax.plt.set_array(phase.T.ravel())
	185	self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))	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	class RTIPlot(Plot):	645	class RTIPlot(Plot):
	189	'''	646	'''
	190	Plot for RTI data	647	Plot for RTI data
	@@ -202,7 +659,7 class RTIPlot(Plot):
	202	self.ylabel = 'Range [km]'	659	self.ylabel = 'Range [km]'
	203	self.xlabel = 'Time'	660	self.xlabel = 'Time'
	204	self.cb_label = 'dB'	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	self.titles = ['{} Channel {}'.format(	663	self.titles = ['{} Channel {}'.format(
	207	self.CODE.upper(), x) for x in range(self.nrows)]	664	self.CODE.upper(), x) for x in range(self.nrows)]
	208		665
	@@ -210,6 +667,78 class RTIPlot(Plot):
	210	self.x = self.data.times	667	self.x = self.data.times
	211	self.y = self.data.heights	668	self.y = self.data.heights
	212	self.z = self.data[self.CODE]	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	self.z = numpy.ma.masked_invalid(self.z)	742	self.z = numpy.ma.masked_invalid(self.z)
	214		743
	215	if self.decimation is None:	744	if self.decimation is None:
	@@ -280,7 +809,7 class PhasePlot(CoherencePlot):
	280		809
	281	class NoisePlot(Plot):	810	class NoisePlot(Plot):
	282	'''	811	'''
	283	Plot for noise	812	Plot for noise
	284	'''	813	'''
	285		814
	286	CODE = 'noise'	815	CODE = 'noise'
	@@ -296,6 +825,7 class NoisePlot(Plot):
	296	self.xlabel = 'Time'	825	self.xlabel = 'Time'
	297	self.titles = ['Noise']	826	self.titles = ['Noise']
	298	self.colorbar = False	827	self.colorbar = False
			828	self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.17, 'right':0.95})
	299		829
	300	def plot(self):	830	def plot(self):
	301		831
	@@ -315,7 +845,7 class NoisePlot(Plot):
	315	self.axes[0].lines[ch].set_data(x, y)	845	self.axes[0].lines[ch].set_data(x, y)
	316		846
	317	self.ymin = numpy.nanmin(Y) - 5	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	class PowerProfilePlot(Plot):	851	class PowerProfilePlot(Plot):
	@@ -342,10 +872,10 class PowerProfilePlot(Plot):
	342	self.y = y	872	self.y = y
	343		873
	344	x = self.data['spcprofile']	874	x = self.data['spcprofile']
	345		875
	346	if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9	876	if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
	347	if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1	877	if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
	348		878
	349	if self.axes[0].firsttime:	879	if self.axes[0].firsttime:
	350	for ch in self.data.channels:	880	for ch in self.data.channels:
	351	self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))	881	self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
	@@ -498,7 +1028,7 class BeaconPhase(Plot):
	498	server=None, folder=None, username=None, password=None,	1028	server=None, folder=None, username=None, password=None,
	499	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):	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	return dataOut	1032	return dataOut
	503		1033
	504	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):	1034	if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
	@@ -638,4 +1168,4 class BeaconPhase(Plot):
	638	thisDatetime=thisDatetime,	1168	thisDatetime=thisDatetime,
	639	update_figfile=update_figfile)	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

@@ -21,4 +21,9 from .jroIO_mira35c import *
21	from .julIO_param import *	21	from .julIO_param import *
22		22
23	from .pxIO_param import *	23	from .pxIO_param import *
24	from .jroIO_simulator import * No newline at end of file	24	from .jroIO_simulator import *
		25
		26	############DP############
		27	from .jroIO_dat import *
		28
		29	############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:
                     year = int(file[1:5])
                     doy = int(file[5:8])
                     set = int(file[8:11])
             def getDateFromRadarFile(file):
                 try:
                     year = int(file[1:5])
                     doy = int(file[5:8])
                     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
                 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(',')
                                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()
                             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,
                                expLabel='', last=False):
                     for path in folders:
                         files = glob.glob1(path, '*{}'.format(ext))
                         files.sort()
                         if last:
                             if files:
                                 fo = files[-1]
                                 try:
                                     dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
                                     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()
+                                #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
                 def searchFilesOffLine(self, path, startDate, endDate,
                                        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)
                 def searchFilesOnLine(self, path, startDate, endDate,
                                       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()
+                            #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 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
                         if fullfilename is not None:
                             break
-                        self.nTries = 1
+                        #self.nTries = 1
                         nextFile = True
                         if nFiles == (self.nFiles - 1):
                             log.log('Trying with next day...', self.name)
                             nextDay = True
                             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:
                         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
+                    #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
                 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:
                         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
                     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.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.name)
                                 time.sleep(self.delay)
                             if not(fullpath):
                                 raise schainpy.admin.SchainError(
                                     '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
                         else:
                             log.log("Searching files in {}".format(self.path), self.name)
                             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   :
                         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
                     utc = numpy.floor(self.dataOut.utctime)
                     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)
                     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.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.counter_records = 0
                     self.nrecords = None
                     self.flagNoMoreFiles = 0
                     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.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.name)
                             time.sleep(self.delay)
                         if not(fullpath):
                             raise schainpy.admin.SchainError(
                                 '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.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
                     self.setNextFile()
                 def readFirstHeader(self):
                     self.parseHeader()
                     self.parseData()
                     self.blockIndex = 0
                     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.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)
                             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.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()):
                         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
                         if self.ext == '.txt':
                             x = self.parameters.index(param.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:
                                 setattr(self.dataOut, value, dummy.reshape(1,-1))
                         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.paramInterval = min(self.intervals)
                     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()
                     return 1
                 def setup(self, path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs):
                     '''
                     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)
                     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]
                     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) :
                         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)
+                            #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
                     )
+                    #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 key in out:
                             rec.set2D(key, nrec, out[key][n])
                         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.createCatalog(**self.catalog)
                     header.createHeader(**self.header)
                     header.write()
                 def putData(self):
                     if self.dataOut.flagNoData:
                         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
                 def close(self):
                     if self.counter > 0:
-                        self.setHeader()
  No newline at end of file
+                        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.
                 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)
                     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.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.name)
                             time.sleep(self.delay)
                         if not(fullpath):
                             raise schainpy.admin.SchainError(
                                 '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
                     else:
                         log.log("Searching files in {}".format(self.path), self.name)
                         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.__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
                     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
                 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,
                                    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.blocksPerFile, ),
                                 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.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.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:
                             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

@@ -14,3 +14,9 from .jroproc_spectra_lags import *
14	from .jroproc_spectra_acf import *	14	from .jroproc_spectra_acf import *
15	from .bltrproc_parameters import *	15	from .bltrproc_parameters import *
16	from .pxproc_parameters import *	16	from .pxproc_parameters import *
		17
		18
		19	###########DP###########
		20	from .jroproc_voltage_lags import *
		21	###########DP###########
		22	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	dataOut.spcparam_range[0]=FrecRange	302	dataOut.spcparam_range[0]=FrecRange
	303	return dataOut	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	class GaussianFit(Operation):	311	class GaussianFit(Operation):
	306		312
	307	'''	313	'''
	@@ -321,135 +327,198 class GaussianFit(Operation):
	321	self.i=0	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	"""This routine will find a couple of generalized Gaussians to a power spectrum	332	"""This routine will find a couple of generalized Gaussians to a power spectrum
			333	methods: generalized, squared
	326	input: spc	334	input: spc
	327	output:	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	self.spc = dataOut.data_pre[0].copy()	339	self.spc = dataOut.data_pre[0].copy()
	332	self.Num_Hei = self.spc.shape[2]	340	self.Num_Hei = self.spc.shape[2]
	333	self.Num_Bin = self.spc.shape[1]	341	self.Num_Bin = self.spc.shape[1]
	334	self.Num_Chn = self.spc.shape[0]	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	start_time = time.time()	344	start_time = time.time()
	345		345
	346	noise_ = dataOut.spc_noise[0].copy()
	347
	348
	349	pool = Pool(processes=self.Num_Chn)	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	objs = [self for __ in range(self.Num_Chn)]	348	objs = [self for __ in range(self.Num_Chn)]
	352	attrs = list(zip(objs, args))	349	attrs = list(zip(objs, args))
	353	~~gauSPC~~ = pool.map(target, attrs)	350	DGauFitParam = pool.map(target, attrs)
	354	dataOut.SPCparam = numpy.asarray(SPCparam)	351	# Parameters:
	355		352	# 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
	356	''' Parameters:	353	dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
	357	1. Amplitude	354
	358	2. Shift	355	# Double Gaussian Curves
	359	3. Width	356	gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
	360	4. Power	357	gau0[:] = numpy.NaN
	361	'''	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	def FitGau(self, X):	396	def FitGau(self, X):
	364		397	# print('Entering FitGau')
	365	Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X	398	# Assigning the variables
	366		399	Vrange, ch, wnoise, num_intg, SNRlimit = X
	367	SPCparam = []	400	# Noise Limits
	368	SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])	401	noisebl = wnoise * 0.9
	369	SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])	402	noisebh = wnoise * 1.1
	370	SPC_ch1[:] = 0#numpy.NaN	403	# Radar Velocity
	371	SPC_ch2[:] = 0#numpy.NaN	404	Va = max(Vrange)
	372		405	deltav = Vrange[1] - Vrange[0]
	373		406	x = numpy.arange(self.Num_Bin)
	374		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	for ht in range(self.Num_Hei):	420	for ht in range(self.Num_Hei):
	376		421	# print (ht)
	377		422	# print ('stop 2')
			423	# Spectra at each range
	378	spc = numpy.asarray(self.spc)[ch,:,ht]	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	# normalizing spc and noise	439	# normalizing spc and noise
	382	# This part differs from gg1	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	#spc = spc / spc_norm_max	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	fatspectra=1.0	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	#wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used	451	#wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
	392	#if wnoise>1.1*pnoise: # to be tested later	452	#if wnoise>1.1*pnoise: # to be tested later
	393	# wnoise=pnoise	453	# wnoise=pnoise
	394	noisebl=wnoise*0.9;	454	# noisebl = wnoise*0.9
	395	noisebh=wnoise*1.1	455	# noisebh = wnoise*1.1
	396	spc=spc-wnoise	456	spc = spc - wnoise # signal
	397		457
	398	minx=numpy.argmin(spc)	458	# print ('stop 2.2')
			459	minx = numpy.argmin(spc)
	399	#spcs=spc.copy()	460	#spcs=spc.copy()
	400	spcs=numpy.roll(spc,-minx)	461	spcs = numpy.roll(spc,-minx)
	401	cum=numpy.cumsum(spcs)	462	cum = numpy.cumsum(spcs)
	402	tot_noise=wnoise * self.Num_Bin #64;	463	# tot_noise = wnoise * self.Num_Bin #64;
	403		464
	404	snr = sum(spcs)/tot_noise	465	# print ('stop 2.3')
	405	snrdB=10.*numpy.log10(snr)	466	# snr = sum(spcs) / tot_noise
	406		467	# snrdB = 10.*numpy.log10(snr)
	407	if snrdB < SNRlimit :	468	#print ('stop 3')
	408	snr = numpy.NaN	469	# if snrdB < SNRlimit :
	409	~~SPC_ch1[:,~~ht] = 0#numpy.NaN	470	# snr = numpy.NaN
	410	SPC_ch1[:,ht] = 0#numpy.NaN	471	# SPC_ch1[:,ht] = 0#numpy.NaN
	411	SPCparam = (SPC_ch1,SPC_ch2)	472	# SPC_ch1[:,ht] = 0#numpy.NaN
	412	continue	473	# SPCparam = (SPC_ch1,SPC_ch2)
			474	# print ('SNR less than SNRth')
			475	# continue
	413		476
	414		477
	415	#if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:	478	#if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
	416	# return [None,]4,[None,]4,None,snrdB,None,None,[None,]5,[None,]9,None	479	# return [None,]4,[None,]4,None,snrdB,None,None,[None,]5,[None,]9,None
	417		480	# print ('stop 4')
	418	cummax=max(cum);	481	cummax = max(cum)
	419	epsi=0.08*fatspectra # cumsum to narrow down the energy region	482	epsi = 0.08 * fatspectra # cumsum to narrow down the energy region
	420	cumlo=cummax*epsi;	483	cumlo = cummax * epsi
	421	cumhi=cummax*(1-epsi)	484	cumhi = cummax * (1-epsi)
	422	powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])	485	powerindex = numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
	423		486
	424		487	# print ('stop 5')
	425	if len(powerindex) < 1:# case for powerindex 0	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	continue	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	498	# print ('stop 6')
	432	~~secondpeak~~=~~powerhi~~-powerwidth/10.#~~second~~ gaussian energy location	499	firstpeak = powerlo + powerwidth/10.# first gaussian energy location
	433	midpeak=(firstpeak+secondpeak)/2.	500	secondpeak = powerhi - powerwidth/10. #second gaussian energy location
	434	firstamp=spcs[int(firstpeak)]	501	midpeak = (firstpeak + secondpeak)/2.
	435	~~secondamp~~=spcs[int(~~second~~peak)]	502	firstamp = spcs[int(firstpeak)]
	436	~~midamp~~=spcs[int(midpeak)]	503	secondamp = spcs[int(secondpeak)]
			504	midamp = spcs[int(midpeak)]
	437		505
	438	x=numpy.arange( self.Num_Bin )	506	y_data = spc + wnoise
	439	y_data=spc+wnoise
	440		507
	441	''' single Gaussian '''	508	''' single Gaussian '''
	442	shift0=numpy.mod(midpeak+minx, self.Num_Bin )	509	shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
	443	width0=powerwidth/4.#Initialization entire power of spectrum divided by 4	510	width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
	444	power0=2.	511	power0 = 2.
	445	amplitude0=midamp	512	amplitude0 = midamp
	446	state0=[shift0,width0,amplitude0,power0,wnoise]	513	state0 = [shift0,width0,amplitude0,power0,wnoise]
	447	bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))	514	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)	515	lsq1 = fmin_l_bfgs_b(self.misfit1, state0, args=(y_data,x,num_intg), bounds=bnds, approx_grad=True)
	449		516	# print ('stop 7.1')
	450	chiSq1=lsq1[1];	517	# print (bnds)
	451		518
	452		519	chiSq1=lsq1[1]
			520
			521	# print ('stop 8')
	453	if fatspectra<1.0 and powerwidth<4:	522	if fatspectra<1.0 and powerwidth<4:
	454	choice=0	523	choice=0
	455	Amplitude0=lsq1[0][2]	524	Amplitude0=lsq1[0][2]
	@@ -464,133 +533,378 class GaussianFit(Operation):
	464	#return (numpy.array([shift0,width0,Amplitude0,p0]),	533	#return (numpy.array([shift0,width0,Amplitude0,p0]),
	465	# numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)	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	#shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)	538	#shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
	469	shift0=numpy.mod(firstpeak+minx, self.Num_Bin );	539	shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
	470	shift1=numpy.mod(secondpeak+minx, self.Num_Bin )	540	shift1 = numpy.mod(secondpeak+minx, self.Num_Bin )
	471	width0=powerwidth/6.;	541	width0 = powerwidth/6.
	472	width1=width0	542	width1 = width0
	473	power0=2.;	543	power0 = 2.
	474	power1=power0	544	power1 = power0
	475	amplitude0=firstamp;	545	amplitude0 = firstamp
	476	amplitude1=secondamp	546	amplitude1 = secondamp
	477	state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]	547	state0 = [shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
	478	#bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))	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	#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))	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	lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )	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];	558	# print ('stop 12')
	486
	487
	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	if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error	563	if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
	492	if oneG:	564	if oneG:
	493	choice=0	565	choice = 0
	494	else:	566	else:
	495	w1=lsq2[0][1]; w2=lsq2[0][5]	567	w1 = lsq2[0][1]; w2 = lsq2[0][5]
	496	a1=lsq2[0][2]; a2=lsq2[0][6]	568	a1 = lsq2[0][2]; a2 = lsq2[0][6]
	497	p1=lsq2[0][3]; p2=lsq2[0][7]	569	p1 = lsq2[0][3]; p2 = lsq2[0][7]
	498	s1=(2*(1+1./p1))scipy.special.gamma(1./p1)/p1;	570	s1 = (2*(1+1./p1))scipy.special.gamma(1./p1)/p1
	499	s2=(2*(1+1./p2))scipy.special.gamma(1./p2)/p2;	571	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	572	gp1 = a1w1s1; gp2 = a2w2s2 # power content of each ggaussian with proper p scaling
	501		573
	502	if gp1>gp2:	574	if gp1>gp2:
	503	if a1>0.7*a2:	575	if a1>0.7*a2:
	504	choice=1	576	choice = 1
	505	else:	577	else:
	506	choice=2	578	choice = 2
	507	elif gp2>gp1:	579	elif gp2>gp1:
	508	if a2>0.7*a1:	580	if a2>0.7*a1:
	509	choice=2	581	choice = 2
	510	else:	582	else:
	511	choice=1	583	choice = 1
	512	else:	584	else:
	513	choice=numpy.argmax([a1,a2])+1	585	choice = numpy.argmax([a1,a2])+1
	514	#else:	586	#else:
	515	#choice=argmin([std2a,std2b])+1	587	#choice=argmin([std2a,std2b])+1
	516		588
	517	else: # with low SNR go to the most energetic peak	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]])	590	choice = numpy.argmax([lsq1[0][2]lsq1[0][1],lsq2[0][2]lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
	519		591
	520		592	# print ('stop 14')
	521	shift0=lsq2[0][0];	593	shift0 = lsq2[0][0]
	522	vel0=Vrange[0] + shift0*(~~Vrange~~[1]-~~Vrange~~[0])	594	vel0 = Vrange[0] + shift0 * deltav
	523	shift1=lsq2[0][4];	595	shift1 = lsq2[0][4]
	524	vel1=Vrange[0] + shift1*(~~Vrange~~[1]-~~Vrange~~[0])	596	# vel1=Vrange[0] + shift1 * deltav
	525		597
	526	max_vel = 1.0	598	# max_vel = 1.0
	527		599	# Va = max(Vrange)
			600	# deltav = Vrange[1]-Vrange[0]
			601	# print ('stop 15')
	528	#first peak will be 0, second peak will be 1	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	603	# if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.Scipión 19.03.2021
	530	shift0=lsq2[0][0]	604	if vel0 > -Va and vel0 < Va : #first peak is in the correct range
	531	~~width0~~=lsq2[0][1]	605	shift0 = lsq2[0][0]
	532	~~Amplitude0~~=lsq2[0][2]	606	width0 = lsq2[0][1]
	533	p0=lsq2[0][3]	607	Amplitude0 = lsq2[0][2]
	534		608	p0 = lsq2[0][3]
	535	shift1=lsq2[0][4]	609
	536	~~width1~~=lsq2[0][5]	610	shift1 = lsq2[0][4]
	537	~~Amplitude1~~=lsq2[0][6]	611	width1 = lsq2[0][5]
	538	p1=lsq2[0][7]	612	Amplitude1 = lsq2[0][6]
	539	~~noise~~=lsq2[0][8]	613	p1 = lsq2[0][7]
			614	noise = lsq2[0][8]
	540	else:	615	else:
	541	shift1=lsq2[0][0]	616	shift1 = lsq2[0][0]
	542	width1=lsq2[0][1]	617	width1 = lsq2[0][1]
	543	Amplitude1=lsq2[0][2]	618	Amplitude1 = lsq2[0][2]
	544	p1=lsq2[0][3]	619	p1 = lsq2[0][3]
	545		620
	546	shift0=lsq2[0][4]	621	shift0 = lsq2[0][4]
	547	width0=lsq2[0][5]	622	width0 = lsq2[0][5]
	548	Amplitude0=lsq2[0][6]	623	Amplitude0 = lsq2[0][6]
	549	p0=lsq2[0][7]	624	p0 = lsq2[0][7]
	550	noise=lsq2[0][8]	625	noise = lsq2[0][8]
	551		626
	552	if Amplitude0<0.05: # in case the peak is noise	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	if Amplitude1<0.05:	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	661	def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
	559	SPC_ch2[:,ht] = noise + Amplitude1numpy.exp(-0.5(abs(x-shift1))/width1)**p1	662	shift0, width0, amplitude0, power0, shift1, width1, amplitude1, power1, noise = state
	560	SPCparam = (SPC_ch1,SPC_ch2)	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):	674	return num_intgsum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))*2)#/(64-5.) # /(64-5.) can be commented
	566	shift0,width0,amplitude0,power0,noise=state
	567	model0=amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)
	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	680	class Oblique_Gauss_Fit(Operation):
	575	shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state	681
	576	model0=amplitude0numpy.exp(-0.5abs((x-shift0)/width0)**power0)	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	# error[indInvalid1] = 13	4312	# error[indInvalid1] = 13
	3999	#	4313	#
	4000	# return heights, error	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

@@ -736,4 +736,4 class SpectraLagsProc(ProcessingUnit):
736		736
737	self.dataOut.noise_estimation = noise.copy()	737	self.dataOut.noise_estimation = noise.copy()
738		738
739	return 1 No newline at end of file	739	return 1

schainpy/model/proc/jroproc_voltage.py 0 0

	1	NO CONTENT: modified file			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			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			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