diff --git a/.gitignore b/.gitignore index c78f88d..e0bc00e 100644 --- a/.gitignore +++ b/.gitignore @@ -112,5 +112,3 @@ schainpy/scripts/ .vscode trash *.log -schainpy/scripts/testDigitalRF.py -schainpy/scripts/testDigitalRFWriter.py diff --git a/schainpy/CHANGELOG.md b/schainpy/CHANGELOG.md index 4f4d159..c54d3f0 100644 --- a/schainpy/CHANGELOG.md +++ b/schainpy/CHANGELOG.md @@ -1,5 +1,12 @@ ## CHANGELOG: +### 3.0 +* Python 3.x compatible +* New architecture with multiprocessing and IPC communication +* Add @MPDecorator for multiprocessing Units and Operations +* Added new type of operation `external` for non-locking operations +* New plotting architecture with buffering/throttle capabilities to speed up plots + ### 2.3 * Added support for Madrigal formats (reading/writing). * Added support for reading BLTR parameters (*.sswma). diff --git a/schainpy/admin.py b/schainpy/admin.py index 7dc2984..0bb49fa 100644 --- a/schainpy/admin.py +++ b/schainpy/admin.py @@ -24,7 +24,7 @@ from email.mime.multipart import MIMEMultipart import schainpy from schainpy.utils import log -from schainpy.model.graphics.jroplot_data import popup +from schainpy.model.graphics.jroplot_base import popup def get_path(): ''' diff --git a/schainpy/controller.py b/schainpy/controller.py index 8e39197..42ad19d 100644 --- a/schainpy/controller.py +++ b/schainpy/controller.py @@ -97,9 +97,7 @@ def wait(context): receiver = c.socket(zmq.SUB) receiver.connect('ipc:///tmp/schain_{}_pub'.format(self.id)) receiver.setsockopt(zmq.SUBSCRIBE, self.id.encode()) - log.error('startinggg') msg = receiver.recv_multipart()[1] - #log.error(msg) context.terminate() class ParameterConf(): @@ -1245,7 +1243,7 @@ class Project(Process): try: zmq.proxy(xpub, xsub) - except zmq.ContextTerminated: + except: # zmq.ContextTerminated: xpub.close() xsub.close() @@ -1260,6 +1258,6 @@ class Project(Process): # Iniciar todos los procesos .start(), monitoreo de procesos. ELiminar lo de abajo - log.success('{} finished (time: {}s)'.format( + log.success('{} Done (time: {}s)'.format( self.name, time.time()-self.start_time)) diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index daa22c1..bebece8 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -7,7 +7,9 @@ $Id: JROData.py 173 2012-11-20 15:06:21Z murco $ import copy import numpy import datetime +import json +from schainpy.utils import log from .jroheaderIO import SystemHeader, RadarControllerHeader @@ -79,7 +81,7 @@ def hildebrand_sekhon(data, navg): j = 0 cont = 1 - while((cont==1)and(j (rtest*sump**2)): j = j - 1 - sump = sump - sortdata[j] - sumq = sumq - sortdata[j]**2 + sump = sump - sortdata[j] + sumq = sumq - sortdata[j]**2 cont = 0 j += 1 - lnoise = sump /j + lnoise = sump / j return lnoise @@ -147,85 +149,48 @@ class JROData(GenericData): # m_ProcessingHeader = ProcessingHeader() systemHeaderObj = SystemHeader() - radarControllerHeaderObj = RadarControllerHeader() - # data = None - type = None - datatype = None # dtype but in string - # dtype = None - # nChannels = None - # nHeights = None - nProfiles = None - heightList = None - channelList = None - flagDiscontinuousBlock = False - useLocalTime = False - utctime = None - timeZone = None - dstFlag = None - errorCount = None - blocksize = None - # nCode = None -# # nBaud = None -# # code = None - flagDecodeData = False # asumo q la data no esta decodificada - flagDeflipData = False # asumo q la data no esta sin flip - flagShiftFFT = False - # ippSeconds = None - # timeInterval = None - nCohInt = None - # noise = None - windowOfFilter = 1 - # Speed of ligth C = 3e8 - frequency = 49.92e6 - realtime = False - beacon_heiIndexList = None - last_block = None - blocknow = None - azimuth = None - zenith = None - beam = Beam() - profileIndex = None - - error = (0, '') + error = None + data = None + nmodes = None def __str__(self): @@ -395,53 +360,29 @@ class Voltage(JROData): ''' self.useLocalTime = True - self.radarControllerHeaderObj = RadarControllerHeader() - self.systemHeaderObj = SystemHeader() - self.type = "Voltage" - self.data = None - # self.dtype = None - # self.nChannels = 0 - # self.nHeights = 0 - self.nProfiles = None - - self.heightList = None - + self.heightList = Non self.channelList = None - # self.channelIndexList = None - self.flagNoData = True - self.flagDiscontinuousBlock = False - self.utctime = None - self.timeZone = None - self.dstFlag = None - self.errorCount = None - self.nCohInt = None - self.blocksize = None - self.flagDecodeData = False # asumo q la data no esta decodificada - self.flagDeflipData = False # asumo q la data no esta sin flip - self.flagShiftFFT = False - self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil - self.profileIndex = 0 def getNoisebyHildebrand(self, channel=None): @@ -505,32 +446,20 @@ class Spectra(JROData): # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas) data_spc = None - # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas) data_cspc = None - # data dc es un numpy array de 2 dmensiones (canales, alturas) data_dc = None - # data power data_pwr = None - nFFTPoints = None - # nPairs = None - pairsList = None - nIncohInt = None - wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia - nCohInt = None # se requiere para determinar el valor de timeInterval - ippFactor = None - profileIndex = 0 - plotting = "spectra" def __init__(self): @@ -539,59 +468,32 @@ class Spectra(JROData): ''' self.useLocalTime = True - self.radarControllerHeaderObj = RadarControllerHeader() - self.systemHeaderObj = SystemHeader() - self.type = "Spectra" - # self.data = None - # self.dtype = None - # self.nChannels = 0 - # self.nHeights = 0 - self.nProfiles = None - self.heightList = None - self.channelList = None - # self.channelIndexList = None - self.pairsList = None - self.flagNoData = True - self.flagDiscontinuousBlock = False - self.utctime = None - self.nCohInt = None - self.nIncohInt = None - self.blocksize = None - self.nFFTPoints = None - self.wavelength = None - self.flagDecodeData = False # asumo q la data no esta decodificada - self.flagDeflipData = False # asumo q la data no esta sin flip - self.flagShiftFFT = False - self.ippFactor = 1 - #self.noise = None - self.beacon_heiIndexList = [] - self.noise_estimation = None def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None): @@ -652,9 +554,12 @@ class Spectra(JROData): deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor) velrange = deltav * (numpy.arange(self.nFFTPoints + - extrapoints) - self.nFFTPoints / 2.) # - deltav/2 + extrapoints) - self.nFFTPoints / 2.) - return velrange + if self.nmodes: + return velrange/self.nmodes + else: + return velrange def getNPairs(self): @@ -692,7 +597,8 @@ class Spectra(JROData): def getTimeInterval(self): - timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor + timeInterval = self.ippSeconds * self.nCohInt * \ + self.nIncohInt * self.nProfiles * self.ippFactor return timeInterval @@ -755,19 +661,12 @@ class Spectra(JROData): class SpectraHeis(Spectra): data_spc = None - data_cspc = None - data_dc = None - nFFTPoints = None - # nPairs = None - pairsList = None - nCohInt = None - nIncohInt = None def __init__(self): @@ -830,36 +729,21 @@ class SpectraHeis(Spectra): class Fits(JROData): heightList = None - channelList = None - flagNoData = True - flagDiscontinuousBlock = False - useLocalTime = False - utctime = None - timeZone = None - # ippSeconds = None - # timeInterval = None - nCohInt = None - nIncohInt = None - noise = None - windowOfFilter = 1 - # Speed of ligth C = 3e8 - frequency = 49.92e6 - realtime = False def __init__(self): @@ -978,33 +862,19 @@ class Fits(JROData): class Correlation(JROData): noise = None - SNR = None - #-------------------------------------------------- - mode = None - split = False - data_cf = None - lags = None - lagRange = None - pairsList = None - normFactor = None - #-------------------------------------------------- - # calculateVelocity = None - nLags = None - nPairs = None - nAvg = None def __init__(self): @@ -1068,7 +938,8 @@ class Correlation(JROData): ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc if ind_vel[0] < 0: - ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof + ind_vel[list(range(0, 1))] = ind_vel[list( + range(0, 1))] + self.num_prof if mode == 1: jspectra[:, freq_dc, :] = ( @@ -1154,55 +1025,30 @@ class Correlation(JROData): class Parameters(Spectra): experimentInfo = None # Information about the experiment - # Information from previous data - inputUnit = None # Type of data to be processed - operation = None # Type of operation to parametrize - # normFactor = None #Normalization Factor - groupList = None # List of Pairs, Groups, etc - # Parameters - data_param = None # Parameters obtained - data_pre = None # Data Pre Parametrization - data_SNR = None # Signal to Noise Ratio - # heightRange = None #Heights - abscissaList = None # Abscissa, can be velocities, lags or time - # noise = None #Noise Potency - utctimeInit = None # Initial UTC time - paramInterval = None # Time interval to calculate Parameters in seconds - useLocalTime = True - # Fitting - data_error = None # Error of the estimation - constants = None - library = None - # Output signal - outputInterval = None # Time interval to calculate output signal in seconds - data_output = None # Out signal - nAvg = None - noise_estimation = None - GauSPC = None # Fit gaussian SPC def __init__(self): @@ -1248,4 +1094,252 @@ class Parameters(Spectra): return self.spc_noise timeInterval = property(getTimeInterval) - noise = property(getNoise, setValue, "I'm the 'Noise' property.") \ No newline at end of file + noise = property(getNoise, setValue, "I'm the 'Noise' property.") + + +class PlotterData(object): + ''' + Object to hold data to be plotted + ''' + + MAXNUMX = 100 + MAXNUMY = 100 + + def __init__(self, code, throttle_value, exp_code, buffering=True): + + self.throttle = throttle_value + self.exp_code = exp_code + self.buffering = buffering + self.ready = False + self.localtime = False + self.data = {} + self.meta = {} + self.__times = [] + self.__heights = [] + + if 'snr' in code: + self.plottypes = ['snr'] + elif code == 'spc': + self.plottypes = ['spc', 'noise', 'rti'] + elif code == 'rti': + self.plottypes = ['noise', 'rti'] + else: + self.plottypes = [code] + + for plot in self.plottypes: + self.data[plot] = {} + + def __str__(self): + dum = ['{}{}'.format(key, self.shape(key)) for key in self.data] + return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times)) + + def __len__(self): + return len(self.__times) + + def __getitem__(self, key): + + if key not in self.data: + raise KeyError(log.error('Missing key: {}'.format(key))) + if 'spc' in key or not self.buffering: + ret = self.data[key] + else: + ret = numpy.array([self.data[key][x] for x in self.times]) + if ret.ndim > 1: + ret = numpy.swapaxes(ret, 0, 1) + return ret + + def __contains__(self, key): + return key in self.data + + def setup(self): + ''' + Configure object + ''' + + self.type = '' + self.ready = False + self.data = {} + self.__times = [] + self.__heights = [] + self.__all_heights = set() + for plot in self.plottypes: + if 'snr' in plot: + plot = 'snr' + self.data[plot] = {} + + if 'spc' in self.data or 'rti' in self.data: + self.data['noise'] = {} + if 'noise' not in self.plottypes: + self.plottypes.append('noise') + + def shape(self, key): + ''' + Get the shape of the one-element data for the given key + ''' + + if len(self.data[key]): + if 'spc' in key or not self.buffering: + return self.data[key].shape + return self.data[key][self.__times[0]].shape + return (0,) + + def update(self, dataOut, tm): + ''' + Update data object with new dataOut + ''' + + if tm in self.__times: + return + + self.type = dataOut.type + self.parameters = getattr(dataOut, 'parameters', []) + if hasattr(dataOut, 'pairsList'): + self.pairs = dataOut.pairsList + if hasattr(dataOut, 'meta'): + self.meta = dataOut.meta + self.channels = dataOut.channelList + self.interval = dataOut.getTimeInterval() + self.localtime = dataOut.useLocalTime + if 'spc' in self.plottypes or 'cspc' in self.plottypes: + self.xrange = (dataOut.getFreqRange(1)/1000., + dataOut.getAcfRange(1), dataOut.getVelRange(1)) + self.__heights.append(dataOut.heightList) + self.__all_heights.update(dataOut.heightList) + self.__times.append(tm) + + for plot in self.plottypes: + if plot == 'spc': + z = dataOut.data_spc/dataOut.normFactor + buffer = 10*numpy.log10(z) + if plot == 'cspc': + buffer = dataOut.data_cspc + if plot == 'noise': + buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + if plot == 'rti': + buffer = dataOut.getPower() + if plot == 'snr_db': + buffer = dataOut.data_SNR + if plot == 'snr': + buffer = 10*numpy.log10(dataOut.data_SNR) + if plot == 'dop': + buffer = 10*numpy.log10(dataOut.data_DOP) + if plot == 'mean': + buffer = dataOut.data_MEAN + if plot == 'std': + buffer = dataOut.data_STD + if plot == 'coh': + buffer = dataOut.getCoherence() + if plot == 'phase': + buffer = dataOut.getCoherence(phase=True) + if plot == 'output': + buffer = dataOut.data_output + if plot == 'param': + buffer = dataOut.data_param + + if 'spc' in plot: + self.data[plot] = buffer + else: + if self.buffering: + self.data[plot][tm] = buffer + else: + self.data[plot] = buffer + + def normalize_heights(self): + ''' + Ensure same-dimension of the data for different heighList + ''' + + H = numpy.array(list(self.__all_heights)) + H.sort() + for key in self.data: + shape = self.shape(key)[:-1] + H.shape + for tm, obj in list(self.data[key].items()): + h = self.__heights[self.__times.index(tm)] + if H.size == h.size: + continue + index = numpy.where(numpy.in1d(H, h))[0] + dummy = numpy.zeros(shape) + numpy.nan + if len(shape) == 2: + dummy[:, index] = obj + else: + dummy[index] = obj + self.data[key][tm] = dummy + + self.__heights = [H for tm in self.__times] + + def jsonify(self, decimate=False): + ''' + Convert data to json + ''' + + data = {} + tm = self.times[-1] + dy = int(self.heights.size/self.MAXNUMY) + 1 + for key in self.data: + if key in ('spc', 'cspc') or not self.buffering: + dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1 + data[key] = self.roundFloats( + self.data[key][::, ::dx, ::dy].tolist()) + else: + data[key] = self.roundFloats(self.data[key][tm].tolist()) + + ret = {'data': data} + ret['exp_code'] = self.exp_code + ret['time'] = float(tm) + ret['interval'] = float(self.interval) + ret['localtime'] = self.localtime + ret['yrange'] = self.roundFloats(self.heights[::dy].tolist()) + if 'spc' in self.data or 'cspc' in self.data: + ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist()) + else: + ret['xrange'] = [] + if hasattr(self, 'pairs'): + ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs] + else: + ret['pairs'] = [] + + for key, value in list(self.meta.items()): + ret[key] = value + + return json.dumps(ret) + + @property + def times(self): + ''' + Return the list of times of the current data + ''' + + ret = numpy.array(self.__times) + ret.sort() + return ret + + @property + def min_time(self): + ''' + Return the minimun time value + ''' + + return self.times[0] + + @property + def max_time(self): + ''' + Return the maximun time value + ''' + + return self.times[-1] + + @property + def heights(self): + ''' + Return the list of heights of the current data + ''' + + return numpy.array(self.__heights[-1]) + + @staticmethod + def roundFloats(obj): + if isinstance(obj, list): + return list(map(PlotterData.roundFloats, obj)) + elif isinstance(obj, float): + return round(obj, 2) diff --git a/schainpy/model/graphics/__init__.py b/schainpy/model/graphics/__init__.py index 98b5033..439b190 100644 --- a/schainpy/model/graphics/__init__.py +++ b/schainpy/model/graphics/__init__.py @@ -4,4 +4,3 @@ from .jroplot_heispectra import * from .jroplot_correlation import * from .jroplot_parameters import * from .jroplot_data import * -from .jroplotter import * diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py new file mode 100644 index 0000000..fd6c1ec --- /dev/null +++ b/schainpy/model/graphics/jroplot_base.py @@ -0,0 +1,803 @@ + +import os +import sys +import zmq +import time +import datetime +from functools import wraps +import numpy +import matplotlib + +if 'BACKEND' in os.environ: + matplotlib.use(os.environ['BACKEND']) +elif 'linux' in sys.platform: + matplotlib.use("TkAgg") +elif 'darwin' in sys.platform: + matplotlib.use('TkAgg') +else: + from schainpy.utils import log + log.warning('Using default Backend="Agg"', 'INFO') + matplotlib.use('Agg') + +import matplotlib.pyplot as plt +from matplotlib.patches import Polygon +from mpl_toolkits.axes_grid1 import make_axes_locatable +from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator + +from schainpy.model.data.jrodata import PlotterData +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator +from schainpy.utils import log + +jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90] +blu_values = matplotlib.pyplot.get_cmap( + 'seismic_r', 20)(numpy.arange(20))[10:15] +ncmap = matplotlib.colors.LinearSegmentedColormap.from_list( + 'jro', numpy.vstack((blu_values, jet_values))) +matplotlib.pyplot.register_cmap(cmap=ncmap) + +CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis', + 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')] + +EARTH_RADIUS = 6.3710e3 + + +def ll2xy(lat1, lon1, lat2, lon2): + + p = 0.017453292519943295 + a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \ + numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2 + r = 12742 * numpy.arcsin(numpy.sqrt(a)) + theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p) + * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p)) + theta = -theta + numpy.pi/2 + return r*numpy.cos(theta), r*numpy.sin(theta) + + +def km2deg(km): + ''' + Convert distance in km to degrees + ''' + + return numpy.rad2deg(km/EARTH_RADIUS) + + +def figpause(interval): + backend = plt.rcParams['backend'] + if backend in matplotlib.rcsetup.interactive_bk: + figManager = matplotlib._pylab_helpers.Gcf.get_active() + if figManager is not None: + canvas = figManager.canvas + if canvas.figure.stale: + canvas.draw() + try: + canvas.start_event_loop(interval) + except: + pass + return + + +def popup(message): + ''' + ''' + + fig = plt.figure(figsize=(12, 8), facecolor='r') + text = '\n'.join([s.strip() for s in message.split(':')]) + fig.text(0.01, 0.5, text, ha='left', va='center', + size='20', weight='heavy', color='w') + fig.show() + figpause(1000) + + +class Throttle(object): + ''' + Decorator that prevents a function from being called more than once every + time period. + To create a function that cannot be called more than once a minute, but + will sleep until it can be called: + @Throttle(minutes=1) + def foo(): + pass + + for i in range(10): + foo() + print "This function has run %s times." % i + ''' + + def __init__(self, seconds=0, minutes=0, hours=0): + self.throttle_period = datetime.timedelta( + seconds=seconds, minutes=minutes, hours=hours + ) + + self.time_of_last_call = datetime.datetime.min + + def __call__(self, fn): + @wraps(fn) + def wrapper(*args, **kwargs): + coerce = kwargs.pop('coerce', None) + if coerce: + self.time_of_last_call = datetime.datetime.now() + return fn(*args, **kwargs) + else: + now = datetime.datetime.now() + time_since_last_call = now - self.time_of_last_call + time_left = self.throttle_period - time_since_last_call + + if time_left > datetime.timedelta(seconds=0): + return + + self.time_of_last_call = datetime.datetime.now() + return fn(*args, **kwargs) + + return wrapper + +def apply_throttle(value): + + @Throttle(seconds=value) + def fnThrottled(fn): + fn() + + return fnThrottled + +@MPDecorator +class Plotter(ProcessingUnit): + ''' + Proccessing unit to handle plot operations + ''' + + def __init__(self): + + ProcessingUnit.__init__(self) + + def setup(self, **kwargs): + + self.connections = 0 + self.web_address = kwargs.get('web_server', False) + self.realtime = kwargs.get('realtime', False) + self.localtime = kwargs.get('localtime', True) + self.buffering = kwargs.get('buffering', True) + self.throttle = kwargs.get('throttle', 2) + self.exp_code = kwargs.get('exp_code', None) + self.set_ready = apply_throttle(self.throttle) + self.dates = [] + self.data = PlotterData( + self.plots, self.throttle, self.exp_code, self.buffering) + self.isConfig = True + + def ready(self): + ''' + Set dataOut ready + ''' + + self.data.ready = True + self.dataOut.data_plt = self.data + + def run(self, realtime=True, localtime=True, buffering=True, + throttle=2, exp_code=None, web_server=None): + + if not self.isConfig: + self.setup(realtime=realtime, localtime=localtime, + buffering=buffering, throttle=throttle, exp_code=exp_code, + web_server=web_server) + + if self.web_address: + log.success( + 'Sending to web: {}'.format(self.web_address), + self.name + ) + self.context = zmq.Context() + self.sender_web = self.context.socket(zmq.REQ) + self.sender_web.connect(self.web_address) + self.poll = zmq.Poller() + self.poll.register(self.sender_web, zmq.POLLIN) + time.sleep(1) + + # t = Thread(target=self.event_monitor, args=(monitor,)) + # t.start() + + self.dataOut = self.dataIn + self.data.ready = False + + if self.dataOut.flagNoData: + coerce = True + else: + coerce = False + + if self.dataOut.type == 'Parameters': + tm = self.dataOut.utctimeInit + else: + tm = self.dataOut.utctime + if self.dataOut.useLocalTime: + if not self.localtime: + tm += time.timezone + dt = datetime.datetime.fromtimestamp(tm).date() + else: + if self.localtime: + tm -= time.timezone + dt = datetime.datetime.utcfromtimestamp(tm).date() + if dt not in self.dates: + if self.data: + self.ready() + self.data.setup() + self.dates.append(dt) + + self.data.update(self.dataOut, tm) + + if False: # TODO check when publishers ends + self.connections -= 1 + if self.connections == 0 and dt in self.dates: + self.data.ended = True + self.ready() + time.sleep(1) + else: + if self.realtime: + self.ready() + if self.web_address: + retries = 5 + while True: + self.sender_web.send(self.data.jsonify()) + socks = dict(self.poll.poll(5000)) + if socks.get(self.sender_web) == zmq.POLLIN: + reply = self.sender_web.recv_string() + if reply == 'ok': + log.log("Response from server ok", self.name) + break + else: + log.warning( + "Malformed reply from server: {}".format(reply), self.name) + + else: + log.warning( + "No response from server, retrying...", self.name) + self.sender_web.setsockopt(zmq.LINGER, 0) + self.sender_web.close() + self.poll.unregister(self.sender_web) + retries -= 1 + if retries == 0: + log.error( + "Server seems to be offline, abandoning", self.name) + self.sender_web = self.context.socket(zmq.REQ) + self.sender_web.connect(self.web_address) + self.poll.register(self.sender_web, zmq.POLLIN) + time.sleep(1) + break + self.sender_web = self.context.socket(zmq.REQ) + self.sender_web.connect(self.web_address) + self.poll.register(self.sender_web, zmq.POLLIN) + time.sleep(1) + else: + self.set_ready(self.ready, coerce=coerce) + + return + + def close(self): + pass + + +@MPDecorator +class Plot(Operation): + ''' + Base class for Schain plotting operations + ''' + + CODE = 'Figure' + colormap = 'jro' + bgcolor = 'white' + __missing = 1E30 + + __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax', + 'zlimits', 'xlabel', 'ylabel', 'xaxis', 'cb_label', 'title', + 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure', + 'showprofile', 'decimation', 'pause'] + + def __init__(self): + + Operation.__init__(self) + self.isConfig = False + self.isPlotConfig = False + + def __fmtTime(self, x, pos): + ''' + ''' + + return '{}'.format(self.getDateTime(x).strftime('%H:%M')) + + def __setup(self, **kwargs): + ''' + Initialize variables + ''' + + self.figures = [] + self.axes = [] + self.cb_axes = [] + self.localtime = kwargs.pop('localtime', True) + self.show = kwargs.get('show', True) + self.save = kwargs.get('save', False) + self.ftp = kwargs.get('ftp', False) + self.colormap = kwargs.get('colormap', self.colormap) + self.colormap_coh = kwargs.get('colormap_coh', 'jet') + self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r') + self.colormaps = kwargs.get('colormaps', None) + self.bgcolor = kwargs.get('bgcolor', self.bgcolor) + self.showprofile = kwargs.get('showprofile', False) + self.title = kwargs.get('wintitle', self.CODE.upper()) + self.cb_label = kwargs.get('cb_label', None) + self.cb_labels = kwargs.get('cb_labels', None) + self.labels = kwargs.get('labels', None) + self.xaxis = kwargs.get('xaxis', 'frequency') + self.zmin = kwargs.get('zmin', None) + self.zmax = kwargs.get('zmax', None) + self.zlimits = kwargs.get('zlimits', None) + self.xmin = kwargs.get('xmin', None) + self.xmax = kwargs.get('xmax', None) + self.xrange = kwargs.get('xrange', 12) + self.xscale = kwargs.get('xscale', None) + self.ymin = kwargs.get('ymin', None) + self.ymax = kwargs.get('ymax', None) + self.yscale = kwargs.get('yscale', None) + self.xlabel = kwargs.get('xlabel', None) + self.decimation = kwargs.get('decimation', None) + self.showSNR = kwargs.get('showSNR', False) + self.oneFigure = kwargs.get('oneFigure', True) + self.width = kwargs.get('width', None) + self.height = kwargs.get('height', None) + self.colorbar = kwargs.get('colorbar', True) + self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1]) + self.channels = kwargs.get('channels', None) + self.titles = kwargs.get('titles', []) + self.polar = False + self.grid = kwargs.get('grid', False) + self.pause = kwargs.get('pause', False) + self.save_labels = kwargs.get('save_labels', None) + self.realtime = kwargs.get('realtime', True) + self.buffering = kwargs.get('buffering', True) + self.throttle = kwargs.get('throttle', 2) + self.exp_code = kwargs.get('exp_code', None) + self.__throttle_plot = apply_throttle(self.throttle) + self.data = PlotterData( + self.CODE, self.throttle, self.exp_code, self.buffering) + + def __setup_plot(self): + ''' + Common setup for all figures, here figures and axes are created + ''' + + self.setup() + + self.time_label = 'LT' if self.localtime else 'UTC' + if self.data.localtime: + self.getDateTime = datetime.datetime.fromtimestamp + else: + self.getDateTime = datetime.datetime.utcfromtimestamp + + if self.width is None: + self.width = 8 + + self.figures = [] + self.axes = [] + self.cb_axes = [] + self.pf_axes = [] + self.cmaps = [] + + size = '15%' if self.ncols == 1 else '30%' + pad = '4%' if self.ncols == 1 else '8%' + + if self.oneFigure: + if self.height is None: + self.height = 1.4 * self.nrows + 1 + fig = plt.figure(figsize=(self.width, self.height), + edgecolor='k', + facecolor='w') + self.figures.append(fig) + for n in range(self.nplots): + ax = fig.add_subplot(self.nrows, self.ncols, + n + 1, polar=self.polar) + ax.tick_params(labelsize=8) + ax.firsttime = True + ax.index = 0 + ax.press = None + self.axes.append(ax) + if self.showprofile: + cax = self.__add_axes(ax, size=size, pad=pad) + cax.tick_params(labelsize=8) + self.pf_axes.append(cax) + else: + if self.height is None: + self.height = 3 + for n in range(self.nplots): + fig = plt.figure(figsize=(self.width, self.height), + edgecolor='k', + facecolor='w') + ax = fig.add_subplot(1, 1, 1, polar=self.polar) + ax.tick_params(labelsize=8) + ax.firsttime = True + ax.index = 0 + ax.press = None + self.figures.append(fig) + self.axes.append(ax) + if self.showprofile: + cax = self.__add_axes(ax, size=size, pad=pad) + cax.tick_params(labelsize=8) + self.pf_axes.append(cax) + + for n in range(self.nrows): + if self.colormaps is not None: + cmap = plt.get_cmap(self.colormaps[n]) + else: + cmap = plt.get_cmap(self.colormap) + cmap.set_bad(self.bgcolor, 1.) + self.cmaps.append(cmap) + + for fig in self.figures: + fig.canvas.mpl_connect('key_press_event', self.OnKeyPress) + fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll) + fig.canvas.mpl_connect('button_press_event', self.onBtnPress) + fig.canvas.mpl_connect('motion_notify_event', self.onMotion) + fig.canvas.mpl_connect('button_release_event', self.onBtnRelease) + if self.show: + fig.show() + + def OnKeyPress(self, event): + ''' + Event for pressing keys (up, down) change colormap + ''' + ax = event.inaxes + if ax in self.axes: + if event.key == 'down': + ax.index += 1 + elif event.key == 'up': + ax.index -= 1 + if ax.index < 0: + ax.index = len(CMAPS) - 1 + elif ax.index == len(CMAPS): + ax.index = 0 + cmap = CMAPS[ax.index] + ax.cbar.set_cmap(cmap) + ax.cbar.draw_all() + ax.plt.set_cmap(cmap) + ax.cbar.patch.figure.canvas.draw() + self.colormap = cmap.name + + def OnBtnScroll(self, event): + ''' + Event for scrolling, scale figure + ''' + cb_ax = event.inaxes + if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]: + ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0] + pt = ax.cbar.ax.bbox.get_points()[:, 1] + nrm = ax.cbar.norm + vmin, vmax, p0, p1, pS = ( + nrm.vmin, nrm.vmax, pt[0], pt[1], event.y) + scale = 2 if event.step == 1 else 0.5 + point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0) + ax.cbar.norm.vmin = point - scale * (point - vmin) + ax.cbar.norm.vmax = point - scale * (point - vmax) + ax.plt.set_norm(ax.cbar.norm) + ax.cbar.draw_all() + ax.cbar.patch.figure.canvas.draw() + + def onBtnPress(self, event): + ''' + Event for mouse button press + ''' + cb_ax = event.inaxes + if cb_ax is None: + return + + if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]: + cb_ax.press = event.x, event.y + else: + cb_ax.press = None + + def onMotion(self, event): + ''' + Event for move inside colorbar + ''' + cb_ax = event.inaxes + if cb_ax is None: + return + if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]: + return + if cb_ax.press is None: + return + + ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0] + xprev, yprev = cb_ax.press + dx = event.x - xprev + dy = event.y - yprev + cb_ax.press = event.x, event.y + scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin + perc = 0.03 + + if event.button == 1: + ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy) + ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy) + elif event.button == 3: + ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy) + ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy) + + ax.cbar.draw_all() + ax.plt.set_norm(ax.cbar.norm) + ax.cbar.patch.figure.canvas.draw() + + def onBtnRelease(self, event): + ''' + Event for mouse button release + ''' + cb_ax = event.inaxes + if cb_ax is not None: + cb_ax.press = None + + def __add_axes(self, ax, size='30%', pad='8%'): + ''' + Add new axes to the given figure + ''' + divider = make_axes_locatable(ax) + nax = divider.new_horizontal(size=size, pad=pad) + ax.figure.add_axes(nax) + return nax + + def setup(self): + ''' + This method should be implemented in the child class, the following + attributes should be set: + + self.nrows: number of rows + 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 + + def fill_gaps(self, x_buffer, y_buffer, z_buffer): + ''' + Create a masked array for missing data + ''' + if x_buffer.shape[0] < 2: + return x_buffer, y_buffer, z_buffer + + deltas = x_buffer[1:] - x_buffer[0:-1] + x_median = numpy.median(deltas) + + index = numpy.where(deltas > 5 * x_median) + + if len(index[0]) != 0: + z_buffer[::, index[0], ::] = self.__missing + z_buffer = numpy.ma.masked_inside(z_buffer, + 0.99 * self.__missing, + 1.01 * self.__missing) + + return x_buffer, y_buffer, z_buffer + + def decimate(self): + + # dx = int(len(self.x)/self.__MAXNUMX) + 1 + dy = int(len(self.y) / self.decimation) + 1 + + # x = self.x[::dx] + x = self.x + y = self.y[::dy] + z = self.z[::, ::, ::dy] + + return x, y, z + + def format(self): + ''' + Set min and max values, labels, ticks and titles + ''' + + if self.xmin is None: + xmin = self.data.min_time + else: + if self.xaxis is 'time': + dt = self.getDateTime(self.data.min_time) + xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) - + datetime.datetime(1970, 1, 1)).total_seconds() + if self.data.localtime: + xmin += time.timezone + else: + xmin = self.xmin + + if self.xmax is None: + xmax = xmin + self.xrange * 60 * 60 + else: + if self.xaxis is 'time': + dt = self.getDateTime(self.data.max_time) + xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) - + datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds() + if self.data.localtime: + xmax += time.timezone + else: + xmax = self.xmax + + ymin = self.ymin if self.ymin else numpy.nanmin(self.y) + ymax = self.ymax if self.ymax else numpy.nanmax(self.y) + + Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]) + i = 1 if numpy.where( + abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0] + ystep = Y[i] / 10. + + if self.xaxis is not 'time': + X = numpy.array([1, 2, 5, 10, 20, 50, 100, + 200, 500, 1000, 2000, 5000])/2. + i = 1 if numpy.where( + abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0] + xstep = X[i] / 10. + + for n, ax in enumerate(self.axes): + if ax.firsttime: + ax.set_facecolor(self.bgcolor) + ax.yaxis.set_major_locator(MultipleLocator(ystep)) + if self.xscale: + ax.xaxis.set_major_formatter(FuncFormatter( + lambda x, pos: '{0:g}'.format(x*self.xscale))) + if self.xscale: + ax.yaxis.set_major_formatter(FuncFormatter( + lambda x, pos: '{0:g}'.format(x*self.yscale))) + if self.xaxis is 'time': + ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime)) + ax.xaxis.set_major_locator(LinearLocator(9)) + else: + ax.xaxis.set_major_locator(MultipleLocator(xstep)) + if self.xlabel is not None: + ax.set_xlabel(self.xlabel) + ax.set_ylabel(self.ylabel) + ax.firsttime = False + if self.showprofile: + self.pf_axes[n].set_ylim(ymin, ymax) + self.pf_axes[n].set_xlim(self.zmin, self.zmax) + self.pf_axes[n].set_xlabel('dB') + self.pf_axes[n].grid(b=True, axis='x') + [tick.set_visible(False) + for tick in self.pf_axes[n].get_yticklabels()] + if self.colorbar: + ax.cbar = plt.colorbar( + ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10) + ax.cbar.ax.tick_params(labelsize=8) + ax.cbar.ax.press = None + if self.cb_label: + ax.cbar.set_label(self.cb_label, size=8) + elif self.cb_labels: + ax.cbar.set_label(self.cb_labels[n], size=8) + else: + ax.cbar = None + if self.grid: + ax.grid(True) + + if not self.polar: + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + ax.set_title('{} {} {}'.format( + self.titles[n], + self.getDateTime(self.data.max_time).strftime( + '%Y-%m-%dT%H:%M:%S'), + self.time_label), + size=8) + else: + ax.set_title('{}'.format(self.titles[n]), size=8) + ax.set_ylim(0, 90) + ax.set_yticks(numpy.arange(0, 90, 20)) + ax.yaxis.labelpad = 40 + + def clear_figures(self): + ''' + Reset axes for redraw plots + ''' + + for ax in self.axes: + ax.clear() + ax.firsttime = True + if ax.cbar: + ax.cbar.remove() + + def __plot(self): + ''' + Main function to plot, format and save figures + ''' + + #try: + self.plot() + self.format() + #except Exception as e: + # log.warning('{} Plot could not be updated... check data'.format( + # self.CODE), self.name) + # log.error(str(e), '') + # return + + 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.tight_layout() + 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: + + if self.save_labels: + labels = self.save_labels + else: + labels = list(range(self.nrows)) + + if self.oneFigure: + label = '' + else: + label = '-{}'.format(labels[n]) + figname = os.path.join( + self.save, + self.CODE, + '{}{}_{}.png'.format( + self.CODE, + label, + self.getDateTime(self.data.max_time).strftime( + '%Y%m%d_%H%M%S'), + ) + ) + log.log('Saving figure: {}'.format(figname), self.name) + if not os.path.isdir(os.path.dirname(figname)): + os.makedirs(os.path.dirname(figname)) + fig.savefig(figname) + + def plot(self): + ''' + Must be defined in the child class + ''' + raise NotImplementedError + + def run(self, dataOut, **kwargs): + + if dataOut.flagNoData and not dataOut.error: + return dataOut + + if dataOut.error: + coerce = True + else: + coerce = False + + if self.isConfig is False: + self.__setup(**kwargs) + self.data.setup() + self.isConfig = True + + if dataOut.type == 'Parameters': + tm = dataOut.utctimeInit + else: + tm = dataOut.utctime + + if dataOut.useLocalTime: + if not self.localtime: + tm += time.timezone + else: + if self.localtime: + tm -= time.timezone + + if self.data and (tm - self.data.min_time) >= self.xrange*60*60: + self.__plot() + self.data.setup() + self.clear_figures() + + self.data.update(dataOut, tm) + + if self.isPlotConfig is False: + self.__setup_plot() + self.isPlotConfig = True + + if self.realtime: + self.__plot() + else: + self.__throttle_plot(self.__plot, coerce=coerce) + + figpause(0.001) + + def close(self): + + if self.data and self.pause: + figpause(10) + diff --git a/schainpy/model/graphics/jroplot_correlation.py b/schainpy/model/graphics/jroplot_correlation.py index 9dea381..37fdc5c 100644 --- a/schainpy/model/graphics/jroplot_correlation.py +++ b/schainpy/model/graphics/jroplot_correlation.py @@ -5,7 +5,7 @@ import copy from schainpy.model import * from .figure import Figure, isRealtime -class CorrelationPlot(Figure): +class CorrelationPlot_(Figure): isConfig = None __nsubplots = None diff --git a/schainpy/model/graphics/jroplot_data.py b/schainpy/model/graphics/jroplot_data.py index f5d5bd9..e85be46 100644 --- a/schainpy/model/graphics/jroplot_data.py +++ b/schainpy/model/graphics/jroplot_data.py @@ -1,41 +1,32 @@ +''' +New Plots Operations + +@author: juan.espinoza@jro.igp.gob.pe +''' + -import os import time -import glob import datetime -from multiprocessing import Process - -import zmq import numpy -import matplotlib -import matplotlib.pyplot as plt -from matplotlib.patches import Polygon -from mpl_toolkits.axes_grid1 import make_axes_locatable -from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator -from schainpy.model.proc.jroproc_base import Operation +from schainpy.model.graphics.jroplot_base import Plot, plt from schainpy.utils import log -jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90] -blu_values = matplotlib.pyplot.get_cmap( - 'seismic_r', 20)(numpy.arange(20))[10:15] -ncmap = matplotlib.colors.LinearSegmentedColormap.from_list( - 'jro', numpy.vstack((blu_values, jet_values))) -matplotlib.pyplot.register_cmap(cmap=ncmap) - -CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis', 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')] - EARTH_RADIUS = 6.3710e3 + def ll2xy(lat1, lon1, lat2, lon2): p = 0.017453292519943295 - a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2 + a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \ + numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2 r = 12742 * numpy.arcsin(numpy.sqrt(a)) - theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)*numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p)) + theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p) + * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p)) theta = -theta + numpy.pi/2 return r*numpy.cos(theta), r*numpy.sin(theta) + def km2deg(km): ''' Convert distance in km to degrees @@ -43,536 +34,8 @@ def km2deg(km): return numpy.rad2deg(km/EARTH_RADIUS) -def figpause(interval): - backend = plt.rcParams['backend'] - if backend in matplotlib.rcsetup.interactive_bk: - figManager = matplotlib._pylab_helpers.Gcf.get_active() - if figManager is not None: - canvas = figManager.canvas - if canvas.figure.stale: - canvas.draw() - try: - canvas.start_event_loop(interval) - except: - pass - return - -def popup(message): - ''' - ''' - - fig = plt.figure(figsize=(12, 8), facecolor='r') - text = '\n'.join([s.strip() for s in message.split(':')]) - fig.text(0.01, 0.5, text, ha='left', va='center', size='20', weight='heavy', color='w') - fig.show() - figpause(1000) - - -class PlotData(Operation, Process): - ''' - Base class for Schain plotting operations - ''' - - CODE = 'Figure' - colormap = 'jro' - bgcolor = 'white' - CONFLATE = False - __missing = 1E30 - - __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax', - 'zlimits', 'xlabel', 'ylabel', 'xaxis','cb_label', 'title', - 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure', - 'showprofile', 'decimation', 'ftp'] - - def __init__(self, **kwargs): - - Operation.__init__(self, plot=True, **kwargs) - Process.__init__(self) - - self.kwargs['code'] = self.CODE - self.mp = False - self.data = None - self.isConfig = False - self.figures = [] - self.axes = [] - self.cb_axes = [] - self.localtime = kwargs.pop('localtime', True) - self.show = kwargs.get('show', True) - self.save = kwargs.get('save', False) - self.ftp = kwargs.get('ftp', False) - self.colormap = kwargs.get('colormap', self.colormap) - self.colormap_coh = kwargs.get('colormap_coh', 'jet') - self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r') - self.colormaps = kwargs.get('colormaps', None) - self.bgcolor = kwargs.get('bgcolor', self.bgcolor) - self.showprofile = kwargs.get('showprofile', False) - self.title = kwargs.get('wintitle', self.CODE.upper()) - self.cb_label = kwargs.get('cb_label', None) - self.cb_labels = kwargs.get('cb_labels', None) - self.labels = kwargs.get('labels', None) - self.xaxis = kwargs.get('xaxis', 'frequency') - self.zmin = kwargs.get('zmin', None) - self.zmax = kwargs.get('zmax', None) - self.zlimits = kwargs.get('zlimits', None) - self.xmin = kwargs.get('xmin', None) - self.xmax = kwargs.get('xmax', None) - self.xrange = kwargs.get('xrange', 24) - self.xscale = kwargs.get('xscale', None) - self.ymin = kwargs.get('ymin', None) - self.ymax = kwargs.get('ymax', None) - self.yscale = kwargs.get('yscale', None) - self.xlabel = kwargs.get('xlabel', None) - self.decimation = kwargs.get('decimation', None) - self.showSNR = kwargs.get('showSNR', False) - self.oneFigure = kwargs.get('oneFigure', True) - self.width = kwargs.get('width', None) - self.height = kwargs.get('height', None) - self.colorbar = kwargs.get('colorbar', True) - self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1]) - self.channels = kwargs.get('channels', None) - self.titles = kwargs.get('titles', []) - self.polar = False - self.grid = kwargs.get('grid', False) - - def __fmtTime(self, x, pos): - ''' - ''' - - return '{}'.format(self.getDateTime(x).strftime('%H:%M')) - - def __setup(self): - ''' - Common setup for all figures, here figures and axes are created - ''' - - if self.CODE not in self.data: - raise ValueError(log.error('Missing data for {}'.format(self.CODE), - self.name)) - - self.setup() - - self.time_label = 'LT' if self.localtime else 'UTC' - if self.data.localtime: - self.getDateTime = datetime.datetime.fromtimestamp - else: - self.getDateTime = datetime.datetime.utcfromtimestamp - - if self.width is None: - self.width = 8 - - self.figures = [] - self.axes = [] - self.cb_axes = [] - self.pf_axes = [] - self.cmaps = [] - - size = '15%' if self.ncols == 1 else '30%' - pad = '4%' if self.ncols == 1 else '8%' - - if self.oneFigure: - if self.height is None: - self.height = 1.4 * self.nrows + 1 - fig = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - self.figures.append(fig) - for n in range(self.nplots): - ax = fig.add_subplot(self.nrows, self.ncols, - n + 1, polar=self.polar) - ax.tick_params(labelsize=8) - ax.firsttime = True - ax.index = 0 - ax.press = None - self.axes.append(ax) - if self.showprofile: - cax = self.__add_axes(ax, size=size, pad=pad) - cax.tick_params(labelsize=8) - self.pf_axes.append(cax) - else: - if self.height is None: - self.height = 3 - for n in range(self.nplots): - fig = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - ax = fig.add_subplot(1, 1, 1, polar=self.polar) - ax.tick_params(labelsize=8) - ax.firsttime = True - ax.index = 0 - ax.press = None - self.figures.append(fig) - self.axes.append(ax) - if self.showprofile: - cax = self.__add_axes(ax, size=size, pad=pad) - cax.tick_params(labelsize=8) - self.pf_axes.append(cax) - - for n in range(self.nrows): - if self.colormaps is not None: - cmap = plt.get_cmap(self.colormaps[n]) - else: - cmap = plt.get_cmap(self.colormap) - cmap.set_bad(self.bgcolor, 1.) - self.cmaps.append(cmap) - - for fig in self.figures: - fig.canvas.mpl_connect('key_press_event', self.OnKeyPress) - fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll) - fig.canvas.mpl_connect('button_press_event', self.onBtnPress) - fig.canvas.mpl_connect('motion_notify_event', self.onMotion) - fig.canvas.mpl_connect('button_release_event', self.onBtnRelease) - if self.show: - fig.show() - - def OnKeyPress(self, event): - ''' - Event for pressing keys (up, down) change colormap - ''' - ax = event.inaxes - if ax in self.axes: - if event.key == 'down': - ax.index += 1 - elif event.key == 'up': - ax.index -= 1 - if ax.index < 0: - ax.index = len(CMAPS) - 1 - elif ax.index == len(CMAPS): - ax.index = 0 - cmap = CMAPS[ax.index] - ax.cbar.set_cmap(cmap) - ax.cbar.draw_all() - ax.plt.set_cmap(cmap) - ax.cbar.patch.figure.canvas.draw() - self.colormap = cmap.name - - def OnBtnScroll(self, event): - ''' - Event for scrolling, scale figure - ''' - cb_ax = event.inaxes - if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]: - ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0] - pt = ax.cbar.ax.bbox.get_points()[:, 1] - nrm = ax.cbar.norm - vmin, vmax, p0, p1, pS = ( - nrm.vmin, nrm.vmax, pt[0], pt[1], event.y) - scale = 2 if event.step == 1 else 0.5 - point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0) - ax.cbar.norm.vmin = point - scale * (point - vmin) - ax.cbar.norm.vmax = point - scale * (point - vmax) - ax.plt.set_norm(ax.cbar.norm) - ax.cbar.draw_all() - ax.cbar.patch.figure.canvas.draw() - - def onBtnPress(self, event): - ''' - Event for mouse button press - ''' - cb_ax = event.inaxes - if cb_ax is None: - return - - if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]: - cb_ax.press = event.x, event.y - else: - cb_ax.press = None - - def onMotion(self, event): - ''' - Event for move inside colorbar - ''' - cb_ax = event.inaxes - if cb_ax is None: - return - if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]: - return - if cb_ax.press is None: - return - - ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0] - xprev, yprev = cb_ax.press - dx = event.x - xprev - dy = event.y - yprev - cb_ax.press = event.x, event.y - scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin - perc = 0.03 - - if event.button == 1: - ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy) - ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy) - elif event.button == 3: - ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy) - ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy) - - ax.cbar.draw_all() - ax.plt.set_norm(ax.cbar.norm) - ax.cbar.patch.figure.canvas.draw() - - def onBtnRelease(self, event): - ''' - Event for mouse button release - ''' - cb_ax = event.inaxes - if cb_ax is not None: - cb_ax.press = None - - def __add_axes(self, ax, size='30%', pad='8%'): - ''' - Add new axes to the given figure - ''' - divider = make_axes_locatable(ax) - nax = divider.new_horizontal(size=size, pad=pad) - ax.figure.add_axes(nax) - return nax - - self.setup() - - def setup(self): - ''' - This method should be implemented in the child class, the following - attributes should be set: - - self.nrows: number of rows - 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 - - def fill_gaps(self, x_buffer, y_buffer, z_buffer): - ''' - Create a masked array for missing data - ''' - if x_buffer.shape[0] < 2: - return x_buffer, y_buffer, z_buffer - - deltas = x_buffer[1:] - x_buffer[0:-1] - x_median = numpy.median(deltas) - - index = numpy.where(deltas > 5 * x_median) - - if len(index[0]) != 0: - z_buffer[::, index[0], ::] = self.__missing - z_buffer = numpy.ma.masked_inside(z_buffer, - 0.99 * self.__missing, - 1.01 * self.__missing) - - return x_buffer, y_buffer, z_buffer - - def decimate(self): - - # dx = int(len(self.x)/self.__MAXNUMX) + 1 - dy = int(len(self.y) / self.decimation) + 1 - - # x = self.x[::dx] - x = self.x - y = self.y[::dy] - z = self.z[::, ::, ::dy] - - return x, y, z - def format(self): - ''' - Set min and max values, labels, ticks and titles - ''' - - if self.xmin is None: - xmin = self.min_time - else: - if self.xaxis is 'time': - dt = self.getDateTime(self.min_time) - xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) - - datetime.datetime(1970, 1, 1)).total_seconds() - if self.data.localtime: - xmin += time.timezone - else: - xmin = self.xmin - - if self.xmax is None: - xmax = xmin + self.xrange * 60 * 60 - else: - if self.xaxis is 'time': - dt = self.getDateTime(self.max_time) - xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) - - datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds() - if self.data.localtime: - xmax += time.timezone - else: - xmax = self.xmax - - ymin = self.ymin if self.ymin else numpy.nanmin(self.y) - ymax = self.ymax if self.ymax else numpy.nanmax(self.y) - - Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]) - i = 1 if numpy.where(abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0] - ystep = Y[i] / 10. - - if self.xaxis is not 'time': - X = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])/2. - i = 1 if numpy.where(abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0] - xstep = X[i] / 10. - - for n, ax in enumerate(self.axes): - if ax.firsttime: - ax.set_facecolor(self.bgcolor) - ax.yaxis.set_major_locator(MultipleLocator(ystep)) - if self.xscale: - ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{0:g}'.format(x*self.xscale))) - if self.xscale: - ax.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{0:g}'.format(x*self.yscale))) - if self.xaxis is 'time': - ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime)) - ax.xaxis.set_major_locator(LinearLocator(9)) - else: - ax.xaxis.set_major_locator(MultipleLocator(xstep)) - if self.xlabel is not None: - ax.set_xlabel(self.xlabel) - ax.set_ylabel(self.ylabel) - ax.firsttime = False - if self.showprofile: - self.pf_axes[n].set_ylim(ymin, ymax) - self.pf_axes[n].set_xlim(self.zmin, self.zmax) - self.pf_axes[n].set_xlabel('dB') - self.pf_axes[n].grid(b=True, axis='x') - [tick.set_visible(False) - for tick in self.pf_axes[n].get_yticklabels()] - if self.colorbar: - ax.cbar = plt.colorbar( - ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10) - ax.cbar.ax.tick_params(labelsize=8) - ax.cbar.ax.press = None - if self.cb_label: - ax.cbar.set_label(self.cb_label, size=8) - elif self.cb_labels: - ax.cbar.set_label(self.cb_labels[n], size=8) - else: - ax.cbar = None - if self.grid: - ax.grid(True) - - if not self.polar: - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) - ax.set_title('{} {} {}'.format( - self.titles[n], - self.getDateTime(self.max_time).strftime('%Y-%m-%dT%H:%M:%S'), - self.time_label), - size=8) - else: - ax.set_title('{}'.format(self.titles[n]), size=8) - ax.set_ylim(0, 90) - ax.set_yticks(numpy.arange(0, 90, 20)) - ax.yaxis.labelpad = 40 - - def __plot(self): - ''' - ''' - log.log('Plotting', self.name) - - try: - self.plot() - self.format() - except Exception as e: - log.warning('{} Plot could not be updated... check data'.format(self.CODE), self.name) - log.error(str(e), '') - return - - 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.tight_layout() - fig.canvas.manager.set_window_title('{} - {}'.format(self.title, - self.getDateTime(self.max_time).strftime('%Y/%m/%d'))) - fig.canvas.draw() - - if self.save and (self.data.ended or not self.data.buffering): - - if self.save_labels: - labels = self.save_labels - else: - labels = list(range(self.nrows)) - - if self.oneFigure: - label = '' - else: - label = '-{}'.format(labels[n]) - figname = os.path.join( - self.save, - self.CODE, - '{}{}_{}.png'.format( - self.CODE, - label, - self.getDateTime(self.saveTime).strftime( - '%Y%m%d_%H%M%S'), - ) - ) - log.log('Saving figure: {}'.format(figname), self.name) - if not os.path.isdir(os.path.dirname(figname)): - os.makedirs(os.path.dirname(figname)) - fig.savefig(figname) - - def plot(self): - ''' - ''' - raise NotImplementedError - - def run(self): - - log.log('Starting', self.name) - - context = zmq.Context() - receiver = context.socket(zmq.SUB) - receiver.setsockopt(zmq.SUBSCRIBE, '') - receiver.setsockopt(zmq.CONFLATE, self.CONFLATE) - - if 'server' in self.kwargs['parent']: - receiver.connect( - 'ipc:///tmp/{}.plots'.format(self.kwargs['parent']['server'])) - else: - receiver.connect("ipc:///tmp/zmq.plots") - - while True: - try: - self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK) - if self.data.localtime and self.localtime: - self.times = self.data.times - elif self.data.localtime and not self.localtime: - self.times = self.data.times + time.timezone - elif not self.data.localtime and self.localtime: - self.times = self.data.times - time.timezone - else: - self.times = self.data.times - - self.min_time = self.times[0] - self.max_time = self.times[-1] - - if self.isConfig is False: - self.__setup() - self.isConfig = True - - self.__plot() - - except zmq.Again as e: - if self.data and self.data.ended: - break - log.log('Waiting for data...') - if self.data: - figpause(self.data.throttle) - else: - time.sleep(2) - - def close(self): - if self.data: - self.__plot() - - -class PlotSpectraData(PlotData): +class SpectraPlot(Plot): ''' Plot for Spectra data ''' @@ -644,10 +107,9 @@ class PlotSpectraData(PlotData): ax.plt_mean.set_data(mean, y) self.titles.append('CH {}: {:3.2f}dB'.format(n, noise)) - self.saveTime = self.max_time -class PlotCrossSpectraData(PlotData): +class CrossSpectraPlot(Plot): CODE = 'cspc' zmin_coh = None @@ -741,10 +203,8 @@ class PlotCrossSpectraData(PlotData): ax.plt.set_array(phase.T.ravel()) self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1])) - self.saveTime = self.max_time - -class PlotSpectraMeanData(PlotSpectraData): +class SpectraMeanPlot(SpectraPlot): ''' Plot for Spectra and Mean ''' @@ -752,7 +212,7 @@ class PlotSpectraMeanData(PlotSpectraData): colormap = 'jro' -class PlotRTIData(PlotData): +class RTIPlot(Plot): ''' Plot for RTI data ''' @@ -771,7 +231,7 @@ class PlotRTIData(PlotData): self.CODE.upper(), x) for x in range(self.nrows)] def plot(self): - self.x = self.times + self.x = self.data.times self.y = self.data.heights self.z = self.data[self.CODE] self.z = numpy.ma.masked_invalid(self.z) @@ -781,7 +241,7 @@ class PlotRTIData(PlotData): else: x, y, z = self.fill_gaps(*self.decimate()) - for n, ax in enumerate(self.axes): + for n, ax in enumerate(self.axes): self.zmin = self.zmin if self.zmin else numpy.min(self.z) self.zmax = self.zmax if self.zmax else numpy.max(self.z) if ax.firsttime: @@ -807,10 +267,8 @@ class PlotRTIData(PlotData): ax.plot_noise.set_data(numpy.repeat( self.data['noise'][n][-1], len(self.y)), self.y) - self.saveTime = self.min_time - -class PlotCOHData(PlotRTIData): +class CoherencePlot(RTIPlot): ''' Plot for Coherence data ''' @@ -833,7 +291,7 @@ class PlotCOHData(PlotRTIData): 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs] -class PlotPHASEData(PlotCOHData): +class PhasePlot(CoherencePlot): ''' Plot for Phase map data ''' @@ -842,7 +300,7 @@ class PlotPHASEData(PlotCOHData): colormap = 'seismic' -class PlotNoiseData(PlotData): +class NoisePlot(Plot): ''' Plot for noise ''' @@ -860,8 +318,8 @@ class PlotNoiseData(PlotData): def plot(self): - x = self.times - xmin = self.min_time + x = self.data.times + xmin = self.data.min_time xmax = xmin + self.xrange * 60 * 60 Y = self.data[self.CODE] @@ -877,10 +335,9 @@ class PlotNoiseData(PlotData): self.ymin = numpy.nanmin(Y) - 5 self.ymax = numpy.nanmax(Y) + 5 - self.saveTime = self.min_time -class PlotSNRData(PlotRTIData): +class SnrPlot(RTIPlot): ''' Plot for SNR Data ''' @@ -889,7 +346,7 @@ class PlotSNRData(PlotRTIData): colormap = 'jet' -class PlotDOPData(PlotRTIData): +class DopplerPlot(RTIPlot): ''' Plot for DOPPLER Data ''' @@ -898,7 +355,7 @@ class PlotDOPData(PlotRTIData): colormap = 'jet' -class PlotSkyMapData(PlotData): +class SkyMapPlot(Plot): ''' Plot for meteors detection data ''' @@ -938,16 +395,15 @@ class PlotSkyMapData(PlotData): else: ax.plot.set_data(x, y) - dt1 = self.getDateTime(self.min_time).strftime('%y/%m/%d %H:%M:%S') - dt2 = self.getDateTime(self.max_time).strftime('%y/%m/%d %H:%M:%S') + dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S') + dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S') title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1, dt2, len(x)) self.titles[0] = title - self.saveTime = self.max_time -class PlotParamData(PlotRTIData): +class ParametersPlot(RTIPlot): ''' Plot for data_param object ''' @@ -973,7 +429,7 @@ class PlotParamData(PlotRTIData): def plot(self): self.data.normalize_heights() - self.x = self.times + self.x = self.data.times self.y = self.data.heights if self.showSNR: self.z = numpy.concatenate( @@ -990,7 +446,7 @@ class PlotParamData(PlotRTIData): x, y, z = self.fill_gaps(*self.decimate()) for n, ax in enumerate(self.axes): - + self.zmax = self.zmax if self.zmax is not None else numpy.max( self.z[n]) self.zmin = self.zmin if self.zmin is not None else numpy.min( @@ -1015,10 +471,8 @@ class PlotParamData(PlotRTIData): cmap=self.cmaps[n] ) - self.saveTime = self.min_time - -class PlotOutputData(PlotParamData): +class OutputPlot(ParametersPlot): ''' Plot data_output object ''' @@ -1027,9 +481,9 @@ class PlotOutputData(PlotParamData): colormap = 'seismic' -class PlotPolarMapData(PlotData): +class PolarMapPlot(Plot): ''' - Plot for meteors detection data + Plot for weather radar ''' CODE = 'param' @@ -1058,20 +512,24 @@ class PlotPolarMapData(PlotData): self.cb_labels = self.data.meta['units'] self.lat = self.data.meta['latitude'] self.lon = self.data.meta['longitude'] - self.xmin, self.xmax = float(km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon) - self.ymin, self.ymax = float(km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat) + self.xmin, self.xmax = float( + km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon) + self.ymin, self.ymax = float( + km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat) # self.polar = True - def plot(self): - + def plot(self): + for n, ax in enumerate(self.axes): data = self.data['param'][self.channels[n]] - - zeniths = numpy.linspace(0, self.data.meta['max_range'], data.shape[1]) - if self.mode == 'E': + + zeniths = numpy.linspace( + 0, self.data.meta['max_range'], data.shape[1]) + if self.mode == 'E': azimuths = -numpy.radians(self.data.heights)+numpy.pi/2 r, theta = numpy.meshgrid(zeniths, azimuths) - x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])) + x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin( + theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])) x = km2deg(x) + self.lon y = km2deg(y) + self.lat else: @@ -1083,35 +541,36 @@ class PlotPolarMapData(PlotData): if ax.firsttime: if self.zlimits is not None: self.zmin, self.zmax = self.zlimits[n] - ax.plt = ax.pcolormesh(#r, theta, numpy.ma.array(data, mask=numpy.isnan(data)), - x, y, numpy.ma.array(data, mask=numpy.isnan(data)), - vmin=self.zmin, - vmax=self.zmax, - cmap=self.cmaps[n]) + ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)), + x, y, numpy.ma.array(data, mask=numpy.isnan(data)), + vmin=self.zmin, + vmax=self.zmax, + cmap=self.cmaps[n]) else: if self.zlimits is not None: self.zmin, self.zmax = self.zlimits[n] ax.collections.remove(ax.collections[0]) - ax.plt = ax.pcolormesh(# r, theta, numpy.ma.array(data, mask=numpy.isnan(data)), - x, y, numpy.ma.array(data, mask=numpy.isnan(data)), - vmin=self.zmin, - vmax=self.zmax, - cmap=self.cmaps[n]) + ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)), + x, y, numpy.ma.array(data, mask=numpy.isnan(data)), + vmin=self.zmin, + vmax=self.zmax, + cmap=self.cmaps[n]) if self.mode == 'A': continue - + # plot district names f = open('/data/workspace/schain_scripts/distrito.csv') for line in f: label, lon, lat = [s.strip() for s in line.split(',') if s] lat = float(lat) - lon = float(lon) + lon = float(lon) # ax.plot(lon, lat, '.b', ms=2) - ax.text(lon, lat, label.decode('utf8'), ha='center', va='bottom', size='8', color='black') - + ax.text(lon, lat, label.decode('utf8'), ha='center', + va='bottom', size='8', color='black') + # plot limites - limites =[] + limites = [] tmp = [] for line in open('/data/workspace/schain_scripts/lima.csv'): if '#' in line: @@ -1122,7 +581,8 @@ class PlotPolarMapData(PlotData): values = line.strip().split(',') tmp.append((float(values[0]), float(values[1]))) for points in limites: - ax.add_patch(Polygon(points, ec='k', fc='none', ls='--', lw=0.5)) + ax.add_patch( + Polygon(points, ec='k', fc='none', ls='--', lw=0.5)) # plot Cuencas for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'): @@ -1133,22 +593,21 @@ class PlotPolarMapData(PlotData): # plot grid for r in (15, 30, 45, 60): - ax.add_artist(plt.Circle((self.lon, self.lat), km2deg(r), color='0.6', fill=False, lw=0.2)) + ax.add_artist(plt.Circle((self.lon, self.lat), + km2deg(r), color='0.6', fill=False, lw=0.2)) ax.text( self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180), self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180), - '{}km'.format(r), + '{}km'.format(r), ha='center', va='bottom', size='8', color='0.6', weight='heavy') - + if self.mode == 'E': title = 'El={}$^\circ$'.format(self.data.meta['elevation']) label = 'E{:02d}'.format(int(self.data.meta['elevation'])) else: title = 'Az={}$^\circ$'.format(self.data.meta['azimuth']) label = 'A{:02d}'.format(int(self.data.meta['azimuth'])) - - self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels] - self.titles = ['{} {}'.format(self.data.parameters[x], title) for x in self.channels] - self.saveTime = self.max_time - \ No newline at end of file + self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels] + self.titles = ['{} {}'.format( + self.data.parameters[x], title) for x in self.channels] diff --git a/schainpy/model/graphics/jroplot_heispectra.py b/schainpy/model/graphics/jroplot_heispectra.py index bfc6e18..2bca309 100644 --- a/schainpy/model/graphics/jroplot_heispectra.py +++ b/schainpy/model/graphics/jroplot_heispectra.py @@ -10,7 +10,7 @@ import numpy from .figure import Figure, isRealtime from .plotting_codes import * -class SpectraHeisScope(Figure): +class SpectraHeisScope_(Figure): isConfig = None @@ -173,7 +173,7 @@ class SpectraHeisScope(Figure): wr_period=wr_period, thisDatetime=thisDatetime) -class RTIfromSpectraHeis(Figure): +class RTIfromSpectraHeis_(Figure): isConfig = None __nsubplots = None diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index f3259ac..b6d0046 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -7,14 +7,190 @@ from .plotting_codes import * from schainpy.model.proc.jroproc_base import MPDecorator from schainpy.utils import log -class FitGauPlot(Figure): +class ParamLine_(Figure): + + isConfig = None + + def __init__(self): + + self.isConfig = False + self.WIDTH = 300 + self.HEIGHT = 200 + self.counter_imagwr = 0 + + def getSubplots(self): + + nrow = self.nplots + ncol = 3 + return nrow, ncol + + def setup(self, id, nplots, wintitle, show): + + self.nplots = nplots + + self.createFigure(id=id, + wintitle=wintitle, + show=show) + + nrow,ncol = self.getSubplots() + colspan = 3 + rowspan = 1 + + for i in range(nplots): + self.addAxes(nrow, ncol, i, 0, colspan, rowspan) + + def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + yreal = y[channelIndexList,:].real + yimag = y[channelIndexList,:].imag + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity - IQ" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag)) + if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag)) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + + axes.pline(x, yreal[i,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2) + + def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) + yreal = y.real + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(yreal) + if ymax == None: ymax = numpy.nanmax(yreal) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = yreal[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, save=False, + figpath='./', figfile=None, show=True, wr_period=1, + ftp=False, server=None, folder=None, username=None, password=None): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + """ + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError("Channel %d is not in dataOut.channelList" % channel) + channelIndexList.append(dataOut.channelList.index(channel)) + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) + + y = dataOut.RR + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = y[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + self.draw() + + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex) + figfile = self.getFilename(name = str_datetime) + + self.save(figpath=figpath, + figfile=figfile, + save=save, + ftp=ftp, + wr_period=wr_period, + thisDatetime=thisDatetime) + + + +class SpcParamPlot_(Figure): isConfig = None __nsubplots = None WIDTHPROF = None HEIGHTPROF = None - PREFIX = 'fitgau' + PREFIX = 'SpcParam' def __init__(self, **kwargs): Figure.__init__(self, **kwargs) @@ -83,7 +259,7 @@ class FitGauPlot(Figure): 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, realtime=False, - xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1): + xaxis="frequency", colormap='jet', normFactor=None , Selector = 0): """ @@ -119,23 +295,22 @@ class FitGauPlot(Figure): # else: # factor = normFactor if xaxis == "frequency": - x = dataOut.spc_range[0] + x = dataOut.spcparam_range[0] xlabel = "Frequency (kHz)" elif xaxis == "time": - x = dataOut.spc_range[1] + x = dataOut.spcparam_range[1] xlabel = "Time (ms)" - else: - x = dataOut.spc_range[2] + else: + x = dataOut.spcparam_range[2] xlabel = "Velocity (m/s)" - ylabel = "Range (Km)" + ylabel = "Range (km)" y = dataOut.getHeiRange() - z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor - print('GausSPC', z[0,32,10:40]) + z = dataOut.SPCparam[Selector] /1966080.0#/ dataOut.normFactor#GauSelector] #dataOut.data_spc/factor z = numpy.where(numpy.isfinite(z), z, numpy.NAN) zdB = 10*numpy.log10(z) @@ -218,7 +393,7 @@ class FitGauPlot(Figure): -class MomentsPlot(Figure): +class MomentsPlot_(Figure): isConfig = None __nsubplots = None @@ -405,7 +580,7 @@ class MomentsPlot(Figure): thisDatetime=thisDatetime) -class SkyMapPlot(Figure): +class SkyMapPlot_(Figure): __isConfig = None __nsubplots = None @@ -561,7 +736,7 @@ class SkyMapPlot(Figure): -class WindProfilerPlot(Figure): +class WindProfilerPlot_(Figure): __isConfig = None __nsubplots = None @@ -657,7 +832,7 @@ class WindProfilerPlot(Figure): # tmax = None x = dataOut.getTimeRange1(dataOut.paramInterval) - y = dataOut.heightList + y = dataOut.heightList z = dataOut.data_output.copy() nplots = z.shape[0] #Number of wind dimensions estimated nplotsw = nplots @@ -666,13 +841,14 @@ class WindProfilerPlot(Figure): #If there is a SNR function defined if dataOut.data_SNR is not None: nplots += 1 - SNR = dataOut.data_SNR - SNRavg = numpy.average(SNR, axis=0) + SNR = dataOut.data_SNR[0] + SNRavg = SNR#numpy.average(SNR, axis=0) SNRdB = 10*numpy.log10(SNR) SNRavgdB = 10*numpy.log10(SNRavg) - if SNRthresh == None: SNRthresh = -5.0 + if SNRthresh == None: + SNRthresh = -5.0 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0] for i in range(nplotsw): @@ -741,8 +917,7 @@ class WindProfilerPlot(Figure): axes = self.axesList[i*self.__nsubplots] z1 = z[i,:].reshape((1,-1))*windFactor[i] - #z1=numpy.ma.masked_where(z1==0.,z1) - + axes.pcolorbuffer(x, y, z1, xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i], xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, @@ -774,7 +949,7 @@ class WindProfilerPlot(Figure): update_figfile = True @MPDecorator -class ParametersPlot(Figure): +class ParametersPlot_(Figure): __isConfig = None __nsubplots = None @@ -792,8 +967,8 @@ class ParametersPlot(Figure): self.isConfig = False self.__nsubplots = 1 - self.WIDTH = 800 - self.HEIGHT = 180 + self.WIDTH = 300 + self.HEIGHT = 550 self.WIDTHPROF = 120 self.HEIGHTPROF = 0 self.counter_imagwr = 0 @@ -905,7 +1080,7 @@ class ParametersPlot(Figure): # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y")) xlabel = "" - ylabel = "Range (Km)" + ylabel = "Range (km)" update_figfile = False @@ -949,24 +1124,81 @@ class ParametersPlot(Figure): self.setWinTitle(title) - for i in range(self.nchan): - index = channelIndexList[i] - title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[i*self.plotFact] - z1 = z[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, z1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, - ticksize=9, cblabel='', cbsize="1%",colormap=colormap) - - if showSNR: - title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[i*self.plotFact + 1] - SNRdB1 = SNRdB[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, SNRdB1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, - ticksize=9, cblabel='', cbsize="1%",colormap='jet') +# for i in range(self.nchan): +# index = channelIndexList[i] +# title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) +# axes = self.axesList[i*self.plotFact] +# z1 = z[i,:].reshape((1,-1)) +# axes.pcolorbuffer(x, y, z1, +# xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, +# xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, +# ticksize=9, cblabel='', cbsize="1%",colormap=colormap) +# +# if showSNR: +# title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) +# axes = self.axesList[i*self.plotFact + 1] +# SNRdB1 = SNRdB[i,:].reshape((1,-1)) +# axes.pcolorbuffer(x, y, SNRdB1, +# xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, +# xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, +# ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=0 + index = channelIndexList[i] + title = "Factor de reflectividad Z [dBZ]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap=colormap) + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=1 + index = channelIndexList[i] + title = "Velocidad vertical Doppler [m/s]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-10, zmax=10, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='seismic_r') + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=2 + index = channelIndexList[i] + title = "Intensidad de lluvia [mm/h]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=40, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='ocean_r') + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') self.draw() @@ -986,7 +1218,7 @@ class ParametersPlot(Figure): return dataOut @MPDecorator -class Parameters1Plot(Figure): +class Parameters1Plot_(Figure): __isConfig = None __nsubplots = None @@ -1067,9 +1299,8 @@ class Parameters1Plot(Figure): save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - #print inspect.getargspec(self.run).args - """ + """ Input: dataOut : id : @@ -1237,7 +1468,7 @@ class Parameters1Plot(Figure): update_figfile=False) return dataOut -class SpectralFittingPlot(Figure): +class SpectralFittingPlot_(Figure): __isConfig = None __nsubplots = None @@ -1415,7 +1646,7 @@ class SpectralFittingPlot(Figure): thisDatetime=thisDatetime) -class EWDriftsPlot(Figure): +class EWDriftsPlot_(Figure): __isConfig = None __nsubplots = None @@ -1621,7 +1852,7 @@ class EWDriftsPlot(Figure): -class PhasePlot(Figure): +class PhasePlot_(Figure): __isConfig = None __nsubplots = None @@ -1785,7 +2016,7 @@ class PhasePlot(Figure): -class NSMeteorDetection1Plot(Figure): +class NSMeteorDetection1Plot_(Figure): isConfig = None __nsubplots = None @@ -1969,7 +2200,7 @@ class NSMeteorDetection1Plot(Figure): thisDatetime=thisDatetime) -class NSMeteorDetection2Plot(Figure): +class NSMeteorDetection2Plot_(Figure): isConfig = None __nsubplots = None diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index d6ea2b6..1f6e1f3 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -14,7 +14,7 @@ from schainpy.model.proc.jroproc_base import MPDecorator from schainpy.utils import log @MPDecorator -class SpectraPlot(Figure): +class SpectraPlot_(Figure): isConfig = None __nsubplots = None @@ -42,6 +42,8 @@ class SpectraPlot(Figure): self.__xfilter_ena = False self.__yfilter_ena = False + + self.indice=1 def getSubplots(self): @@ -139,10 +141,9 @@ class SpectraPlot(Figure): x = dataOut.getVelRange(1) xlabel = "Velocity (m/s)" - ylabel = "Range (Km)" + ylabel = "Range (km)" y = dataOut.getHeiRange() - z = dataOut.data_spc/factor z = numpy.where(numpy.isfinite(z), z, numpy.NAN) zdB = 10*numpy.log10(z) @@ -155,6 +156,7 @@ class SpectraPlot(Figure): thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Spectra" + if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith) @@ -223,10 +225,11 @@ class SpectraPlot(Figure): ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime) + return dataOut @MPDecorator -class CrossSpectraPlot(Figure): +class CrossSpectraPlot_(Figure): isConfig = None __nsubplots = None @@ -252,6 +255,8 @@ class CrossSpectraPlot(Figure): self.EXP_CODE = None self.SUB_EXP_CODE = None self.PLOT_POS = None + + self.indice=0 def getSubplots(self): @@ -396,6 +401,7 @@ class CrossSpectraPlot(Figure): self.isConfig = True self.setWinTitle(title) + for i in range(self.nplots): pair = dataOut.pairsList[pairsIndexList[i]] @@ -420,7 +426,7 @@ class CrossSpectraPlot(Figure): xlabel=xlabel, ylabel=ylabel, title=title, ticksize=9, colormap=power_cmap, cblabel='') - coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:]) + coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] ) coherence = numpy.abs(coherenceComplex) # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi @@ -439,8 +445,6 @@ class CrossSpectraPlot(Figure): xlabel=xlabel, ylabel=ylabel, title=title, ticksize=9, colormap=phase_cmap, cblabel='') - - self.draw() self.save(figpath=figpath, @@ -453,7 +457,7 @@ class CrossSpectraPlot(Figure): return dataOut @MPDecorator -class RTIPlot(Figure): +class RTIPlot_(Figure): __isConfig = None __nsubplots = None @@ -470,7 +474,7 @@ class RTIPlot(Figure): self.__nsubplots = 1 self.WIDTH = 800 - self.HEIGHT = 180 + self.HEIGHT = 250 self.WIDTHPROF = 120 self.HEIGHTPROF = 0 self.counter_imagwr = 0 @@ -667,7 +671,7 @@ class RTIPlot(Figure): return dataOut @MPDecorator -class CoherenceMap(Figure): +class CoherenceMap_(Figure): isConfig = None __nsubplots = None @@ -878,7 +882,7 @@ class CoherenceMap(Figure): return dataOut @MPDecorator -class PowerProfilePlot(Figure): +class PowerProfilePlot_(Figure): isConfig = None __nsubplots = None @@ -1008,7 +1012,7 @@ class PowerProfilePlot(Figure): return dataOut @MPDecorator -class SpectraCutPlot(Figure): +class SpectraCutPlot_(Figure): isConfig = None __nsubplots = None @@ -1145,7 +1149,7 @@ class SpectraCutPlot(Figure): return dataOut @MPDecorator -class Noise(Figure): +class Noise_(Figure): isConfig = None __nsubplots = None @@ -1352,7 +1356,7 @@ class Noise(Figure): return dataOut @MPDecorator -class BeaconPhase(Figure): +class BeaconPhase_(Figure): __isConfig = None __nsubplots = None @@ -1497,9 +1501,6 @@ class BeaconPhase(Figure): avgcoherenceComplex = ccf/numpy.sqrt(powa*powb) phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi - #print "Phase %d%d" %(pair[0], pair[1]) - #print phase[dataOut.beacon_heiIndexList] - if dataOut.beacon_heiIndexList: phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList]) else: diff --git a/schainpy/model/graphics/jroplot_voltage.py b/schainpy/model/graphics/jroplot_voltage.py index 5544143..ae4a972 100644 --- a/schainpy/model/graphics/jroplot_voltage.py +++ b/schainpy/model/graphics/jroplot_voltage.py @@ -12,7 +12,7 @@ from .figure import Figure @MPDecorator -class Scope(Figure): +class Scope_(Figure): isConfig = None diff --git a/schainpy/model/graphics/jroplotter.py b/schainpy/model/graphics/jroplotter.py deleted file mode 100644 index d222d1c..0000000 --- a/schainpy/model/graphics/jroplotter.py +++ /dev/null @@ -1,240 +0,0 @@ -''' -Created on Jul 9, 2014 - -@author: roj-idl71 -''' -import os, sys -import datetime -import numpy -import traceback - -from time import sleep -from threading import Lock -# from threading import Thread - -import schainpy -import schainpy.admin - -from schainpy.model.proc.jroproc_base import Operation -from schainpy.model.serializer.data import obj2Dict, dict2Obj -from .jroplot_correlation import * -from .jroplot_heispectra import * -from .jroplot_parameters import * -from .jroplot_spectra import * -from .jroplot_voltage import * - - -class Plotter(Operation): - - isConfig = None - name = None - __queue = None - - def __init__(self, plotter_name, plotter_queue=None, **kwargs): - - Operation.__init__(self, **kwargs) - - self.isConfig = False - self.name = plotter_name - self.__queue = plotter_queue - - def getSubplots(self): - - nrow = self.nplots - ncol = 1 - return nrow, ncol - - def setup(self, **kwargs): - - print("Initializing ...") - - - def run(self, dataOut, id=None, **kwargs): - - """ - - Input: - dataOut : - id : - """ - - packDict = {} - - packDict['id'] = id - packDict['name'] = self.name - packDict['kwargs'] = kwargs - -# packDict['data'] = obj2Dict(dataOut) - packDict['data'] = dataOut - - self.__queue.put(packDict) - -# class PlotManager(Thread): -class PlotManager(): - - __err = False - __stop = False - __realtime = False - - controllerThreadObj = None - - plotterList = ['Scope', - 'SpectraPlot', 'RTIPlot', - 'SpectraCutPlot', - 'CrossSpectraPlot', 'CoherenceMap', - 'PowerProfilePlot', 'Noise', 'BeaconPhase', - 'CorrelationPlot', - 'SpectraHeisScope', 'RTIfromSpectraHeis'] - - def __init__(self, plotter_queue): - -# Thread.__init__(self) -# self.setDaemon(True) - - self.__queue = plotter_queue - self.__lock = Lock() - - self.plotInstanceDict = {} - - self.__err = False - self.__stop = False - self.__realtime = False - - def __handleError(self, name="", send_email=False): - - err = traceback.format_exception(sys.exc_info()[0], - sys.exc_info()[1], - sys.exc_info()[2]) - - print("***** Error occurred in PlotManager *****") - print("***** [%s]: %s" %(name, err[-1])) - - message = "\nError ocurred in %s:\n" %name - message += "".join(err) - - sys.stderr.write(message) - - if not send_email: - return - - import socket - - subject = "SChain v%s: Error running %s\n" %(schainpy.__version__, name) - - subtitle = "%s:\n" %(name) - subtitle += "Hostname: %s\n" %socket.gethostbyname(socket.gethostname()) - subtitle += "Working directory: %s\n" %os.path.abspath("./") - # subtitle += "Configuration file: %s\n" %self.filename - subtitle += "Time: %s\n" %str(datetime.datetime.now()) - - adminObj = schainpy.admin.SchainNotify() - adminObj.sendAlert(message=message, - subject=subject, - subtitle=subtitle) - - def run(self): - - if self.__queue.empty(): - return - - if self.__err: - serial_data = self.__queue.get() - self.__queue.task_done() - return - - self.__lock.acquire() - -# if self.__queue.full(): -# for i in range(int(self.__queue.qsize()/2)): -# serial_data = self.__queue.get() -# self.__queue.task_done() - - n = int(self.__queue.qsize()/3 + 1) - - for i in range(n): - - if self.__queue.empty(): - break - - serial_data = self.__queue.get() - self.__queue.task_done() - - plot_id = serial_data['id'] - plot_name = serial_data['name'] - kwargs = serial_data['kwargs'] -# dataDict = serial_data['data'] -# -# dataPlot = dict2Obj(dataDict) - - dataPlot = serial_data['data'] - - if plot_id not in list(self.plotInstanceDict.keys()): - className = eval(plot_name) - self.plotInstanceDict[plot_id] = className(**kwargs) - - plotter = self.plotInstanceDict[plot_id] - try: - plotter.run(dataPlot, plot_id, **kwargs) - except: - self.__err = True - self.__handleError(plot_name, send_email=True) - break - - self.__lock.release() - - def isEmpty(self): - - return self.__queue.empty() - - def stop(self): - - self.__lock.acquire() - - self.__stop = True - - self.__lock.release() - - def close(self): - - self.__lock.acquire() - - for plot_id in list(self.plotInstanceDict.keys()): - plotter = self.plotInstanceDict[plot_id] - plotter.close() - - self.__lock.release() - - def setController(self, controllerThreadObj): - - self.controllerThreadObj = controllerThreadObj - - def start(self): - - if not self.controllerThreadObj.isRunning(): - raise RuntimeError("controllerThreadObj has not been initialized. Use controllerThreadObj.start() before call this method") - - self.join() - - def join(self): - - #Execute plotter while controller is running - while self.controllerThreadObj.isRunning(): - self.run() - - self.controllerThreadObj.stop() - - #Wait until plotter queue is empty - while not self.isEmpty(): - self.run() - - self.close() - - def isErrorDetected(self): - - self.__lock.acquire() - - err = self.__err - - self.__lock.release() - - return err \ No newline at end of file diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index 536daa2..9b4b1b6 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -434,7 +434,6 @@ def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='' def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''): ax = iplot.axes - printLabels(ax, xlabel, ylabel, title) for i in range(len(ax.lines)): diff --git a/schainpy/model/io/bltrIO_param.py b/schainpy/model/io/bltrIO_param.py index 7048346..70a996c 100644 --- a/schainpy/model/io/bltrIO_param.py +++ b/schainpy/model/io/bltrIO_param.py @@ -13,7 +13,7 @@ import datetime import numpy -from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator +from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator from schainpy.model.data.jrodata import Parameters from schainpy.model.io.jroIO_base import JRODataReader, isNumber from schainpy.utils import log @@ -111,7 +111,6 @@ class BLTRParamReader(JRODataReader, ProcessingUnit): timezone=0, status_value=0, **kwargs): - self.path = path self.startDate = startDate self.endDate = endDate @@ -185,7 +184,6 @@ class BLTRParamReader(JRODataReader, ProcessingUnit): file_id = self.fileIndex if file_id == len(self.fileList): - log.success('No more files in the folder', 'BLTRParamReader') self.flagNoMoreFiles = 1 return 0 @@ -240,7 +238,7 @@ class BLTRParamReader(JRODataReader, ProcessingUnit): pointer = self.fp.tell() header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1) - self.nchannels = header_rec['nchan'][0] / 2 + self.nchannels = int(header_rec['nchan'][0] / 2) self.kchan = header_rec['nrxs'][0] self.nmodes = header_rec['nmodes'][0] self.nranges = header_rec['nranges'][0] @@ -358,8 +356,7 @@ class BLTRParamReader(JRODataReader, ProcessingUnit): ''' if self.flagNoMoreFiles: self.dataOut.flagNoData = True - log.success('No file left to process', 'BLTRParamReader') - return 0 + self.dataOut.error = (1, 'No More files to read') if not self.readNextBlock(): self.dataOut.flagNoData = True diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py index 2a6c9f4..941ea92 100644 --- a/schainpy/model/io/jroIO_base.py +++ b/schainpy/model/io/jroIO_base.py @@ -1815,7 +1815,7 @@ class JRODataWriter(JRODataIO): return 1 - def run(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs): + def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs): if not(self.isConfig): diff --git a/schainpy/model/io/jroIO_bltr.py b/schainpy/model/io/jroIO_bltr.py index dd29481..593d3cb 100644 --- a/schainpy/model/io/jroIO_bltr.py +++ b/schainpy/model/io/jroIO_bltr.py @@ -63,9 +63,6 @@ class Header(object): if attr: message += "%s = %s" % ("size", attr) + "\n" - # print message - - FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes ('FileMgcNumber', ' endFp: sys.stderr.write( "Warning %s: Size value read from System Header is lower than it has to be\n" % fp) @@ -537,9 +440,6 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa FileList.append(IndexFile) nFiles += 1 - # print 'Files2Read' - # print 'Existen '+str(nFiles)+' archivos .fdt' - self.filenameList = FileList # List of files from least to largest by names def run(self, **kwargs): @@ -553,7 +453,6 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.isConfig = True self.getData() - # print 'running' def setup(self, path=None, startDate=None, @@ -590,22 +489,19 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa if self.flagNoMoreFiles: self.dataOut.flagNoData = True - print('NoData se vuelve true') return 0 self.fp = self.path self.Files2Read(self.fp) self.readFile(self.fp) self.dataOut.data_spc = self.data_spc - self.dataOut.data_cspc = self.data_cspc - self.dataOut.data_output = self.data_output - - print('self.dataOut.data_output', shape(self.dataOut.data_output)) - - # self.removeDC() - return self.dataOut.data_spc - - def readFile(self, fp): + self.dataOut.data_cspc =self.data_cspc + self.dataOut.data_output=self.data_output + + return self.dataOut.data_spc + + + def readFile(self,fp): ''' You must indicate if you are reading in Online or Offline mode and load the The parameters for this file reading mode. @@ -615,23 +511,18 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa 1. Get the BLTR FileHeader. 2. Start reading the first block. ''' - - # The address of the folder is generated the name of the .fdt file that will be read - print("File: ", self.fileSelector + 1) - + if self.fileSelector < len(self.filenameList): self.fpFile = str(fp) + '/' + \ str(self.filenameList[self.fileSelector]) - # print self.fpFile fheader = FileHeaderBLTR() fheader.FHread(self.fpFile) # Bltr FileHeader Reading self.nFDTdataRecors = fheader.nFDTdataRecors self.readBlock() # Block reading else: - print('readFile FlagNoData becomes true') - self.flagNoMoreFiles = True + self.flagNoMoreFiles=True self.dataOut.flagNoData = True return 0 @@ -658,12 +549,11 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa 2. Fill the buffer with the current block number. ''' - - if self.BlockCounter < self.nFDTdataRecors - 2: - print(self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') - if self.ReadMode == 1: - rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1) - elif self.ReadMode == 0: + + if self.BlockCounter < self.nFDTdataRecors-1: + if self.ReadMode==1: + rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1) + elif self.ReadMode==0: rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter) rheader.RHread(self.fpFile) # Bltr FileHeader Reading @@ -683,31 +573,26 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.nRdPairs = len(self.dataOut.pairsList) self.dataOut.nRdPairs = self.nRdPairs - - self.__firstHeigth = rheader.StartRangeSamp - self.__deltaHeigth = rheader.SampResolution - self.dataOut.heightList = self.__firstHeigth + \ - numpy.array(list(range(self.nHeights))) * self.__deltaHeigth - self.dataOut.channelList = list(range(self.nChannels)) - self.dataOut.nProfiles = rheader.nProfiles - self.dataOut.nIncohInt = rheader.nIncohInt - self.dataOut.nCohInt = rheader.nCohInt - self.dataOut.ippSeconds = 1 / float(rheader.PRFhz) - self.dataOut.PRF = rheader.PRFhz - self.dataOut.nFFTPoints = rheader.nProfiles - self.dataOut.utctime = rheader.nUtime - self.dataOut.timeZone = 0 - self.dataOut.normFactor = self.dataOut.nProfiles * \ - self.dataOut.nIncohInt * self.dataOut.nCohInt - self.dataOut.outputInterval = self.dataOut.ippSeconds * \ - self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles - - self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN - print('self.data_output', shape(self.data_output)) - self.dataOut.velocityX = [] - self.dataOut.velocityY = [] - self.dataOut.velocityV = [] - + self.__firstHeigth=rheader.StartRangeSamp + self.__deltaHeigth=rheader.SampResolution + self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth + self.dataOut.channelList = range(self.nChannels) + self.dataOut.nProfiles=rheader.nProfiles + self.dataOut.nIncohInt=rheader.nIncohInt + self.dataOut.nCohInt=rheader.nCohInt + self.dataOut.ippSeconds= 1/float(rheader.PRFhz) + self.dataOut.PRF=rheader.PRFhz + self.dataOut.nFFTPoints=rheader.nProfiles + self.dataOut.utctime=rheader.nUtime + self.dataOut.timeZone=0 + self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt + self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles + + self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN + self.dataOut.velocityX=[] + self.dataOut.velocityY=[] + self.dataOut.velocityV=[] + '''Block Reading, the Block Data is received and Reshape is used to give it shape. ''' @@ -734,18 +619,17 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa y = rho * numpy.sin(phi) return(x, y) - if self.DualModeIndex == self.ReadMode: - - self.data_fft = numpy.fromfile( - startDATA, [('complex', ' 0.0001) : -# -# try: -# popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) -# -# if numpy.amax(popt)>numpy.amax(yMean)*0.3: -# FitGauss=gaus(xSamples,*popt) -# -# else: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# print 'Verificador: Dentro', Height -# except RuntimeError: -# -# try: -# for j in range(len(ySamples[1])): -# yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]])) -# popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma]) -# FitGauss=gaus(xSamples,*popt) -# print 'Verificador: Exepcion1', Height -# except RuntimeError: -# -# try: -# popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma]) -# FitGauss=gaus(xSamples,*popt) -# print 'Verificador: Exepcion2', Height -# except RuntimeError: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# print 'Verificador: Exepcion3', Height -# else: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# #print 'Verificador: Fuera', Height -# -# -# -# Maximun=numpy.amax(yMean) -# eMinus1=Maximun*numpy.exp(-1) -# -# HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) -# HalfWidth= xFrec[HWpos] -# GCpos=Find(FitGauss, numpy.amax(FitGauss)) -# Vpos=Find(FactNorm, numpy.amax(FactNorm)) -# #Vpos=numpy.sum(FactNorm)/len(FactNorm) -# #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) ))) -# #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos -# '''****** Getting Fij ******''' -# -# GaussCenter=xFrec[GCpos] -# if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): -# Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 -# else: -# Fij=abs(GaussCenter-HalfWidth)+0.0000001 -# -# '''****** Getting Frecuency range of significant data ******''' -# -# Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) -# -# if Rangpos5 and len(FrecRange) 0.: -# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag -# #print 'Vmag',Vmag -# else: -# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN) -# -# if abs(Vx)<100 and abs(Vx) > 0.: -# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang -# #print 'Vang',Vang -# else: -# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN) -# -# if abs(GaussCenter)<2: -# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos]) -# -# else: -# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN) -# -# -# # print '********************************************' -# # print 'HalfWidth ', HalfWidth -# # print 'Maximun ', Maximun -# # print 'eMinus1 ', eMinus1 -# # print 'Rangpos ', Rangpos -# # print 'GaussCenter ',GaussCenter -# # print 'E01 ',E01 -# # print 'N01 ',N01 -# # print 'E02 ',E02 -# # print 'N02 ',N02 -# # print 'E12 ',E12 -# # print 'N12 ',N12 -# #print 'self.dataOut.velocityX ', self.dataOut.velocityX -# # print 'Fij ', Fij -# # print 'cC ', cC -# # print 'cF ', cF -# # print 'cG ', cG -# # print 'cA ', cA -# # print 'cB ', cB -# # print 'cH ', cH -# # print 'Vx ', Vx -# # print 'Vy ', Vy -# # print 'Vmag ', Vmag -# # print 'Vang ', Vang*180/numpy.pi -# # print 'PhaseSlope ',PhaseSlope[0] -# # print 'PhaseSlope ',PhaseSlope[1] -# # print 'PhaseSlope ',PhaseSlope[2] -# # print '********************************************' -# #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY) -# -# #print 'self.dataOut.velocityX', len(self.dataOut.velocityX) -# #print 'self.dataOut.velocityY', len(self.dataOut.velocityY) -# #print 'self.dataOut.velocityV', self.dataOut.velocityV -# -# self.data_output[0]=numpy.array(self.dataOut.velocityX) -# self.data_output[1]=numpy.array(self.dataOut.velocityY) -# self.data_output[2]=numpy.array(self.dataOut.velocityV) -# -# prin= self.data_output[0][~numpy.isnan(self.data_output[0])] -# print ' ' -# print 'VmagAverage',numpy.mean(prin) -# print ' ' -# # plt.figure(5) -# # plt.subplot(211) -# # plt.plot(self.dataOut.velocityX,'yo:') -# # plt.subplot(212) -# # plt.plot(self.dataOut.velocityY,'yo:') -# -# # plt.figure(1) -# # # plt.subplot(121) -# # # plt.plot(xFrec,ySamples[0],'k',label='Ch0') -# # # plt.plot(xFrec,ySamples[1],'g',label='Ch1') -# # # plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') -# # # plt.legend() -# # plt.title('DATOS A ALTURA DE 2850 METROS') -# # -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # # plt.subplot(122) -# # # plt.title('Fit for Time Constant') -# # #plt.plot(xFrec,zline) -# # #plt.plot(xFrec,SmoothSPC,'g') -# # plt.plot(xFrec,FactNorm) -# # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # -# # plt.figure(10) -# # # plt.subplot(121) -# # plt.plot(xFrec,ySamples[0],'b',label='Ch0') -# # plt.plot(xFrec,ySamples[1],'y',label='Ch1') -# # plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') -# # plt.legend() -# # plt.title('SELFSPECTRA EN CANALES') -# # -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # # plt.subplot(122) -# # # plt.title('Fit for Time Constant') -# # #plt.plot(xFrec,zline) -# # #plt.plot(xFrec,SmoothSPC,'g') -# # # plt.plot(xFrec,FactNorm) -# # # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # -# # plt.figure(9) -# # -# # -# # plt.title('DATOS SUAVIZADOS') -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # plt.plot(xFrec,SmoothSPC,'g') -# # -# # #plt.plot(xFrec,FactNorm) -# # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # # -# # plt.figure(2) -# # # #plt.subplot(121) -# # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra') -# # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano') -# # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo') -# # # #plt.plot(xFrec,phase) -# # # plt.xlabel('Suavizado, promediado KHz') -# # plt.title('SELFSPECTRA PROMEDIADO') -# # # #plt.subplot(122) -# # # #plt.plot(xSamples,zline) -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # plt.legend() -# # # -# # # plt.figure(3) -# # # plt.subplot(311) -# # # #plt.plot(xFrec,phase[0]) -# # # plt.plot(xFrec,phase[0],'g') -# # # plt.subplot(312) -# # # plt.plot(xFrec,phase[1],'g') -# # # plt.subplot(313) -# # # plt.plot(xFrec,phase[2],'g') -# # # #plt.plot(xFrec,phase[2]) -# # # -# # # plt.figure(4) -# # # -# # # plt.plot(xSamples,coherence[0],'b') -# # # plt.plot(xSamples,coherence[1],'r') -# # # plt.plot(xSamples,coherence[2],'g') -# # plt.show() -# # # -# # # plt.clf() -# # # plt.cla() -# # # plt.close() -# -# print ' ' - - self.BlockCounter += 2 - + self.BlockCounter+=2 + else: - self.fileSelector += 1 - self.BlockCounter = 0 - print("Next File") \ No newline at end of file + self.fileSelector+=1 + self.BlockCounter=0 diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index 2845d8a..322769a 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -179,9 +179,6 @@ class ParamReader(JRODataReader,ProcessingUnit): print("[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)) print() -# for i in range(len(filenameList)): -# print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime()) - self.filenameList = filenameList self.datetimeList = datetimeList @@ -504,20 +501,11 @@ class ParamReader(JRODataReader,ProcessingUnit): def getData(self): -# if self.flagNoMoreFiles: -# self.dataOut.flagNoData = True -# print 'Process finished' -# return 0 -# if self.blockIndex==self.blocksPerFile: if not( self.__setNextFileOffline() ): self.dataOut.flagNoData = True return 0 -# if self.datablock == None: # setear esta condicion cuando no hayan datos por leers -# self.dataOut.flagNoData = True -# return 0 -# self.__readData() self.__setDataOut() self.dataOut.flagNoData = False @@ -637,7 +625,10 @@ class ParamWriter(Operation): dsDict['variable'] = self.dataList[i] #--------------------- Conditionals ------------------------ #There is no data + + if dataAux is None: + return 0 #Not array, just a number @@ -821,7 +812,7 @@ class ParamWriter(Operation): return False def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -1095,7 +1086,6 @@ class ParamWriter(Operation): return self.isConfig = True -# self.putMetadata() self.setNextFile() self.putData() diff --git a/schainpy/model/io/jroIO_spectra.py b/schainpy/model/io/jroIO_spectra.py index ea4b4af..6c523d7 100644 --- a/schainpy/model/io/jroIO_spectra.py +++ b/schainpy/model/io/jroIO_spectra.py @@ -413,9 +413,7 @@ class SpectraWriter(JRODataWriter, Operation): data_dc = None -# dataOut = None - - def __init__(self):#, **kwargs): + def __init__(self): """ Inicializador de la clase SpectraWriter para la escritura de datos de espectros. @@ -429,9 +427,7 @@ class SpectraWriter(JRODataWriter, Operation): Return: None """ - Operation.__init__(self)#, **kwargs) - - #self.isConfig = False + Operation.__init__(self) self.nTotalBlocks = 0 @@ -496,7 +492,7 @@ class SpectraWriter(JRODataWriter, Operation): def writeBlock(self): - """ + """processingHeaderObj Escribe el buffer en el file designado Affected: @@ -519,8 +515,10 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) if self.data_cspc is not None: - data = numpy.zeros( self.shape_cspc_Buffer, self.dtype ) + cspc = numpy.transpose( self.data_cspc, (0,2,1) ) + #data = numpy.zeros( numpy.shape(cspc), self.dtype ) + #print 'data.shape', self.shape_cspc_Buffer if not self.processingHeaderObj.shif_fft: cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones data['real'] = cspc.real @@ -529,8 +527,9 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) if self.data_dc is not None: - data = numpy.zeros( self.shape_dc_Buffer, self.dtype ) + dc = self.data_dc + data = numpy.zeros( numpy.shape(dc), self.dtype ) data['real'] = dc.real data['imag'] = dc.imag data = data.reshape((-1)) diff --git a/schainpy/model/proc/bltrproc_parameters.py b/schainpy/model/proc/bltrproc_parameters.py index a043562..30fcc63 100644 --- a/schainpy/model/proc/bltrproc_parameters.py +++ b/schainpy/model/proc/bltrproc_parameters.py @@ -12,13 +12,11 @@ from time import gmtime from numpy import transpose -from .jroproc_base import ProcessingUnit, MPDecorator, Operation +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator from schainpy.model.data.jrodata import Parameters @MPDecorator -class BLTRParametersProc(ProcessingUnit): - - METHODS = {} +class BLTRParametersProc(ProcessingUnit): ''' Processing unit for BLTR parameters data (winds) @@ -46,9 +44,7 @@ class BLTRParametersProc(ProcessingUnit): Inputs: None ''' ProcessingUnit.__init__(self) - self.setupReq = False self.dataOut = Parameters() - self.isConfig = False def setup(self, mode): ''' diff --git a/schainpy/model/proc/jroproc_base.py b/schainpy/model/proc/jroproc_base.py index 6e2a4e9..5dda991 100644 --- a/schainpy/model/proc/jroproc_base.py +++ b/schainpy/model/proc/jroproc_base.py @@ -11,13 +11,15 @@ Based on: $Author: murco $ $Id: jroproc_base.py 1 2012-11-12 18:56:07Z murco $ ''' -from platform import python_version + import inspect import zmq import time import pickle import os from multiprocessing import Process +from zmq.utils.monitor import recv_monitor_message + from schainpy.utils import log @@ -35,15 +37,6 @@ class ProcessingUnit(object): """ - - METHODS = {} - dataIn = None - dataInList = [] - id = None - inputId = None - dataOut = None - dictProcs = None - isConfig = False def __init__(self): @@ -51,15 +44,15 @@ class ProcessingUnit(object): self.dataOut = None self.isConfig = False self.operations = [] + self.plots = [] def getAllowedArgs(self): if hasattr(self, '__attrs__'): return self.__attrs__ else: return inspect.getargspec(self.run).args - - def addOperation(self, conf, operation): + def addOperation(self, conf, operation): """ This method is used in the controller, and update the dictionary containing the operations to execute. The dict posses the id of the operation process (IPC purposes) @@ -76,7 +69,11 @@ class ProcessingUnit(object): objId : identificador del objeto, necesario para comunicar con master(procUnit) """ - self.operations.append((operation, conf.type, conf.id, conf.getKwargs())) + self.operations.append( + (operation, conf.type, conf.id, conf.getKwargs())) + + if 'plot' in self.name.lower(): + self.plots.append(operation.CODE) def getOperationObj(self, objId): @@ -86,7 +83,6 @@ class ProcessingUnit(object): return self.operations[objId] def operation(self, **kwargs): - """ Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los atributos del objeto dataOut @@ -96,7 +92,7 @@ class ProcessingUnit(object): **kwargs : Diccionario de argumentos de la funcion a ejecutar """ - raise NotImplementedError + raise NotImplementedError def setup(self): @@ -107,9 +103,10 @@ class ProcessingUnit(object): raise NotImplementedError def close(self): - #Close every thread, queue or any other object here is it is neccesary. + return - + + class Operation(object): """ @@ -126,22 +123,15 @@ class Operation(object): Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer) """ - id = None - __buffer = None - dest = None - isConfig = False - readyFlag = None def __init__(self): - self.buffer = None - self.dest = None + self.id = None self.isConfig = False - self.readyFlag = False if not hasattr(self, 'name'): self.name = self.__class__.__name__ - + def getAllowedArgs(self): if hasattr(self, '__attrs__'): return self.__attrs__ @@ -154,9 +144,7 @@ class Operation(object): raise NotImplementedError - def run(self, dataIn, **kwargs): - """ Realiza las operaciones necesarias sobre la dataIn.data y actualiza los atributos del objeto dataIn. @@ -180,18 +168,17 @@ class Operation(object): def close(self): - pass + return def MPDecorator(BaseClass): - """ Multiprocessing class decorator This function add multiprocessing features to a BaseClass. Also, it handle the communication beetween processes (readers, procUnits and operations). """ - + class MPClass(BaseClass, Process): def __init__(self, *args, **kwargs): @@ -203,42 +190,38 @@ def MPDecorator(BaseClass): self.sender = None self.receiver = None self.name = BaseClass.__name__ + self.start_time = time.time() if len(self.args) is 3: self.typeProc = "ProcUnit" - self.id = args[0] - self.inputId = args[1] - self.project_id = args[2] - else: + self.id = args[0] + self.inputId = args[1] + self.project_id = args[2] + elif len(self.args) is 2: self.id = args[0] self.inputId = args[0] self.project_id = args[1] self.typeProc = "Operation" - def getAllowedArgs(self): - - if hasattr(self, '__attrs__'): - return self.__attrs__ - else: - return inspect.getargspec(BaseClass.run).args - def subscribe(self): ''' This function create a socket to receive objects from the topic `inputId`. ''' - + c = zmq.Context() self.receiver = c.socket(zmq.SUB) - self.receiver.connect('ipc:///tmp/schain/{}_pub'.format(self.project_id)) + self.receiver.connect( + 'ipc:///tmp/schain/{}_pub'.format(self.project_id)) self.receiver.setsockopt(zmq.SUBSCRIBE, self.inputId.encode()) def listen(self): ''' This function waits for objects and deserialize using pickle ''' - - data = pickle.loads(self.receiver.recv_multipart()[1]) + + data = pickle.loads(self.receiver.recv_multipart()[1]) + return data def set_publisher(self): @@ -248,14 +231,14 @@ def MPDecorator(BaseClass): time.sleep(1) c = zmq.Context() - self.sender = c.socket(zmq.PUB) - self.sender.connect('ipc:///tmp/schain/{}_sub'.format(self.project_id)) + self.sender = c.socket(zmq.PUB) + self.sender.connect( + 'ipc:///tmp/schain/{}_sub'.format(self.project_id)) - def publish(self, data, id): + def publish(self, data, id): ''' This function publish an object, to a specific topic. ''' - self.sender.send_multipart([str(id).encode(), pickle.dumps(data)]) def runReader(self): @@ -263,28 +246,32 @@ def MPDecorator(BaseClass): Run fuction for read units ''' while True: - - BaseClass.run(self, **self.kwargs) - if self.dataOut.error[0] == -1: - log.error(self.dataOut.error[1]) - self.publish('end', self.id) - #self.sender.send_multipart([str(self.project_id).encode(), 'end'.encode()]) - break + BaseClass.run(self, **self.kwargs) - for op, optype, id, kwargs in self.operations: - if optype=='self': + for op, optype, opId, kwargs in self.operations: + if optype == 'self': op(**kwargs) - elif optype=='other': + elif optype == 'other': self.dataOut = op.run(self.dataOut, **self.kwargs) - elif optype=='external': + elif optype == 'external': self.publish(self.dataOut, opId) - if self.dataOut.flagNoData: + if self.dataOut.flagNoData and self.dataOut.error is None: continue - - self.publish(self.dataOut, self.id) - + + self.publish(self.dataOut, self.id) + + if self.dataOut.error: + if self.dataOut.error[0] == -1: + log.error(self.dataOut.error[1], self.name) + if self.dataOut.error[0] == 1: + log.success(self.dataOut.error[1], self.name) + # self.sender.send_multipart([str(self.project_id).encode(), 'end'.encode()]) + break + + time.sleep(1) + def runProc(self): ''' Run function for proccessing units @@ -293,49 +280,45 @@ def MPDecorator(BaseClass): while True: self.dataIn = self.listen() - if self.dataIn == 'end': - self.publish('end', self.id) - for op, optype, opId, kwargs in self.operations: - if optype == 'external': - self.publish('end', opId) - break - - if self.dataIn.flagNoData: + if self.dataIn.flagNoData and self.dataIn.error is None: continue BaseClass.run(self, **self.kwargs) for op, optype, opId, kwargs in self.operations: - if optype=='self': + if optype == 'self': op(**kwargs) - elif optype=='other': + elif optype == 'other': self.dataOut = op.run(self.dataOut, **kwargs) - elif optype=='external': + elif optype == 'external': self.publish(self.dataOut, opId) - if self.dataOut.flagNoData: - continue - self.publish(self.dataOut, self.id) + if self.dataIn.error: + break + + time.sleep(1) def runOp(self): ''' - Run function for operations + Run function for external operations (this operations just receive data + ex: plots, writers, publishers) ''' while True: dataOut = self.listen() - if dataOut == 'end': - break - BaseClass.run(self, dataOut, **self.kwargs) - + + if dataOut.error: + break + time.sleep(1) + def run(self): - + if self.typeProc is "ProcUnit": - + if self.inputId is not None: self.subscribe() self.set_publisher() @@ -346,22 +329,48 @@ def MPDecorator(BaseClass): self.runReader() elif self.typeProc is "Operation": - + self.subscribe() self.runOp() else: raise ValueError("Unknown type") - print("%s done" % BaseClass.__name__) self.close() + def event_monitor(self, monitor): + + events = {} + + for name in dir(zmq): + if name.startswith('EVENT_'): + value = getattr(zmq, name) + events[value] = name + + while monitor.poll(): + evt = recv_monitor_message(monitor) + if evt['event'] == 32: + self.connections += 1 + if evt['event'] == 512: + pass + + evt.update({'description': events[evt['event']]}) + + if evt['event'] == zmq.EVENT_MONITOR_STOPPED: + break + monitor.close() + print('event monitor thread done!') + def close(self): - + + BaseClass.close(self) + if self.sender: self.sender.close() - + if self.receiver: self.receiver.close() - return MPClass \ No newline at end of file + log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name) + + return MPClass diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py old mode 100644 new mode 100755 index 737ed9e..5c4e2f1 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -10,11 +10,7 @@ import importlib import itertools from multiprocessing import Pool, TimeoutError from multiprocessing.pool import ThreadPool -import types -from functools import partial import time -#from sklearn.cluster import KMeans - from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters from .jroproc_base import ProcessingUnit, Operation, MPDecorator @@ -128,6 +124,7 @@ class ParametersProc(ProcessingUnit): self.dataOut.abscissaList = self.dataIn.getVelRange(1) self.dataOut.spc_noise = self.dataIn.getNoise() self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1)) + # self.dataOut.normFactor = self.dataIn.normFactor self.dataOut.pairsList = self.dataIn.pairsList self.dataOut.groupList = self.dataIn.pairsList self.dataOut.flagNoData = False @@ -136,9 +133,9 @@ class ParametersProc(ProcessingUnit): self.dataOut.ChanDist = self.dataIn.ChanDist else: self.dataOut.ChanDist = None - if hasattr(self.dataIn, 'VelRange'): #Velocities range - self.dataOut.VelRange = self.dataIn.VelRange - else: self.dataOut.VelRange = None + #if hasattr(self.dataIn, 'VelRange'): #Velocities range + # self.dataOut.VelRange = self.dataIn.VelRange + #else: self.dataOut.VelRange = None if hasattr(self.dataIn, 'RadarConst'): #Radar Constant self.dataOut.RadarConst = self.dataIn.RadarConst @@ -184,9 +181,112 @@ class ParametersProc(ProcessingUnit): def target(tups): obj, args = tups - #print 'TARGETTT', obj, args + return obj.FitGau(args) + +class SpectralFilters(Operation): + + '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR + + LimitR : It is the limit in m/s of Rainfall + LimitW : It is the limit in m/s for Winds + + Input: + + self.dataOut.data_pre : SPC and CSPC + self.dataOut.spc_range : To select wind and rainfall velocities + + Affected: + + self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind + self.dataOut.spcparam_range : Used in SpcParamPlot + self.dataOut.SPCparam : Used in PrecipitationProc + + + ''' + + def __init__(self): + Operation.__init__(self) + self.i=0 + + def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5): + + + #Limite de vientos + LimitR = PositiveLimit + LimitN = NegativeLimit + + self.spc = dataOut.data_pre[0].copy() + self.cspc = dataOut.data_pre[1].copy() + + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] + + VelRange = dataOut.spc_range[2] + TimeRange = dataOut.spc_range[1] + FrecRange = dataOut.spc_range[0] + + Vmax= 2*numpy.max(dataOut.spc_range[2]) + Tmax= 2*numpy.max(dataOut.spc_range[1]) + Fmax= 2*numpy.max(dataOut.spc_range[0]) + + Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()] + Breaker1R=numpy.where(VelRange == Breaker1R) + + Delta = self.Num_Bin/2 - Breaker1R[0] + + + '''Reacomodando SPCrange''' + + VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0) + + VelRange[-(self.Num_Bin/2):]+= Vmax + + FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0) + + FrecRange[-(self.Num_Bin/2):]+= Fmax + + TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0) + + TimeRange[-(self.Num_Bin/2):]+= Tmax + + ''' ------------------ ''' + + Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()] + Breaker2R=numpy.where(VelRange == Breaker2R) + + + SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1) + + SPCcut = SPCroll.copy() + for i in range(self.Num_Chn): + + SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i] + SPCcut[i,-int(Delta):,:] = dataOut.noise[i] + + SPCcut[i]=SPCcut[i]- dataOut.noise[i] + SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20 + + SPCroll[i]=SPCroll[i]-dataOut.noise[i] + SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20 + + SPC_ch1 = SPCroll + + SPC_ch2 = SPCcut + + SPCparam = (SPC_ch1, SPC_ch2, self.spc) + dataOut.SPCparam = numpy.asarray(SPCparam) + + + dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1]) + + dataOut.spcparam_range[2]=VelRange + dataOut.spcparam_range[1]=TimeRange + dataOut.spcparam_range[0]=FrecRange + + class GaussianFit(Operation): ''' @@ -198,15 +298,15 @@ class GaussianFit(Operation): self.dataOut.data_pre : SelfSpectra Output: - self.dataOut.GauSPC : SPC_ch1, SPC_ch2 + self.dataOut.SPCparam : SPC_ch1, SPC_ch2 ''' - def __init__(self, **kwargs): - Operation.__init__(self, **kwargs) + def __init__(self): + Operation.__init__(self) self.i=0 - def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points + 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 """This routine will find a couple of generalized Gaussians to a power spectrum input: spc output: @@ -214,31 +314,12 @@ class GaussianFit(Operation): """ self.spc = dataOut.data_pre[0].copy() - - - print('SelfSpectra Shape', numpy.asarray(self.spc).shape) - - - #plt.figure(50) - #plt.subplot(121) - #plt.plot(self.spc,'k',label='spc(66)') - #plt.plot(xFrec,ySamples[1],'g',label='Ch1') - #plt.plot(xFrec,ySamples[2],'r',label='Ch2') - #plt.plot(xFrec,FitGauss,'yo:',label='fit') - #plt.legend() - #plt.title('DATOS A ALTURA DE 7500 METROS') - #plt.show() - self.Num_Hei = self.spc.shape[2] - #self.Num_Bin = len(self.spc) self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] - Vrange = dataOut.abscissaList - #print 'self.spc2', numpy.asarray(self.spc).shape - - GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei]) + GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei]) SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch1[:] = numpy.NaN @@ -250,272 +331,12 @@ class GaussianFit(Operation): noise_ = dataOut.spc_noise[0].copy() - pool = Pool(processes=self.Num_Chn) args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)] objs = [self for __ in range(self.Num_Chn)] attrs = list(zip(objs, args)) gauSPC = pool.map(target, attrs) - dataOut.GauSPC = numpy.asarray(gauSPC) -# ret = [] -# for n in range(self.Num_Chn): -# self.FitGau(args[n]) -# dataOut.GauSPC = ret - - - -# for ch in range(self.Num_Chn): -# -# for ht in range(self.Num_Hei): -# #print (numpy.asarray(self.spc).shape) -# spc = numpy.asarray(self.spc)[ch,:,ht] -# -# ############################################# -# # normalizing spc and noise -# # This part differs from gg1 -# spc_norm_max = max(spc) -# spc = spc / spc_norm_max -# pnoise = pnoise / spc_norm_max -# ############################################# -# -# if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's -# fatspectra=1.0 -# else: -# fatspectra=0.5 -# -# wnoise = noise_ / spc_norm_max -# #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise -# #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used -# #if wnoise>1.1*pnoise: # to be tested later -# # wnoise=pnoise -# noisebl=wnoise*0.9; noisebh=wnoise*1.1 -# spc=spc-wnoise -# -# minx=numpy.argmin(spc) -# spcs=numpy.roll(spc,-minx) -# cum=numpy.cumsum(spcs) -# tot_noise=wnoise * self.Num_Bin #64; -# #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' -# #snr=tot_signal/tot_noise -# #snr=cum[-1]/tot_noise -# -# #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise -# -# snr = sum(spcs)/tot_noise -# snrdB=10.*numpy.log10(snr) -# -# #if snrdB < -9 : -# # snrdB = numpy.NaN -# # continue -# -# #print 'snr',snrdB # , sum(spcs) , tot_noise -# -# -# #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: -# # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None -# -# cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region -# cumlo=cummax*epsi; -# cumhi=cummax*(1-epsi) -# powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-9: # when SNR is strong pick the peak with least shift (LOS velocity) error -# if oneG: -# choice=0 -# else: -# w1=lsq2[0][1]; w2=lsq2[0][5] -# a1=lsq2[0][2]; a2=lsq2[0][6] -# p1=lsq2[0][3]; p2=lsq2[0][7] -# s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; -# gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling -# -# if gp1>gp2: -# if a1>0.7*a2: -# choice=1 -# else: -# choice=2 -# elif gp2>gp1: -# if a2>0.7*a1: -# choice=2 -# else: -# choice=1 -# else: -# choice=numpy.argmax([a1,a2])+1 -# #else: -# #choice=argmin([std2a,std2b])+1 -# -# else: # with low SNR go to the most energetic peak -# choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) -# -# #print 'choice',choice -# -# if choice==0: # pick the single gaussian fit -# Amplitude0=lsq1[0][2] -# shift0=lsq1[0][0] -# width0=lsq1[0][1] -# p0=lsq1[0][3] -# Amplitude1=0. -# shift1=0. -# width1=0. -# p1=0. -# noise=lsq1[0][4] -# elif choice==1: # take the first one of the 2 gaussians fitted -# Amplitude0 = lsq2[0][2] -# shift0 = lsq2[0][0] -# width0 = lsq2[0][1] -# p0 = lsq2[0][3] -# Amplitude1 = lsq2[0][6] # This is 0 in gg1 -# shift1 = lsq2[0][4] # This is 0 in gg1 -# width1 = lsq2[0][5] # This is 0 in gg1 -# p1 = lsq2[0][7] # This is 0 in gg1 -# noise = lsq2[0][8] -# else: # the second one -# Amplitude0 = lsq2[0][6] -# shift0 = lsq2[0][4] -# width0 = lsq2[0][5] -# p0 = lsq2[0][7] -# Amplitude1 = lsq2[0][2] # This is 0 in gg1 -# shift1 = lsq2[0][0] # This is 0 in gg1 -# width1 = lsq2[0][1] # This is 0 in gg1 -# p1 = lsq2[0][3] # This is 0 in gg1 -# noise = lsq2[0][8] -# -# #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) -# SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 -# SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 -# #print 'SPC_ch1.shape',SPC_ch1.shape -# #print 'SPC_ch2.shape',SPC_ch2.shape -# #dataOut.data_param = SPC_ch1 -# GauSPC[0] = SPC_ch1 -# GauSPC[1] = SPC_ch2 - -# #plt.gcf().clear() -# plt.figure(50+self.i) -# self.i=self.i+1 -# #plt.subplot(121) -# plt.plot(self.spc,'k')#,label='spc(66)') -# plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1') -# #plt.plot(SPC_ch2,'r')#,label='gg2') -# #plt.plot(xFrec,ySamples[1],'g',label='Ch1') -# #plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# #plt.plot(xFrec,FitGauss,'yo:',label='fit') -# plt.legend() -# plt.title('DATOS A ALTURA DE 7500 METROS') -# plt.show() -# print 'shift0', shift0 -# print 'Amplitude0', Amplitude0 -# print 'width0', width0 -# print 'p0', p0 -# print '========================' -# print 'shift1', shift1 -# print 'Amplitude1', Amplitude1 -# print 'width1', width1 -# print 'p1', p1 -# print 'noise', noise -# print 's_noise', wnoise - - print('========================================================') - print('total_time: ', time.time()-start_time) - - # re-normalizing spc and noise - # This part differs from gg1 - - + dataOut.SPCparam = numpy.asarray(SPCparam) ''' Parameters: 1. Amplitude @@ -524,16 +345,11 @@ class GaussianFit(Operation): 4. Power ''' - - ############################################################################### def FitGau(self, X): Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X - #print 'VARSSSS', ch, pnoise, noise, num_intg - - #print 'HEIGHTS', self.Num_Hei - - GauSPC = [] + + SPCparam = [] SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch1[:] = 0#numpy.NaN @@ -542,10 +358,6 @@ class GaussianFit(Operation): for ht in range(self.Num_Hei): - #print (numpy.asarray(self.spc).shape) - - #print 'TTTTT', ch , ht - #print self.spc.shape spc = numpy.asarray(self.spc)[ch,:,ht] @@ -554,27 +366,26 @@ class GaussianFit(Operation): # normalizing spc and noise # This part differs from gg1 spc_norm_max = max(spc) - spc = spc / spc_norm_max - pnoise = pnoise / spc_norm_max + #spc = spc / spc_norm_max + pnoise = pnoise #/ spc_norm_max ############################################# - + fatspectra=1.0 - wnoise = noise_ / spc_norm_max + wnoise = noise_ #/ spc_norm_max #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used #if wnoise>1.1*pnoise: # to be tested later # wnoise=pnoise - noisebl=wnoise*0.9; noisebh=wnoise*1.1 + noisebl=wnoise*0.9; + noisebh=wnoise*1.1 spc=spc-wnoise - # print 'wnoise', noise_[0], spc_norm_max, wnoise + minx=numpy.argmin(spc) + #spcs=spc.copy() spcs=numpy.roll(spc,-minx) cum=numpy.cumsum(spcs) tot_noise=wnoise * self.Num_Bin #64; - #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise - #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' - #snr=tot_signal/tot_noise - #snr=cum[-1]/tot_noise + snr = sum(spcs)/tot_noise snrdB=10.*numpy.log10(snr) @@ -582,16 +393,15 @@ class GaussianFit(Operation): snr = numpy.NaN SPC_ch1[:,ht] = 0#numpy.NaN SPC_ch1[:,ht] = 0#numpy.NaN - GauSPC = (SPC_ch1,SPC_ch2) + SPCparam = (SPC_ch1,SPC_ch2) continue - #print 'snr',snrdB #, sum(spcs) , tot_noise - #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None - cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region + cummax=max(cum); + epsi=0.08*fatspectra # cumsum to narrow down the energy region cumlo=cummax*epsi; cumhi=cummax*(1-epsi) powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-6: # when SNR is strong pick the peak with least shift (LOS velocity) error + + if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error if oneG: choice=0 else: @@ -690,7 +483,7 @@ class GaussianFit(Operation): s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling - + if gp1>gp2: if a1>0.7*a2: choice=1 @@ -710,13 +503,15 @@ class GaussianFit(Operation): choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) - shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0]) - shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0]) + shift0=lsq2[0][0]; + vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0]) + shift1=lsq2[0][4]; + vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0]) - max_vel = 20 + max_vel = 1.0 #first peak will be 0, second peak will be 1 - if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range + if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range shift0=lsq2[0][0] width0=lsq2[0][1] Amplitude0=lsq2[0][2] @@ -739,120 +534,18 @@ class GaussianFit(Operation): p0=lsq2[0][7] noise=lsq2[0][8] - if Amplitude0<0.1: # in case the peak is noise - shift0,width0,Amplitude0,p0 = 4*[numpy.NaN] - if Amplitude1<0.1: - shift1,width1,Amplitude1,p1 = 4*[numpy.NaN] - - -# if choice==0: # pick the single gaussian fit -# Amplitude0=lsq1[0][2] -# shift0=lsq1[0][0] -# width0=lsq1[0][1] -# p0=lsq1[0][3] -# Amplitude1=0. -# shift1=0. -# width1=0. -# p1=0. -# noise=lsq1[0][4] -# elif choice==1: # take the first one of the 2 gaussians fitted -# Amplitude0 = lsq2[0][2] -# shift0 = lsq2[0][0] -# width0 = lsq2[0][1] -# p0 = lsq2[0][3] -# Amplitude1 = lsq2[0][6] # This is 0 in gg1 -# shift1 = lsq2[0][4] # This is 0 in gg1 -# width1 = lsq2[0][5] # This is 0 in gg1 -# p1 = lsq2[0][7] # This is 0 in gg1 -# noise = lsq2[0][8] -# else: # the second one -# Amplitude0 = lsq2[0][6] -# shift0 = lsq2[0][4] -# width0 = lsq2[0][5] -# p0 = lsq2[0][7] -# Amplitude1 = lsq2[0][2] # This is 0 in gg1 -# shift1 = lsq2[0][0] # This is 0 in gg1 -# width1 = lsq2[0][1] # This is 0 in gg1 -# p1 = lsq2[0][3] # This is 0 in gg1 -# noise = lsq2[0][8] - - #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) + if Amplitude0<0.05: # in case the peak is noise + shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN] + if Amplitude1<0.05: + shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN] + + SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 - #print 'SPC_ch1.shape',SPC_ch1.shape - #print 'SPC_ch2.shape',SPC_ch2.shape - #dataOut.data_param = SPC_ch1 - GauSPC = (SPC_ch1,SPC_ch2) - #GauSPC[1] = SPC_ch2 - -# print 'shift0', shift0 -# print 'Amplitude0', Amplitude0 -# print 'width0', width0 -# print 'p0', p0 -# print '========================' -# print 'shift1', shift1 -# print 'Amplitude1', Amplitude1 -# print 'width1', width1 -# print 'p1', p1 -# print 'noise', noise -# print 's_noise', wnoise + SPCparam = (SPC_ch1,SPC_ch2) - return GauSPC - - - def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination. - y_model=self.y_model1(x,state) - s0,w0,a0,p0,n=state - e0=((x-s0)/w0)**2; - - e0u=((x-s0-self.Num_Bin)/w0)**2; - - e0d=((x-s0+self.Num_Bin)/w0)**2 - m0=numpy.exp(-0.5*e0**(p0/2.)); - m0u=numpy.exp(-0.5*e0u**(p0/2.)); - m0d=numpy.exp(-0.5*e0d**(p0/2.)) - JA=m0+m0u+m0d - JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d) - - JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0) - - JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2 - jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model]) - return jack1.T - - def y_jacobian2(self,x,state): - y_model=self.y_model2(x,state) - s0,w0,a0,p0,s1,w1,a1,p1,n=state - e0=((x-s0)/w0)**2; - - e0u=((x-s0- self.Num_Bin )/w0)**2; - - e0d=((x-s0+ self.Num_Bin )/w0)**2 - e1=((x-s1)/w1)**2; - e1u=((x-s1- self.Num_Bin )/w1)**2; - - e1d=((x-s1+ self.Num_Bin )/w1)**2 - m0=numpy.exp(-0.5*e0**(p0/2.)); - m0u=numpy.exp(-0.5*e0u**(p0/2.)); - m0d=numpy.exp(-0.5*e0d**(p0/2.)) - m1=numpy.exp(-0.5*e1**(p1/2.)); - m1u=numpy.exp(-0.5*e1u**(p1/2.)); - m1d=numpy.exp(-0.5*e1d**(p1/2.)) - JA=m0+m0u+m0d - JA1=m1+m1u+m1d - JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d) - JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d) - - JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0) - - JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1) - - JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2 - - JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2 - jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model]) - return jack2.T + return GauSPC def y_model1(self,x,state): shift0,width0,amplitude0,power0,noise=state @@ -884,6 +577,7 @@ class GaussianFit(Operation): def misfit2(self,state,y_data,x,num_intg): return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.) + class PrecipitationProc(Operation): @@ -900,24 +594,61 @@ class PrecipitationProc(Operation): Parameters affected: ''' - - def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None, - tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None): + def __init__(self): + Operation.__init__(self) + self.i=0 + + + def gaus(self,xSamples,Amp,Mu,Sigma): + return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) + + + + def Moments(self, ySamples, xSamples): + Pot = numpy.nansum( ySamples ) # Potencia, momento 0 + yNorm = ySamples / Pot + + Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento + Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento + Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral + + return numpy.array([Pot, Vr, Desv]) + + def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118, + tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350): - self.spc = dataOut.data_pre[0].copy() - self.Num_Hei = self.spc.shape[2] - self.Num_Bin = self.spc.shape[1] - self.Num_Chn = self.spc.shape[0] - Velrange = dataOut.abscissaList + Velrange = dataOut.spcparam_range[2] + FrecRange = dataOut.spcparam_range[0] + + dV= Velrange[1]-Velrange[0] + dF= FrecRange[1]-FrecRange[0] if radar == "MIRA35C" : + self.spc = dataOut.data_pre[0].copy() + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] Ze = self.dBZeMODE2(dataOut) else: + self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() # + + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" + + self.spc[:,:,0:7]= numpy.NaN + + """##########################################""" + + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] + + ''' Se obtiene la constante del RADAR ''' + self.Pt = Pt self.Gt = Gt self.Gr = Gr @@ -927,48 +658,101 @@ class PrecipitationProc(Operation): self.ThetaT = ThetaT self.ThetaR = ThetaR - RadarConstant = GetRadarConstant() - SPCmean = numpy.mean(self.spc,0) - ETA = numpy.zeros(self.Num_Hei) - Pr = numpy.sum(SPCmean,0) + Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) + Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR) + RadarConstant = 5e-26 * Numerator / Denominator # - #for R in range(self.Num_Hei): - # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) - - D_range = numpy.zeros(self.Num_Hei) - EqSec = numpy.zeros(self.Num_Hei) + ''' ============================= ''' + + self.spc[0] = (self.spc[0]-dataOut.noise[0]) + self.spc[1] = (self.spc[1]-dataOut.noise[1]) + self.spc[2] = (self.spc[2]-dataOut.noise[2]) + + self.spc[ numpy.where(self.spc < 0)] = 0 + + SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise)) + SPCmean[ numpy.where(SPCmean < 0)] = 0 + + ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei]) + ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei]) + ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei]) + + Pr = SPCmean[:,:] + + VelMeteoro = numpy.mean(SPCmean,axis=0) + + D_range = numpy.zeros([self.Num_Bin,self.Num_Hei]) + SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei]) + N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei]) + V_mean = numpy.zeros(self.Num_Hei) del_V = numpy.zeros(self.Num_Hei) + Z = numpy.zeros(self.Num_Hei) + Ze = numpy.zeros(self.Num_Hei) + RR = numpy.zeros(self.Num_Hei) + + Range = dataOut.heightList*1000. for R in range(self.Num_Hei): - ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) - h = R + Altitude #Range from ground to radar pulse altitude + h = Range[R] + Altitude #Range from ground to radar pulse altitude del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity - D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity - SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma) + D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3 + + '''NOTA: ETA(n) dn = ETA(f) df + + dn = 1 Diferencial de muestreo + df = ETA(n) / ETA(f) + + ''' + + ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA) + + ETAv[:,R]=ETAn[:,R]/dV + + ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R]) + + SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma) + + N_dist[:,R] = ETAn[:,R] / SIGMA[:,R] + + DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin]) + + try: + popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments) + except: + popt01=numpy.zeros(3) + popt01[1]= DMoments[1] + + if popt01[1]<0 or popt01[1]>20: + popt01[1]=numpy.NaN + + + V_mean[R]=popt01[1] + + Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18 + + RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate + + Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km) - N_dist[R] = ETA[R] / SIGMA[R] - - Ze = (ETA * Lambda**4) / (numpy.pi * Km) - Z = numpy.sum( N_dist * D_range**6 ) - RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate - RR = (Ze/200)**(1/1.6) + RR2 = (Z/200)**(1/1.6) dBRR = 10*numpy.log10(RR) + dBRR2 = 10*numpy.log10(RR2) dBZe = 10*numpy.log10(Ze) - dataOut.data_output = Ze - dataOut.data_param = numpy.ones([2,self.Num_Hei]) - dataOut.channelList = [0,1] - print('channelList', dataOut.channelList) - dataOut.data_param[0]=dBZe - dataOut.data_param[1]=dBRR - print('RR SHAPE', dBRR.shape) - print('Ze SHAPE', dBZe.shape) - print('dataOut.data_param SHAPE', dataOut.data_param.shape) - + dBZ = 10*numpy.log10(Z) + + dataOut.data_output = RR[8] + dataOut.data_param = numpy.ones([3,self.Num_Hei]) + dataOut.channelList = [0,1,2] + + dataOut.data_param[0]=dBZ + dataOut.data_param[1]=V_mean + dataOut.data_param[2]=RR + def dBZeMODE2(self, dataOut): # Processing for MIRA35C @@ -983,7 +767,7 @@ class PrecipitationProc(Operation): data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN ETA = numpy.sum(SNR,1) - print('ETA' , ETA) + ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN) Ze = numpy.ones([self.Num_Chn, self.Num_Hei] ) @@ -995,26 +779,27 @@ class PrecipitationProc(Operation): return Ze - def GetRadarConstant(self): - - """ - Constants: - - Pt: Transmission Power dB - Gt: Transmission Gain dB - Gr: Reception Gain dB - Lambda: Wavelenght m - aL: Attenuation loses dB - tauW: Width of transmission pulse s - ThetaT: Transmission antenna bean angle rad - ThetaR: Reception antenna beam angle rad - - """ - Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) - Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) - RadarConstant = Numerator / Denominator - - return RadarConstant +# def GetRadarConstant(self): +# +# """ +# Constants: +# +# Pt: Transmission Power dB 5kW 5000 +# Gt: Transmission Gain dB 24.7 dB 295.1209 +# Gr: Reception Gain dB 18.5 dB 70.7945 +# Lambda: Wavelenght m 0.6741 m 0.6741 +# aL: Attenuation loses dB 4dB 2.5118 +# tauW: Width of transmission pulse s 4us 4e-6 +# ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317 +# ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087 +# +# """ +# +# Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) +# Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) +# RadarConstant = Numerator / Denominator +# +# return RadarConstant @@ -1037,10 +822,20 @@ class FullSpectralAnalysis(Operation): Parameters affected: Winds, height range, SNR """ - def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7): + def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7): + + self.indice=int(numpy.random.rand()*1000) spc = dataOut.data_pre[0].copy() - cspc = dataOut.data_pre[1].copy() + cspc = dataOut.data_pre[1] + + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" + + SNRspc = spc.copy() + SNRspc[:,:,0:7]= numpy.NaN + + """##########################################""" + nChannel = spc.shape[0] nProfiles = spc.shape[1] @@ -1050,14 +845,9 @@ class FullSpectralAnalysis(Operation): if dataOut.ChanDist is not None : ChanDist = dataOut.ChanDist else: - ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]]) - - #print 'ChanDist', ChanDist + ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]]) - if dataOut.VelRange is not None: - VelRange= dataOut.VelRange - else: - VelRange= dataOut.abscissaList + FrecRange = dataOut.spc_range[0] ySamples=numpy.ones([nChannel,nProfiles]) phase=numpy.ones([nChannel,nProfiles]) @@ -1065,113 +855,108 @@ class FullSpectralAnalysis(Operation): coherence=numpy.ones([nChannel,nProfiles]) PhaseSlope=numpy.ones(nChannel) PhaseInter=numpy.ones(nChannel) - dataSNR = dataOut.data_SNR - - + data_SNR=numpy.zeros([nProfiles]) data = dataOut.data_pre noise = dataOut.noise - print('noise',noise) - #SNRdB = 10*numpy.log10(dataOut.data_SNR) - FirstMoment = numpy.average(dataOut.data_param[:,1,:],0) - #SNRdBMean = [] - + dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0] - #for j in range(nHeights): - # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]])) - # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]])) - - data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN + dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20 + + + data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN velocityX=[] velocityY=[] velocityV=[] + PhaseLine=[] - dbSNR = 10*numpy.log10(dataSNR) + dbSNR = 10*numpy.log10(dataOut.data_SNR) dbSNR = numpy.average(dbSNR,0) + for Height in range(nHeights): - - [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit) + + [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit) + PhaseLine = numpy.append(PhaseLine, PhaseSlope) if abs(Vzon)<100. and abs(Vzon)> 0.: velocityX=numpy.append(velocityX, Vzon)#Vmag else: - print('Vzon',Vzon) velocityX=numpy.append(velocityX, numpy.NaN) if abs(Vmer)<100. and abs(Vmer) > 0.: - velocityY=numpy.append(velocityY, Vmer)#Vang + velocityY=numpy.append(velocityY, -Vmer)#Vang else: - print('Vmer',Vmer) velocityY=numpy.append(velocityY, numpy.NaN) if dbSNR[Height] > SNRlimit: - velocityV=numpy.append(velocityV, FirstMoment[Height]) + velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height]) else: velocityV=numpy.append(velocityV, numpy.NaN) - #FirstMoment[Height]= numpy.NaN -# if SNRdBMean[Height] <12: -# FirstMoment[Height] = numpy.NaN -# velocityX[Height] = numpy.NaN -# velocityY[Height] = numpy.NaN - - - data_output[0]=numpy.array(velocityX) - data_output[1]=numpy.array(velocityY) - data_output[2]=-velocityV#FirstMoment - - print(' ') - #print 'FirstMoment' - #print FirstMoment - print('velocityX',data_output[0]) - print(' ') - print('velocityY',data_output[1]) - #print numpy.array(velocityY) - print(' ') - #print 'SNR' - #print 10*numpy.log10(dataOut.data_SNR) - #print numpy.shape(10*numpy.log10(dataOut.data_SNR)) - print(' ') + + + '''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR''' + data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1) + data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1) + data_output[2] = velocityV#FirstMoment + + xFrec=FrecRange[0:spc.shape[1]] dataOut.data_output=data_output + return def moving_average(self,x, N=2): return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] - def gaus(self,xSamples,a,x0,sigma): - return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2)) + def gaus(self,xSamples,Amp,Mu,Sigma): + return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) + + - def Find(self,x,value): - for index in range(len(x)): - if x[index]==value: - return index + def Moments(self, ySamples, xSamples): + Pot = numpy.nansum( ySamples ) # Potencia, momento 0 + yNorm = ySamples / Pot + Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento + Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento + Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral + + return numpy.array([Pot, Vr, Desv]) - def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit): + def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit): + + ySamples=numpy.ones([spc.shape[0],spc.shape[1]]) phase=numpy.ones([spc.shape[0],spc.shape[1]]) CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) coherence=numpy.ones([spc.shape[0],spc.shape[1]]) - PhaseSlope=numpy.ones(spc.shape[0]) + PhaseSlope=numpy.zeros(spc.shape[0]) PhaseInter=numpy.ones(spc.shape[0]) - xFrec=VelRange + xFrec=AbbsisaRange[0][0:spc.shape[1]] + xVel =AbbsisaRange[2][0:spc.shape[1]] + Vv=numpy.empty(spc.shape[2])*0 + SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]# + + SPCmoments = self.Moments(SPCav[:,Height], xVel ) + CSPCmoments = [] + cspcNoise = numpy.empty(3) '''Getting Eij and Nij''' - E01=ChanDist[0][0] - N01=ChanDist[0][1] + Xi01=ChanDist[0][0] + Eta01=ChanDist[0][1] - E02=ChanDist[1][0] - N02=ChanDist[1][1] + Xi02=ChanDist[1][0] + Eta02=ChanDist[1][1] - E12=ChanDist[2][0] - N12=ChanDist[2][1] + Xi12=ChanDist[2][0] + Eta12=ChanDist[2][1] z = spc.copy() z = numpy.where(numpy.isfinite(z), z, numpy.NAN) @@ -1179,176 +964,197 @@ class FullSpectralAnalysis(Operation): for i in range(spc.shape[0]): '''****** Line of Data SPC ******''' - zline=z[i,:,Height] + zline=z[i,:,Height].copy() - noise[i] # Se resta ruido '''****** SPC is normalized ******''' - FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy()) - FactNorm= FactNorm/numpy.sum(FactNorm) + SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido + FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado - SmoothSPC=self.moving_average(FactNorm,N=3) - - xSamples = ar(list(range(len(SmoothSPC)))) - ySamples[i] = SmoothSPC - - #dbSNR=10*numpy.log10(dataSNR) - print(' ') - print(' ') - print(' ') - - #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120] - print('SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20]) - print('noise',noise) - print('zline',zline.shape, zline[0:20]) - print('FactNorm',FactNorm.shape, FactNorm[0:20]) - print('FactNorm suma', numpy.sum(FactNorm)) + xSamples = xFrec # Se toma el rango de frecuncias + ySamples[i] = FactNorm # Se toman los valores de SPC normalizado for i in range(spc.shape[0]): '''****** Line of Data CSPC ******''' - cspcLine=cspc[i,:,Height].copy() + cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido + SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido + cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado - '''****** CSPC is normalized ******''' + '''****** CSPC is normalized with respect to Briggs and Vincent ******''' chan_index0 = pairsList[i][0] chan_index1 = pairsList[i][1] - CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) # - CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor) + CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2 + CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor) CSPCSamples[i] = CSPCNorm + coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor) - coherence[i]= self.moving_average(coherence[i],N=2) + #coherence[i]= self.moving_average(coherence[i],N=1) phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi - print('cspcLine', cspcLine.shape, cspcLine[0:20]) - print('CSPCFactor', CSPCFactor)#, CSPCFactor[0:20] - print(numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i]) - print('CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20]) - print('CSPCNorm suma', numpy.sum(CSPCNorm)) - print('CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20]) + CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples), + self.Moments(numpy.abs(CSPCSamples[1]), xSamples), + self.Moments(numpy.abs(CSPCSamples[2]), xSamples)]) - '''****** Getting fij width ******''' + + popt=[1e-10,0,1e-10] + popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10] + FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0 + + CSPCMask01 = numpy.abs(CSPCSamples[0]) + CSPCMask02 = numpy.abs(CSPCSamples[1]) + CSPCMask12 = numpy.abs(CSPCSamples[2]) + + mask01 = ~numpy.isnan(CSPCMask01) + mask02 = ~numpy.isnan(CSPCMask02) + mask12 = ~numpy.isnan(CSPCMask12) - yMean=[] - yMean2=[] + #mask = ~numpy.isnan(CSPCMask01) + CSPCMask01 = CSPCMask01[mask01] + CSPCMask02 = CSPCMask02[mask02] + CSPCMask12 = CSPCMask12[mask12] + #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01) - for j in range(len(ySamples[1])): - yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]])) - '''******* Getting fitting Gaussian ******''' - meanGauss=sum(xSamples*yMean) / len(xSamples) - sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) - print('****************************') - print('len(xSamples): ',len(xSamples)) - print('yMean: ', yMean.shape, yMean[0:20]) - print('ySamples', ySamples.shape, ySamples[0,0:20]) - print('xSamples: ',xSamples.shape, xSamples[0:20]) + '''***Fit Gauss CSPC01***''' + if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 : + try: + popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0]) + popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1]) + popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2]) + FitGauss01 = self.gaus(xSamples,*popt01) + FitGauss02 = self.gaus(xSamples,*popt02) + FitGauss12 = self.gaus(xSamples,*popt12) + except: + FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0])) + FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1])) + FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2])) + + + CSPCopt = numpy.vstack([popt01,popt02,popt12]) + + '''****** Getting fij width ******''' + + yMean = numpy.average(ySamples, axis=0) # ySamples[0] + + '''******* Getting fitting Gaussian *******''' + meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia) + sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia) - print('meanGauss',meanGauss) - print('sigma',sigma) + yMoments = self.Moments(yMean, xSamples) - #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001): - if dbSNR > SNRlimit : - try: - popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) + if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001: + try: + popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments) + FitGauss=self.gaus(xSamples,*popt) - if numpy.amax(popt)>numpy.amax(yMean)*0.3: - FitGauss=self.gaus(xSamples,*popt) - - else: - FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) - print('Verificador: Dentro', Height) except :#RuntimeError: FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) - + else: FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) - Maximun=numpy.amax(yMean) - eMinus1=Maximun*numpy.exp(-1)#*0.8 - HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) - HalfWidth= xFrec[HWpos] - GCpos=self.Find(FitGauss, numpy.amax(FitGauss)) - Vpos=self.Find(FactNorm, numpy.amax(FactNorm)) - - #Vpos=FirstMoment[] '''****** Getting Fij ******''' + Fijcspc = CSPCopt[:,2]/2*3 + + + GaussCenter = popt[1] #xFrec[GCpos] + #Punto en Eje X de la Gaussiana donde se encuentra el centro + ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()] + PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0] + + #Punto e^-1 hubicado en la Gaussiana + PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1) + FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss" + PointFij = numpy.where(FitGauss==FijClosest)[0][0] - GaussCenter=xFrec[GCpos] - if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): - Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 + if xSamples[PointFij] > xSamples[PointGauCenter]: + Fij = xSamples[PointFij] - xSamples[PointGauCenter] + else: - Fij=abs(GaussCenter-HalfWidth)+0.0000001 + Fij = xSamples[PointGauCenter] - xSamples[PointFij] + + + '''****** Taking frequency ranges from SPCs ******''' - '''****** Getting Frecuency range of significant data ******''' - Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) + #GaussCenter = popt[1] #Primer momento 01 + GauWidth = popt[2] *3/2 #Ancho de banda de Gau01 + Range = numpy.empty(2) + Range[0] = GaussCenter - GauWidth + Range[1] = GaussCenter + GauWidth + #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max) + ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()] + ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()] - if Rangpos5 and len(FrecRange)5 and len(FrecRange)4: + Vver=popt[1] + else: + Vver=numpy.NaN + FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12]) + + + return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC + class SpectralMoments(Operation): ''' @@ -1384,7 +1195,7 @@ class SpectralMoments(Operation): self.dataOut.noise : Noise level per channel Affected: - self.dataOut.data_param : Parameters per channel + self.dataOut.moments : Parameters per channel self.dataOut.data_SNR : SNR per channel ''' @@ -1401,7 +1212,7 @@ class SpectralMoments(Operation): for ind in range(nChannel): data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] ) - dataOut.data_param = data_param[:,1:,:] + dataOut.moments = data_param[:,1:,:] dataOut.data_SNR = data_param[:,0] dataOut.data_DOP = data_param[:,1] dataOut.data_MEAN = data_param[:,2] @@ -1431,6 +1242,8 @@ class SpectralMoments(Operation): vec_fd = numpy.zeros(oldspec.shape[1]) vec_w = numpy.zeros(oldspec.shape[1]) vec_snr = numpy.zeros(oldspec.shape[1]) + + oldspec = numpy.ma.masked_invalid(oldspec) for ind in range(oldspec.shape[1]): @@ -1469,7 +1282,7 @@ class SpectralMoments(Operation): fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power) snr = (spec2.mean()-n0)/n0 - + if (snr < 1.e-20) : snr = 1.e-20 @@ -1477,7 +1290,7 @@ class SpectralMoments(Operation): vec_fd[ind] = fd vec_w[ind] = w vec_snr[ind] = snr - + moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w)) return moments @@ -1675,7 +1488,6 @@ class SpectralFitting(Operation): dataCross = dataCross**2/K for h in range(nHeights): -# print self.dataOut.heightList[h] #Input d = data[:,h] @@ -1734,7 +1546,7 @@ class SpectralFitting(Operation): fm = self.dataOut.library.modelFunction(p, constants) fmp=numpy.dot(LT,fm) - + return dp-fmp def __getSNR(self, z, noise): @@ -1768,8 +1580,8 @@ class WindProfiler(Operation): n = None - def __init__(self, **kwargs): - Operation.__init__(self, **kwargs) + def __init__(self): + Operation.__init__(self) def __calculateCosDir(self, elev, azim): zen = (90 - elev)*numpy.pi/180 @@ -2071,12 +1883,9 @@ class WindProfiler(Operation): Parameters affected: Winds ''' -# print arrayMeteor.shape #Settings nInt = (heightMax - heightMin)/2 -# print nInt nInt = int(nInt) -# print nInt winds = numpy.zeros((2,nInt))*numpy.nan #Filter errors @@ -2475,8 +2284,8 @@ class WindProfiler(Operation): class EWDriftsEstimation(Operation): - def __init__(self, **kwargs): - Operation.__init__(self, **kwargs) + def __init__(self): + Operation.__init__(self) def __correctValues(self, heiRang, phi, velRadial, SNR): listPhi = phi.tolist() diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index db71e0d..5c68c91 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -159,9 +159,7 @@ class SpectraProc(ProcessingUnit): dtype='complex') if self.dataIn.flagDataAsBlock: - # data dimension: [nChannels, nProfiles, nSamples] nVoltProfiles = self.dataIn.data.shape[1] - # nVoltProfiles = self.dataIn.nProfiles if nVoltProfiles == nProfiles: self.buffer = self.dataIn.data.copy() @@ -299,7 +297,57 @@ class SpectraProc(ProcessingUnit): self.__selectPairsByChannel(self.dataOut.channelList) return 1 + + + def selectFFTs(self, minFFT, maxFFT ): + """ + Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango + minFFT<= FFT <= maxFFT + """ + + if (minFFT > maxFFT): + raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT)) + + if (minFFT < self.dataOut.getFreqRange()[0]): + minFFT = self.dataOut.getFreqRange()[0] + + if (maxFFT > self.dataOut.getFreqRange()[-1]): + maxFFT = self.dataOut.getFreqRange()[-1] + + minIndex = 0 + maxIndex = 0 + FFTs = self.dataOut.getFreqRange() + + inda = numpy.where(FFTs >= minFFT) + indb = numpy.where(FFTs <= maxFFT) + + try: + minIndex = inda[0][0] + except: + minIndex = 0 + + try: + maxIndex = indb[0][-1] + except: + maxIndex = len(FFTs) + + self.selectFFTsByIndex(minIndex, maxIndex) + return 1 + + + def setH0(self, h0, deltaHeight = None): + + if not deltaHeight: + deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0] + + nHeights = self.dataOut.nHeights + + newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight + + self.dataOut.heightList = newHeiRange + + def selectHeights(self, minHei, maxHei): """ Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango @@ -316,9 +364,9 @@ class SpectraProc(ProcessingUnit): 1 si el metodo se ejecuto con exito caso contrario devuelve 0 """ + if (minHei > maxHei): - raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % ( - minHei, maxHei)) + raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)) if (minHei < self.dataOut.heightList[0]): minHei = self.dataOut.heightList[0] @@ -344,6 +392,7 @@ class SpectraProc(ProcessingUnit): maxIndex = len(heights) self.selectHeightsByIndex(minIndex, maxIndex) + return 1 @@ -389,6 +438,40 @@ class SpectraProc(ProcessingUnit): return 1 + def selectFFTsByIndex(self, minIndex, maxIndex): + """ + + """ + + if (minIndex < 0) or (minIndex > maxIndex): + raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)) + + if (maxIndex >= self.dataOut.nProfiles): + maxIndex = self.dataOut.nProfiles-1 + + #Spectra + data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:] + + data_cspc = None + if self.dataOut.data_cspc is not None: + data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:] + + data_dc = None + if self.dataOut.data_dc is not None: + data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:] + + self.dataOut.data_spc = data_spc + self.dataOut.data_cspc = data_cspc + self.dataOut.data_dc = data_dc + + self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1]) + self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1] + self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1] + + return 1 + + + def selectHeightsByIndex(self, minIndex, maxIndex): """ Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango @@ -494,7 +577,32 @@ class SpectraProc(ProcessingUnit): return 1 - def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None): + def removeInterference2(self): + + cspc = self.dataOut.data_cspc + spc = self.dataOut.data_spc + Heights = numpy.arange(cspc.shape[2]) + realCspc = numpy.abs(cspc) + + for i in range(cspc.shape[0]): + LinePower= numpy.sum(realCspc[i], axis=0) + Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)] + SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ] + InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 ) + InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)] + InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)] + + + InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) ) + #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax])) + if len(InterferenceRange) datetime.timedelta(seconds=0): - return - - self.time_of_last_call = datetime.datetime.now() - return fn(*args, **kwargs) - - return wrapper - -class Data(object): - ''' - Object to hold data to be plotted - ''' - - def __init__(self, plottypes, throttle_value, exp_code, buffering=True): - self.plottypes = plottypes - self.throttle = throttle_value - self.exp_code = exp_code - self.buffering = buffering - self.ended = False - self.localtime = False - self.meta = {} - self.__times = [] - self.__heights = [] - - def __str__(self): - dum = ['{}{}'.format(key, self.shape(key)) for key in self.data] - return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times)) - - def __len__(self): - return len(self.__times) - - def __getitem__(self, key): - if key not in self.data: - raise KeyError(log.error('Missing key: {}'.format(key))) - - if 'spc' in key or not self.buffering: - ret = self.data[key] - else: - ret = numpy.array([self.data[key][x] for x in self.times]) - if ret.ndim > 1: - ret = numpy.swapaxes(ret, 0, 1) - return ret - - def __contains__(self, key): - return key in self.data - - def setup(self): - ''' - Configure object - ''' - - self.type = '' - self.ended = False - self.data = {} - self.__times = [] - self.__heights = [] - self.__all_heights = set() - for plot in self.plottypes: - if 'snr' in plot: - plot = 'snr' - self.data[plot] = {} - - def shape(self, key): - ''' - Get the shape of the one-element data for the given key - ''' - - if len(self.data[key]): - if 'spc' in key or not self.buffering: - return self.data[key].shape - return self.data[key][self.__times[0]].shape - return (0,) - - def update(self, dataOut, tm): - ''' - Update data object with new dataOut - ''' - - if tm in self.__times: - return - - self.type = dataOut.type - self.parameters = getattr(dataOut, 'parameters', []) - if hasattr(dataOut, 'pairsList'): - self.pairs = dataOut.pairsList - if hasattr(dataOut, 'meta'): - self.meta = dataOut.meta - self.channels = dataOut.channelList - self.interval = dataOut.getTimeInterval() - self.localtime = dataOut.useLocalTime - if 'spc' in self.plottypes or 'cspc' in self.plottypes: - self.xrange = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) - self.__heights.append(dataOut.heightList) - self.__all_heights.update(dataOut.heightList) - self.__times.append(tm) - - for plot in self.plottypes: - if plot == 'spc': - z = dataOut.data_spc/dataOut.normFactor - buffer = 10*numpy.log10(z) - if plot == 'cspc': - buffer = dataOut.data_cspc - if plot == 'noise': - buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) - if plot == 'rti': - buffer = dataOut.getPower() - if plot == 'snr_db': - buffer = dataOut.data_SNR - if plot == 'snr': - buffer = 10*numpy.log10(dataOut.data_SNR) - if plot == 'dop': - buffer = 10*numpy.log10(dataOut.data_DOP) - if plot == 'mean': - buffer = dataOut.data_MEAN - if plot == 'std': - buffer = dataOut.data_STD - if plot == 'coh': - buffer = dataOut.getCoherence() - if plot == 'phase': - buffer = dataOut.getCoherence(phase=True) - if plot == 'output': - buffer = dataOut.data_output - if plot == 'param': - buffer = dataOut.data_param - - if 'spc' in plot: - self.data[plot] = buffer - else: - if self.buffering: - self.data[plot][tm] = buffer - else: - self.data[plot] = buffer - - def normalize_heights(self): - ''' - Ensure same-dimension of the data for different heighList - ''' - - H = numpy.array(list(self.__all_heights)) - H.sort() - for key in self.data: - shape = self.shape(key)[:-1] + H.shape - for tm, obj in list(self.data[key].items()): - h = self.__heights[self.__times.index(tm)] - if H.size == h.size: - continue - index = numpy.where(numpy.in1d(H, h))[0] - dummy = numpy.zeros(shape) + numpy.nan - if len(shape) == 2: - dummy[:, index] = obj - else: - dummy[index] = obj - self.data[key][tm] = dummy - - self.__heights = [H for tm in self.__times] - - def jsonify(self, decimate=False): - ''' - Convert data to json - ''' - - data = {} - tm = self.times[-1] - dy = int(self.heights.size/MAXNUMY) + 1 - for key in self.data: - if key in ('spc', 'cspc') or not self.buffering: - dx = int(self.data[key].shape[1]/MAXNUMX) + 1 - data[key] = roundFloats(self.data[key][::, ::dx, ::dy].tolist()) - else: - data[key] = roundFloats(self.data[key][tm].tolist()) - - ret = {'data': data} - ret['exp_code'] = self.exp_code - ret['time'] = tm - ret['interval'] = self.interval - ret['localtime'] = self.localtime - ret['yrange'] = roundFloats(self.heights[::dy].tolist()) - if 'spc' in self.data or 'cspc' in self.data: - ret['xrange'] = roundFloats(self.xrange[2][::dx].tolist()) - else: - ret['xrange'] = [] - if hasattr(self, 'pairs'): - ret['pairs'] = self.pairs - else: - ret['pairs'] = [] - - for key, value in list(self.meta.items()): - ret[key] = value - - return json.dumps(ret) - - @property - def times(self): - ''' - Return the list of times of the current data - ''' - - ret = numpy.array(self.__times) - ret.sort() - return ret - - @property - def heights(self): - ''' - Return the list of heights of the current data - ''' - - return numpy.array(self.__heights[-1]) class PublishData(Operation): ''' Operation to send data over zmq. ''' - __attrs__ = ['host', 'port', 'delay', 'zeromq', 'mqtt', 'verbose'] + __attrs__ = ['host', 'port', 'delay', 'verbose'] def __init__(self, **kwargs): """Inicio.""" Operation.__init__(self, **kwargs) self.isConfig = False - self.client = None - self.zeromq = None - self.mqtt = None - def on_disconnect(self, client, userdata, rc): - if rc != 0: - log.warning('Unexpected disconnection.') - self.connect() - - def connect(self): - log.warning('trying to connect') - try: - self.client.connect( - host=self.host, - port=self.port, - keepalive=60*10, - bind_address='') - self.client.loop_start() - # self.client.publish( - # self.topic + 'SETUP', - # json.dumps(setup), - # retain=True - # ) - except: - log.error('MQTT Conection error.') - self.client = False - - def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs): + def setup(self, server='zmq.pipe', delay=0, verbose=True, **kwargs): self.counter = 0 - self.topic = kwargs.get('topic', 'schain') self.delay = kwargs.get('delay', 0) - self.plottype = kwargs.get('plottype', 'spectra') - self.host = kwargs.get('host', "10.10.10.82") - self.port = kwargs.get('port', 3000) - self.clientId = clientId self.cnt = 0 - self.zeromq = zeromq - self.mqtt = kwargs.get('plottype', 0) - self.client = None self.verbose = verbose setup = [] - if mqtt is 1: - self.client = mqtt.Client( - client_id=self.clientId + self.topic + 'SCHAIN', - clean_session=True) - self.client.on_disconnect = self.on_disconnect - self.connect() - for plot in self.plottype: - setup.append({ - 'plot': plot, - 'topic': self.topic + plot, - 'title': getattr(self, plot + '_' + 'title', False), - 'xlabel': getattr(self, plot + '_' + 'xlabel', False), - 'ylabel': getattr(self, plot + '_' + 'ylabel', False), - 'xrange': getattr(self, plot + '_' + 'xrange', False), - 'yrange': getattr(self, plot + '_' + 'yrange', False), - 'zrange': getattr(self, plot + '_' + 'zrange', False), - }) - if zeromq is 1: - context = zmq.Context() - self.zmq_socket = context.socket(zmq.PUSH) - server = kwargs.get('server', 'zmq.pipe') - - if 'tcp://' in server: - address = server - else: - address = 'ipc:///tmp/%s' % server - - self.zmq_socket.connect(address) - time.sleep(1) + context = zmq.Context() + self.zmq_socket = context.socket(zmq.PUSH) + server = kwargs.get('server', 'zmq.pipe') + + if 'tcp://' in server: + address = server + else: + address = 'ipc:///tmp/%s' % server + + self.zmq_socket.connect(address) + time.sleep(1) def publish_data(self): self.dataOut.finished = False - if self.mqtt is 1: - yData = self.dataOut.heightList[:2].tolist() - if self.plottype == 'spectra': - data = getattr(self.dataOut, 'data_spc') - z = data/self.dataOut.normFactor - zdB = 10*numpy.log10(z) - xlen, ylen = zdB[0].shape - dx = int(xlen/MAXNUMX) + 1 - dy = int(ylen/MAXNUMY) + 1 - Z = [0 for i in self.dataOut.channelList] - for i in self.dataOut.channelList: - Z[i] = zdB[i][::dx, ::dy].tolist() - payload = { - 'timestamp': self.dataOut.utctime, - 'data': roundFloats(Z), - 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList], - 'interval': self.dataOut.getTimeInterval(), - 'type': self.plottype, - 'yData': yData - } - - elif self.plottype in ('rti', 'power'): - data = getattr(self.dataOut, 'data_spc') - z = data/self.dataOut.normFactor - avg = numpy.average(z, axis=1) - avgdB = 10*numpy.log10(avg) - xlen, ylen = z[0].shape - dy = numpy.floor(ylen/self.__MAXNUMY) + 1 - AVG = [0 for i in self.dataOut.channelList] - for i in self.dataOut.channelList: - AVG[i] = avgdB[i][::dy].tolist() - payload = { - 'timestamp': self.dataOut.utctime, - 'data': roundFloats(AVG), - 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList], - 'interval': self.dataOut.getTimeInterval(), - 'type': self.plottype, - 'yData': yData - } - elif self.plottype == 'noise': - noise = self.dataOut.getNoise()/self.dataOut.normFactor - noisedB = 10*numpy.log10(noise) - payload = { - 'timestamp': self.dataOut.utctime, - 'data': roundFloats(noisedB.reshape(-1, 1).tolist()), - 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList], - 'interval': self.dataOut.getTimeInterval(), - 'type': self.plottype, - 'yData': yData - } - elif self.plottype == 'snr': - data = getattr(self.dataOut, 'data_SNR') - avgdB = 10*numpy.log10(data) - - ylen = data[0].size - dy = numpy.floor(ylen/self.__MAXNUMY) + 1 - AVG = [0 for i in self.dataOut.channelList] - for i in self.dataOut.channelList: - AVG[i] = avgdB[i][::dy].tolist() - payload = { - 'timestamp': self.dataOut.utctime, - 'data': roundFloats(AVG), - 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList], - 'type': self.plottype, - 'yData': yData - } - else: - print("Tipo de grafico invalido") - payload = { - 'data': 'None', - 'timestamp': 'None', - 'type': None - } - - self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0) - - if self.zeromq is 1: - if self.verbose: - log.log( - 'Sending {} - {}'.format(self.dataOut.type, self.dataOut.datatime), - self.name - ) - self.zmq_socket.send_pyobj(self.dataOut) + + if self.verbose: + log.log( + 'Sending {} - {}'.format(self.dataOut.type, self.dataOut.datatime), + self.name + ) + self.zmq_socket.send_pyobj(self.dataOut) def run(self, dataOut, **kwargs): self.dataOut = dataOut @@ -487,15 +109,12 @@ class PublishData(Operation): time.sleep(self.delay) def close(self): - if self.zeromq is 1: - self.dataOut.finished = True - self.zmq_socket.send_pyobj(self.dataOut) - time.sleep(0.1) - self.zmq_socket.close() - if self.client: - self.client.loop_stop() - self.client.disconnect() - + + self.dataOut.finished = True + self.zmq_socket.send_pyobj(self.dataOut) + time.sleep(0.1) + self.zmq_socket.close() + class ReceiverData(ProcessingUnit): @@ -536,185 +155,6 @@ class ReceiverData(ProcessingUnit): 'Receiving') -class PlotterReceiver(ProcessingUnit, Process): - - throttle_value = 5 - __attrs__ = ['server', 'plottypes', 'realtime', 'localtime', 'throttle', - 'exp_code', 'web_server', 'buffering'] - - def __init__(self, **kwargs): - - ProcessingUnit.__init__(self, **kwargs) - Process.__init__(self) - self.mp = False - self.isConfig = False - self.isWebConfig = False - self.connections = 0 - server = kwargs.get('server', 'zmq.pipe') - web_server = kwargs.get('web_server', None) - if 'tcp://' in server: - address = server - else: - address = 'ipc:///tmp/%s' % server - self.address = address - self.web_address = web_server - self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')] - self.realtime = kwargs.get('realtime', False) - self.localtime = kwargs.get('localtime', True) - self.buffering = kwargs.get('buffering', True) - self.throttle_value = kwargs.get('throttle', 5) - self.exp_code = kwargs.get('exp_code', None) - self.sendData = self.initThrottle(self.throttle_value) - self.dates = [] - self.setup() - - def setup(self): - - self.data = Data(self.plottypes, self.throttle_value, self.exp_code, self.buffering) - self.isConfig = True - - def event_monitor(self, monitor): - - events = {} - - for name in dir(zmq): - if name.startswith('EVENT_'): - value = getattr(zmq, name) - events[value] = name - - while monitor.poll(): - evt = recv_monitor_message(monitor) - if evt['event'] == 32: - self.connections += 1 - if evt['event'] == 512: - pass - - evt.update({'description': events[evt['event']]}) - - if evt['event'] == zmq.EVENT_MONITOR_STOPPED: - break - monitor.close() - print('event monitor thread done!') - - def initThrottle(self, throttle_value): - - @throttle(seconds=throttle_value) - def sendDataThrottled(fn_sender, data): - fn_sender(data) - - return sendDataThrottled - - def send(self, data): - log.log('Sending {}'.format(data), self.name) - self.sender.send_pyobj(data) - - def run(self): - - log.log( - 'Starting from {}'.format(self.address), - self.name - ) - - self.context = zmq.Context() - self.receiver = self.context.socket(zmq.PULL) - self.receiver.bind(self.address) - monitor = self.receiver.get_monitor_socket() - self.sender = self.context.socket(zmq.PUB) - if self.web_address: - log.success( - 'Sending to web: {}'.format(self.web_address), - self.name - ) - self.sender_web = self.context.socket(zmq.REQ) - self.sender_web.connect(self.web_address) - self.poll = zmq.Poller() - self.poll.register(self.sender_web, zmq.POLLIN) - time.sleep(1) - - if 'server' in self.kwargs: - self.sender.bind("ipc:///tmp/{}.plots".format(self.kwargs['server'])) - else: - self.sender.bind("ipc:///tmp/zmq.plots") - - time.sleep(2) - - t = Thread(target=self.event_monitor, args=(monitor,)) - t.start() - - while True: - dataOut = self.receiver.recv_pyobj() - if not dataOut.flagNoData: - if dataOut.type == 'Parameters': - tm = dataOut.utctimeInit - else: - tm = dataOut.utctime - if dataOut.useLocalTime: - if not self.localtime: - tm += time.timezone - dt = datetime.datetime.fromtimestamp(tm).date() - else: - if self.localtime: - tm -= time.timezone - dt = datetime.datetime.utcfromtimestamp(tm).date() - coerce = False - if dt not in self.dates: - if self.data: - self.data.ended = True - self.send(self.data) - coerce = True - self.data.setup() - self.dates.append(dt) - - self.data.update(dataOut, tm) - - if dataOut.finished is True: - self.connections -= 1 - if self.connections == 0 and dt in self.dates: - self.data.ended = True - self.send(self.data) - # self.data.setup() - time.sleep(1) - break - else: - if self.realtime: - self.send(self.data) - if self.web_address: - retries = 5 - while True: - self.sender_web.send(self.data.jsonify()) - socks = dict(self.poll.poll(5000)) - if socks.get(self.sender_web) == zmq.POLLIN: - reply = self.sender_web.recv_string() - if reply == 'ok': - log.log("Response from server ok", self.name) - break - else: - log.warning("Malformed reply from server: {}".format(reply), self.name) - - else: - log.warning("No response from server, retrying...", self.name) - self.sender_web.setsockopt(zmq.LINGER, 0) - self.sender_web.close() - self.poll.unregister(self.sender_web) - retries -= 1 - if retries == 0: - log.error("Server seems to be offline, abandoning", self.name) - self.sender_web = self.context.socket(zmq.REQ) - self.sender_web.connect(self.web_address) - self.poll.register(self.sender_web, zmq.POLLIN) - time.sleep(1) - break - self.sender_web = self.context.socket(zmq.REQ) - self.sender_web.connect(self.web_address) - self.poll.register(self.sender_web, zmq.POLLIN) - time.sleep(1) - else: - self.sendData(self.send, self.data, coerce=coerce) - coerce = False - - return - - class SendToFTP(Operation, Process): ''' diff --git a/schainpy/scripts/schain.xml b/schainpy/scripts/schain.xml deleted file mode 100644 index 80c9912..0000000 --- a/schainpy/scripts/schain.xml +++ /dev/null @@ -1 +0,0 @@ -