##// END OF EJS Templates
New GUI added (kivy framework), clean code, delete unused files
New GUI added (kivy framework), clean code, delete unused files

File last commit:

r1285:7f395116f8a5
r1288:8eedfb38dc53
Show More
jroplot_spectra.py
650 lines | 21.2 KiB | text/x-python | PythonLexer
'''
Created on Jul 9, 2014
Modified on May 10, 2020
@author: Juan C. Espinoza
'''
import os
import datetime
import numpy
from schainpy.model.graphics.jroplot_base import Plot, plt
class SpectraPlot(Plot):
'''
Plot for Spectra data
'''
CODE = 'spc'
colormap = 'jet'
plot_name = 'Spectra'
plot_type = 'pcolor'
def setup(self):
self.nplots = len(self.data.channels)
self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
self.height = 3 * self.nrows
self.cb_label = 'dB'
if self.showprofile:
self.width = 4 * self.ncols
else:
self.width = 3.5 * self.ncols
self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
self.ylabel = 'Range [km]'
def plot(self):
if self.xaxis == "frequency":
x = self.data.xrange[0]
self.xlabel = "Frequency (kHz)"
elif self.xaxis == "time":
x = self.data.xrange[1]
self.xlabel = "Time (ms)"
else:
x = self.data.xrange[2]
self.xlabel = "Velocity (m/s)"
if self.CODE == 'spc_moments':
x = self.data.xrange[2]
self.xlabel = "Velocity (m/s)"
self.titles = []
y = self.data.heights
self.y = y
z = self.data['spc']
for n, ax in enumerate(self.axes):
noise = self.data['noise'][n][-1]
if self.CODE == 'spc_moments':
mean = self.data['moments'][n, :, 1, :][-1]
if ax.firsttime:
self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
self.xmin = self.xmin if self.xmin else -self.xmax
self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
ax.plt = ax.pcolormesh(x, y, z[n].T,
vmin=self.zmin,
vmax=self.zmax,
cmap=plt.get_cmap(self.colormap)
)
if self.showprofile:
ax.plt_profile = self.pf_axes[n].plot(
self.data['rti'][n][-1], y)[0]
ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
color="k", linestyle="dashed", lw=1)[0]
if self.CODE == 'spc_moments':
ax.plt_mean = ax.plot(mean, y, color='k')[0]
else:
ax.plt.set_array(z[n].T.ravel())
if self.showprofile:
ax.plt_profile.set_data(self.data['rti'][n][-1], y)
ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
if self.CODE == 'spc_moments':
ax.plt_mean.set_data(mean, y)
self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
class CrossSpectraPlot(Plot):
CODE = 'cspc'
colormap = 'jet'
plot_name = 'CrossSpectra'
plot_type = 'pcolor'
zmin_coh = None
zmax_coh = None
zmin_phase = None
zmax_phase = None
def setup(self):
self.ncols = 4
self.nrows = len(self.data.pairs)
self.nplots = self.nrows * 4
self.width = 3.4 * self.ncols
self.height = 3 * self.nrows
self.ylabel = 'Range [km]'
self.showprofile = False
self.plots_adjust.update({'bottom': 0.08})
def plot(self):
if self.xaxis == "frequency":
x = self.data.xrange[0]
self.xlabel = "Frequency (kHz)"
elif self.xaxis == "time":
x = self.data.xrange[1]
self.xlabel = "Time (ms)"
else:
x = self.data.xrange[2]
self.xlabel = "Velocity (m/s)"
self.titles = []
y = self.data.heights
self.y = y
spc = self.data['spc']
cspc = self.data['cspc']
for n in range(self.nrows):
noise = self.data['noise'][n][-1]
pair = self.data.pairs[n]
ax = self.axes[4 * n]
spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
if ax.firsttime:
self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
self.xmin = self.xmin if self.xmin else -self.xmax
self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
ax.plt = ax.pcolormesh(x , y , spc0.T,
vmin=self.zmin,
vmax=self.zmax,
cmap=plt.get_cmap(self.colormap)
)
else:
ax.plt.set_array(spc0.T.ravel())
self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
ax = self.axes[4 * n + 1]
spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
if ax.firsttime:
ax.plt = ax.pcolormesh(x , y, spc1.T,
vmin=self.zmin,
vmax=self.zmax,
cmap=plt.get_cmap(self.colormap)
)
else:
ax.plt.set_array(spc1.T.ravel())
self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
coh = numpy.abs(out)
phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
ax = self.axes[4 * n + 2]
if ax.firsttime:
ax.plt = ax.pcolormesh(x, y, coh.T,
vmin=0,
vmax=1,
cmap=plt.get_cmap(self.colormap_coh)
)
else:
ax.plt.set_array(coh.T.ravel())
self.titles.append(
'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
ax = self.axes[4 * n + 3]
if ax.firsttime:
ax.plt = ax.pcolormesh(x, y, phase.T,
vmin=-180,
vmax=180,
cmap=plt.get_cmap(self.colormap_phase)
)
else:
ax.plt.set_array(phase.T.ravel())
self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
class RTIPlot(Plot):
'''
Plot for RTI data
'''
CODE = 'rti'
colormap = 'jet'
plot_name = 'RTI'
plot_type = 'pcolorbuffer'
def setup(self):
self.xaxis = 'time'
self.ncols = 1
self.nrows = len(self.data.channels)
self.nplots = len(self.data.channels)
self.ylabel = 'Range [km]'
self.xlabel = 'Time'
self.cb_label = 'dB'
self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
self.titles = ['{} Channel {}'.format(
self.CODE.upper(), x) for x in range(self.nrows)]
def plot(self):
self.x = self.data.times
self.y = self.data.heights
self.z = self.data[self.CODE]
self.z = numpy.ma.masked_invalid(self.z)
if self.decimation is None:
x, y, z = self.fill_gaps(self.x, self.y, self.z)
else:
x, y, z = self.fill_gaps(*self.decimate())
for n, ax in enumerate(self.axes):
self.zmin = self.zmin if self.zmin else numpy.min(self.z)
self.zmax = self.zmax if self.zmax else numpy.max(self.z)
if ax.firsttime:
ax.plt = ax.pcolormesh(x, y, z[n].T,
vmin=self.zmin,
vmax=self.zmax,
cmap=plt.get_cmap(self.colormap)
)
if self.showprofile:
ax.plot_profile = self.pf_axes[n].plot(
self.data['rti'][n][-1], self.y)[0]
ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
color="k", linestyle="dashed", lw=1)[0]
else:
ax.collections.remove(ax.collections[0])
ax.plt = ax.pcolormesh(x, y, z[n].T,
vmin=self.zmin,
vmax=self.zmax,
cmap=plt.get_cmap(self.colormap)
)
if self.showprofile:
ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
ax.plot_noise.set_data(numpy.repeat(
self.data['noise'][n][-1], len(self.y)), self.y)
class CoherencePlot(RTIPlot):
'''
Plot for Coherence data
'''
CODE = 'coh'
plot_name = 'Coherence'
def setup(self):
self.xaxis = 'time'
self.ncols = 1
self.nrows = len(self.data.pairs)
self.nplots = len(self.data.pairs)
self.ylabel = 'Range [km]'
self.xlabel = 'Time'
self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
if self.CODE == 'coh':
self.cb_label = ''
self.titles = [
'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
else:
self.cb_label = 'Degrees'
self.titles = [
'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
class PhasePlot(CoherencePlot):
'''
Plot for Phase map data
'''
CODE = 'phase'
colormap = 'seismic'
plot_name = 'Phase'
class NoisePlot(Plot):
'''
Plot for noise
'''
CODE = 'noise'
plot_name = 'Noise'
plot_type = 'scatterbuffer'
def setup(self):
self.xaxis = 'time'
self.ncols = 1
self.nrows = 1
self.nplots = 1
self.ylabel = 'Intensity [dB]'
self.titles = ['Noise']
self.colorbar = False
def plot(self):
x = self.data.times
xmin = self.data.min_time
xmax = xmin + self.xrange * 60 * 60
Y = self.data[self.CODE]
if self.axes[0].firsttime:
for ch in self.data.channels:
y = Y[ch]
self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
plt.legend()
else:
for ch in self.data.channels:
y = Y[ch]
self.axes[0].lines[ch].set_data(x, y)
self.ymin = numpy.nanmin(Y) - 5
self.ymax = numpy.nanmax(Y) + 5
class PowerProfilePlot(Plot):
CODE = 'spcprofile'
plot_name = 'Power Profile'
plot_type = 'scatter'
buffering = False
def setup(self):
self.ncols = 1
self.nrows = 1
self.nplots = 1
self.height = 4
self.width = 3
self.ylabel = 'Range [km]'
self.xlabel = 'Intensity [dB]'
self.titles = ['Power Profile']
self.colorbar = False
def plot(self):
y = self.data.heights
self.y = y
x = self.data['spcprofile']
if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
if self.axes[0].firsttime:
for ch in self.data.channels:
self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
plt.legend()
else:
for ch in self.data.channels:
self.axes[0].lines[ch].set_data(x[ch], y)
class SpectraCutPlot(Plot):
CODE = 'spc_cut'
plot_name = 'Spectra Cut'
plot_type = 'scatter'
buffering = False
def setup(self):
self.nplots = len(self.data.channels)
self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
self.width = 3.4 * self.ncols + 1.5
self.height = 3 * self.nrows
self.ylabel = 'Power [dB]'
self.colorbar = False
self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
def plot(self):
if self.xaxis == "frequency":
x = self.data.xrange[0][1:]
self.xlabel = "Frequency (kHz)"
elif self.xaxis == "time":
x = self.data.xrange[1]
self.xlabel = "Time (ms)"
else:
x = self.data.xrange[2]
self.xlabel = "Velocity (m/s)"
self.titles = []
y = self.data.heights
#self.y = y
z = self.data['spc_cut']
if self.height_index:
index = numpy.array(self.height_index)
else:
index = numpy.arange(0, len(y), int((len(y))/9))
for n, ax in enumerate(self.axes):
if ax.firsttime:
self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
self.xmin = self.xmin if self.xmin else -self.xmax
self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
ax.plt = ax.plot(x, z[n, :, index].T)
labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
self.figures[0].legend(ax.plt, labels, loc='center right')
else:
for i, line in enumerate(ax.plt):
line.set_data(x, z[n, :, i])
self.titles.append('CH {}'.format(n))
class BeaconPhase(Plot):
__isConfig = None
__nsubplots = None
PREFIX = 'beacon_phase'
def __init__(self):
Plot.__init__(self)
self.timerange = 24*60*60
self.isConfig = False
self.__nsubplots = 1
self.counter_imagwr = 0
self.WIDTH = 800
self.HEIGHT = 400
self.WIDTHPROF = 120
self.HEIGHTPROF = 0
self.xdata = None
self.ydata = None
self.PLOT_CODE = BEACON_CODE
self.FTP_WEI = None
self.EXP_CODE = None
self.SUB_EXP_CODE = None
self.PLOT_POS = None
self.filename_phase = None
self.figfile = None
self.xmin = None
self.xmax = None
def getSubplots(self):
ncol = 1
nrow = 1
return nrow, ncol
def setup(self, id, nplots, wintitle, showprofile=True, show=True):
self.__showprofile = showprofile
self.nplots = nplots
ncolspan = 7
colspan = 6
self.__nsubplots = 2
self.createFigure(id = id,
wintitle = wintitle,
widthplot = self.WIDTH+self.WIDTHPROF,
heightplot = self.HEIGHT+self.HEIGHTPROF,
show=show)
nrow, ncol = self.getSubplots()
self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
def save_phase(self, filename_phase):
f = open(filename_phase,'w+')
f.write('\n\n')
f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
f.close()
def save_data(self, filename_phase, data, data_datetime):
f=open(filename_phase,'a')
timetuple_data = data_datetime.timetuple()
day = str(timetuple_data.tm_mday)
month = str(timetuple_data.tm_mon)
year = str(timetuple_data.tm_year)
hour = str(timetuple_data.tm_hour)
minute = str(timetuple_data.tm_min)
second = str(timetuple_data.tm_sec)
f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
f.close()
def plot(self):
log.warning('TODO: Not yet implemented...')
def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
timerange=None,
save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
server=None, folder=None, username=None, password=None,
ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
if dataOut.flagNoData:
return dataOut
if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
return
if pairsList == None:
pairsIndexList = dataOut.pairsIndexList[:10]
else:
pairsIndexList = []
for pair in pairsList:
if pair not in dataOut.pairsList:
raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
pairsIndexList.append(dataOut.pairsList.index(pair))
if pairsIndexList == []:
return
# if len(pairsIndexList) > 4:
# pairsIndexList = pairsIndexList[0:4]
hmin_index = None
hmax_index = None
if hmin != None and hmax != None:
indexes = numpy.arange(dataOut.nHeights)
hmin_list = indexes[dataOut.heightList >= hmin]
hmax_list = indexes[dataOut.heightList <= hmax]
if hmin_list.any():
hmin_index = hmin_list[0]
if hmax_list.any():
hmax_index = hmax_list[-1]+1
x = dataOut.getTimeRange()
#y = dataOut.getHeiRange()
thisDatetime = dataOut.datatime
title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
xlabel = "Local Time"
ylabel = "Phase (degrees)"
update_figfile = False
nplots = len(pairsIndexList)
#phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
phase_beacon = numpy.zeros(len(pairsIndexList))
for i in range(nplots):
pair = dataOut.pairsList[pairsIndexList[i]]
ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
if dataOut.beacon_heiIndexList:
phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
else:
phase_beacon[i] = numpy.average(phase)
if not self.isConfig:
nplots = len(pairsIndexList)
self.setup(id=id,
nplots=nplots,
wintitle=wintitle,
showprofile=showprofile,
show=show)
if timerange != None:
self.timerange = timerange
self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
if ymin == None: ymin = 0
if ymax == None: ymax = 360
self.FTP_WEI = ftp_wei
self.EXP_CODE = exp_code
self.SUB_EXP_CODE = sub_exp_code
self.PLOT_POS = plot_pos
self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
self.isConfig = True
self.figfile = figfile
self.xdata = numpy.array([])
self.ydata = numpy.array([])
update_figfile = True
#open file beacon phase
path = '%s%03d' %(self.PREFIX, self.id)
beacon_file = os.path.join(path,'%s.txt'%self.name)
self.filename_phase = os.path.join(figpath,beacon_file)
#self.save_phase(self.filename_phase)
#store data beacon phase
#self.save_data(self.filename_phase, phase_beacon, thisDatetime)
self.setWinTitle(title)
title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
axes = self.axesList[0]
self.xdata = numpy.hstack((self.xdata, x[0:1]))
if len(self.ydata)==0:
self.ydata = phase_beacon.reshape(-1,1)
else:
self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
XAxisAsTime=True, grid='both'
)
self.draw()
if dataOut.ltctime >= self.xmax:
self.counter_imagwr = wr_period
self.isConfig = False
update_figfile = True
self.save(figpath=figpath,
figfile=figfile,
save=save,
ftp=ftp,
wr_period=wr_period,
thisDatetime=thisDatetime,
update_figfile=update_figfile)
return dataOut