diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 41ce3ea..fcc4004 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -8,6 +8,7 @@ import re import datetime import copy import sys +import os import importlib import itertools from multiprocessing import Pool, TimeoutError @@ -18,12 +19,13 @@ from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papamete from .jroproc_base import ProcessingUnit, Operation, MPDecorator from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon from schainpy.model.data.jrodata import Spectra -#from scipy import asarray as ar,exp +from scipy import asarray as ar,exp from scipy.optimize import fmin, curve_fit from schainpy.utils import log import warnings from numpy import NaN from scipy.optimize.optimize import OptimizeWarning +from schainpy.model.io.jroIO_madrigal import MADWriter warnings.filterwarnings('ignore') import os @@ -31,6 +33,12 @@ import csv from scipy import signal import matplotlib.pyplot as plt +import json +import mysql.connector +from scipy.signal import savgol_filter +from scipy.signal import argrelextrema +from scipy.signal import find_peaks + SPEED_OF_LIGHT = 299792458 '''solving pickling issue''' @@ -51,7 +59,7 @@ def _unpickle_method(func_name, obj, cls): break return func.__get__(obj, cls) - +# @MPDecorator class ParametersProc(ProcessingUnit): METHODS = {} @@ -171,6 +179,8 @@ class ParametersProc(ProcessingUnit): self.dataOut.spc_noise = self.dataIn.getNoise() self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , 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 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels @@ -224,6 +234,7 @@ class ParametersProc(ProcessingUnit): self.__updateObjFromInput() self.dataOut.utctimeInit = self.dataIn.utctime self.dataOut.paramInterval = self.dataIn.timeInterval + return @@ -285,6 +296,7 @@ class RemoveWideGC(Operation): continue junk3 = numpy.squeeze(numpy.diff(j1index)) junk4 = numpy.squeeze(numpy.diff(j2index)) + valleyindex = j2index[numpy.where(junk4>1)] peakindex = j1index[numpy.where(junk3>1)] @@ -294,6 +306,7 @@ class RemoveWideGC(Operation): if numpy.size(isvalid) >1 : vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir]) isvalid = isvalid[vindex] + # clutter peak gcpeak = peakindex[isvalid] vl = numpy.where(valleyindex < gcpeak) @@ -2655,7 +2668,7 @@ class SpectralMoments(Operation): if (type1 is None): type1 = 0 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1 - if (snrth is None): snrth = -3 #-20.0 + if (snrth is None): snrth = -3 if (dc is None): dc = 0 if (aliasing is None): aliasing = 0 if (oldfd is None): oldfd = 0 @@ -3190,6 +3203,7 @@ class SpectralFitting(Operation): Operation.__init__(self) self.i=0 self.isConfig = False + self.aux = 1 def setup(self,nChan,nProf,nHei,nBlocks): self.__dataReady = False @@ -3212,6 +3226,7 @@ class SpectralFitting(Operation): if (wwauto is None): wwauto = 0 if (n0 < 1.e-20): n0 = 1.e-20 + freq = oldfreq vec_power = numpy.zeros(oldspec.shape[1]) vec_fd = numpy.zeros(oldspec.shape[1]) @@ -3620,7 +3635,7 @@ class SpectralFitting(Operation): junk = tmp tmp = junk*0.0 - tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))] + tmp[list(xpos + (ypos*num_prof))] = junk[list(xpos + (ypos*num_prof))] array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei)) return array @@ -3642,6 +3657,146 @@ class SpectralFitting(Operation): fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0]) smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp) return [fmom,numpy.sqrt(smom)] + + def find_max_excluding_neighbors_and_borders(self,matrix): + n, m = matrix.shape # Dimensions of the matrix + if n < 3 or m < 3: + raise ValueError("Matrix too small to exclude borders and neighbors of NaN.") + + # Create a mask for the borders + border_mask = numpy.zeros_like(matrix, dtype=bool) + border_mask[0, :] = True # Top border + border_mask[-1, :] = True # Bottom border + border_mask[:, 0] = True # Left border + border_mask[:, -1] = True # Right border + + # Exclude elements within 3 indices from the border + three_index_border_mask = numpy.zeros_like(matrix, dtype=bool) + three_index_border_mask[0:3, :] = True # Top 3 rows + three_index_border_mask[-3:, :] = True # Bottom 3 rows + three_index_border_mask[:, 0:3] = True # Left 3 columns + three_index_border_mask[:, -3:] = True # Right 3 columns + + # Create a mask for neighbors of NaNs + nan_mask = numpy.isnan(matrix) + nan_neighbors_mask = numpy.zeros_like(matrix, dtype=bool) + for i in range(1, n - 1): + for j in range(1, m - 1): + if nan_mask[i-1:i+2, j-1:j+2].any(): # Check 3x3 region around (i, j) + nan_neighbors_mask[i, j] = True + + # Combine all masks (borders, 3 indices from border, and NaN neighbors) + combined_mask = border_mask | three_index_border_mask | nan_neighbors_mask + + # Exclude these positions from the matrix + valid_matrix = numpy.where(combined_mask, numpy.nan, matrix) + + # Find the maximum value index in the valid matrix + if numpy.isnan(valid_matrix).all(): + return None # No valid maximum + max_index = numpy.unravel_index(numpy.nanargmax(valid_matrix), matrix.shape) + return max_index, valid_matrix + + def construct_three_phase_path(self,maximas_2, max_index): + """ + Constructs a path through maximas_2 starting from the middle point (max_index) + and going downwards, passing through the starting point, and then going upwards. + + Parameters: + - maximas_2 (numpy.ndarray): 2D array of points. + - max_index (tuple): Starting position as (row_index, starting_value). + + Returns: + - list: A path of points in the three phases: down, through the start, then up. + """ + row_index, starting_value = max_index + path = [] + # Initialize the path with the starting value + path_o = [starting_value] + path_up = [] + path_dw = [] + threshold = 4 + flag_broken = 0 + print(row_index,maximas_2.shape[0]) + # Phase 1: Going downwards (rows below the starting point) + current_value = starting_value + for i in range(row_index + 1, maximas_2.shape[0]): + + if flag_broken == 1: + path_up.append(numpy.nan) + continue + + distances = numpy.abs(maximas_2[i] - current_value) + closest_index = numpy.argmin(distances) + + if numpy.nanmin(distances) > threshold: + path_up.append(numpy.nan) + flag_broken = 1 + continue + + current_value = maximas_2[i, closest_index] + path_up.append(current_value) + + flag_broken = 0 + + # Phase 2: Going upwards (rows above the starting point) + current_value = starting_value + for i in range(row_index - 1, -1, -1): + if flag_broken == 1: + path_dw.append(numpy.nan) + continue + distances = numpy.abs(maximas_2[i] - current_value) + closest_index = numpy.argmin(distances) + + if numpy.nanmin(distances) > threshold: + path_dw.append(numpy.nan) + flag_broken = 1 + continue + + current_value = maximas_2[i, closest_index] + path_dw.append(current_value) + path.append(path_dw[::-1]) + path.append(path_o) + path.append(path_up) + return numpy.concatenate(path) + + def Vr_correction(self,dataset): + + dataset = savgol_filter(dataset, window_length=15, polyorder=3, axis=1) + dataset = savgol_filter(dataset, window_length=15, polyorder=3, axis=1) + dataset = savgol_filter(dataset, window_length=7, polyorder=3, axis=0) + dataset[:12,:] = numpy.nan + dataset[45:, :] = numpy.nan + max_index, _ = self.find_max_excluding_neighbors_and_borders(dataset) + max_index = numpy.array(max_index) + print(max_index) + maximas = argrelextrema(dataset, numpy.greater, axis=1) + heighs = numpy.unique(maximas[0]) + #grouped_arrays = {value: maximas[0][numpy.where(maximas[0] == value)] for value in numpy.unique(maximas[0])} + grouped_arrays = [maximas[0][numpy.where(maximas[0] == value)] for value in numpy.unique(maximas[0])] + + maximas_2_list = [] + num_maximas = 3 + for i, h in enumerate(heighs): + n = len(grouped_arrays[i]) + idx = numpy.where(maximas[0] == h)[0] + maxs_ = numpy.argsort(dataset[h, maximas[1][idx]])[::-1][:num_maximas] + maxs = maximas[1][idx][maxs_] + heigh = numpy.ones(len(maxs)) * h + arr = numpy.array([heigh, maxs]) + if n < num_maximas: + arr = numpy.hstack((arr, numpy.full((arr.shape[0], num_maximas - arr.shape[1]), numpy.nan))) + maximas_2_list.append(arr) + + maximas_2 = numpy.array(maximas_2_list) + _maximas_2_ = maximas_2[:, 1] + # + max_index[0] = numpy.where(heighs == max_index[0])[0][0] + # + path = self.construct_three_phase_path(_maximas_2_, max_index) + return int(numpy.nanmedian(path)) + + def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints): ''' @@ -4010,7 +4165,6 @@ class SpectralFitting(Operation): power = numpy.sum(spectra, axis=1) nPairs = len(crosspairs) absc = dataOut.abscissaList[:-1] - print('para escribir h5 ',dataOut.paramInterval) if not self.isConfig: self.isConfig = True @@ -4024,6 +4178,7 @@ class SpectralFitting(Operation): my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th) clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index) + dataOut.data_spc = incoh_spectra dataOut.data_cspc = incoh_cspectra dataOut.sat_spectra = sat_spectra @@ -4138,7 +4293,7 @@ class SpectralFitting(Operation): dataOut.nIncohInt= int(dataOut.nIncohInt) dataOut.nProfiles = int(dataOut.nProfiles) dataOut.nCohInt = int(dataOut.nCohInt) - print('sale',dataOut.nProfiles,dataOut.nHeights) + #print('sale',dataOut.nProfiles,dataOut.nHeights) #dataOut.nFFTPoints=nprofs #dataOut.normFactor = nprofs dataOut.channelList = channelList @@ -4148,7 +4303,7 @@ class SpectralFitting(Operation): vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds) #dataOut.ippSeconds=ipp absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf - print('sale 2',dataOut.ippSeconds,M,N) + #print('sale 2',dataOut.ippSeconds,M,N) print('Empieza procesamiento offline') if path != None: sys.path.append(path) @@ -4170,6 +4325,8 @@ class SpectralFitting(Operation): dataOut.data_snr1_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan # dataOut.smooth_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan + if self.aux == 1: + dataOut.p03last = numpy.zeros(nGroups, 'float32') for i in range(nGroups): coord = groupArray[i,:] #Input data array @@ -4193,7 +4350,7 @@ class SpectralFitting(Operation): n1= n1i my_noises[2*i+0] = n0 my_noises[2*i+1] = n1 - snrth = -13 #-13.0 # -4 -16 -25 + snrth = -12.0 #-11.0 # -4 -16 -25 snrth = 10**(snrth/10.0) jvelr = numpy.zeros(nHeights, dtype = 'float') #snr0 = numpy.zeros(nHeights, dtype = 'float') @@ -4201,7 +4358,9 @@ class SpectralFitting(Operation): hvalid = [0] coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:]) - + + SIGNAL_0 = [] + SIGNAL_1 = [] for h in range(nHeights): smooth = clean_num_aver[i+1,h] signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth @@ -4228,10 +4387,32 @@ class SpectralFitting(Operation): else: jvelr[h] = absc[0] if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None) #print(maxp0,absc[maxp0],snr0,jvelr[h]) + SIGNAL_0.append(signal0) + SIGNAL_1.append(signal1) if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1 - else: fd0 = numpy.nan + else: fd0 = numpy.nan print(fd0) + fd0n = [fd0] + SIGNAL_0 = numpy.array(SIGNAL_0) + SIGNAL_1 = numpy.array(SIGNAL_1) + ### ---- Signal correction + + diff = numpy.abs(dataOut.p03last[i] - fd0) + if (filec != None) and (self.aux == 0) and (diff > 6): + print('hace correccion de vr') + try: + maxs_0 = self.Vr_correction(SIGNAL_0) + maxs_1 = self.Vr_correction(SIGNAL_1) + fd0_aux = -1*(absc[maxs_0] + absc[maxs_1]) / 2 + if numpy.abs(dataOut.p03last[i] - fd0_aux ) > diff: print("Wrong correction omitted", fd0_aux) + else: fd0 = fd0_aux; print("Changed fd0: ", fd0) + except Exception as e: + print("Error in correction processes skiped: ", e) + + dataOut.p03last[i] = fd0 + ### --- + fd0n = [fd0] for h in range(nHeights): d = data[:,h] smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:] @@ -4294,6 +4475,16 @@ class SpectralFitting(Operation): minp = p0*numpy.nan error0 = numpy.nan error1 = p0*numpy.nan +# s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0)) +# covp = covp*s_sq +# error = [] +# for ip in range(len(minp)): +# try: +# error.append(numpy.absolute(covp[ip][ip])**0.5) +# except: +# error.append( 0.00 ) + #if i==1 and h==11 and index == 139: print(p0, minp,data_spc) + fd0n = numpy.concatenate((fd0n,minp[3]), axis=None) else : data_spc = dataOut.data_spc[coord,:,h] p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants)) @@ -4308,7 +4499,25 @@ class SpectralFitting(Operation): dataOut.data_param[i,:,h] = minp dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0 dataOut.data_snr1_i[i*2+1,h] = numpy.sum(signalpn1/(nProf-1))/n1 - + #dataOut.smooth_i[i*2,h] = clean_num_aver[i+1,h] + #print(fd0,dataOut.data_param[i,3,h]) + #print(fd0,dataOut.data_param[i,3,:]) + fd0n = fd0n[1:] + dvr = abs(dataOut.data_param[i,3,:]-numpy.nanmedian(fd0n)) + indbad=numpy.where(dvr > 25) + esf = 0 + if filec != None: + esf = self.weightf.esffit(esf,tini.tm_year,tini.tm_yday,index) + #print(indbad,'media ',numpy.nanmedian(fd0n), esf) + + #plt.plot(abs(dataOut.data_param[i,3,:]-numpy.nanmedian(fd0n)),'o') + #plt.plot(dataOut.data_param[i,3,:]) + #plt.ylim(-200,200) + if esf == 0 : + dataOut.data_param[i,3,indbad] = numpy.nanmedian(fd0n) + print('corrige') + #plt.plot(dataOut.data_param[i,3,:]) + #plt.show() for ht in range(nHeights-1) : smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam] dataOut.data_paramC[4*i,ht,1] = smooth @@ -4393,6 +4602,7 @@ class SpectralFitting(Operation): sat_fits[1+2*beam,ht] = numpy.sum(signalpn1/(val1_npoints*nProf))/n1 dataOut.sat_fits = sat_fits + self.aux = 0 return dataOut def __residFunction(self, p, dp, LT, constants): @@ -4588,7 +4798,6 @@ class WindProfiler(Operation): return distx, disty, dist, ang #Calculo de Matrices - def __calculateVelVer(self, phase, lagTRange, _lambda): Ts = lagTRange[1] - lagTRange[0] @@ -4715,8 +4924,6 @@ class WindProfiler(Operation): azim = meteorAux[:, 3]*numpy.pi/180 n = numpy.cos(zen) - # m = (1 - n**2)/(1 - numpy.tan(azim)**2) - # l = m*numpy.tan(azim) l = numpy.sin(zen)*numpy.sin(azim) m = numpy.sin(zen)*numpy.cos(azim) @@ -4781,7 +4988,6 @@ class WindProfiler(Operation): for t in uniqueTime: metArray1 = metArray[utctime==t,:] -# phaseDerThresh = numpy.pi/4 #reducir Phase thresh tmet = metArray1[:,1].astype(int) hmet = metArray1[:,2].astype(int) @@ -4912,8 +5118,7 @@ class WindProfiler(Operation): def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs): - param = dataOut.data_param - #if dataOut.abscissaList != None: + param = dataOut.moments ### LA VERSION ANTERIOR trabajaba con param = dataOut.data_param if numpy.any(dataOut.abscissaList): absc = dataOut.abscissaList[:-1] # noise = dataOut.noise @@ -4933,29 +5138,9 @@ class WindProfiler(Operation): elif technique == 'SA': #Parameters -# position_x = kwargs['positionX'] -# position_y = kwargs['positionY'] -# azimuth = kwargs['azimuth'] -# -# if kwargs.has_key('crosspairsList'): -# pairs = kwargs['crosspairsList'] -# else: -# pairs = None -# -# if kwargs.has_key('correctFactor'): -# correctFactor = kwargs['correctFactor'] -# else: -# correctFactor = 1 - -# tau = dataOut.data_param -# _lambda = dataOut.C/dataOut.frequency -# pairsList = dataOut.groupList -# nChannels = dataOut.nChannels - kwargs['groupList'] = dataOut.groupList kwargs['tau'] = dataOut.data_param kwargs['_lambda'] = dataOut.C/dataOut.frequency -# dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) dataOut.data_output = self.techniqueSA(kwargs) dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval @@ -4984,7 +5169,6 @@ class WindProfiler(Operation): dataOut.outputInterval = nHours*3600 if self.__isConfig == False: -# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() @@ -5046,7 +5230,6 @@ class WindProfiler(Operation): if self.__isConfig == False: dataOut.outputInterval = nMins*60 -# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) minuteAux = initime.minute @@ -5077,69 +5260,19 @@ class WindProfiler(Operation): dataOut.flagNoData = False self.__buffer = None - return - -class WindProfiler(Operation): - - __isConfig = False - - __initime = None - __lastdatatime = None - __integrationtime = None - - __buffer = None - - __dataReady = False - - __firstdata = None + return dataOut - n = None +class EWDriftsEstimation(Operation): def __init__(self): Operation.__init__(self) - def __calculateCosDir(self, elev, azim): - zen = (90 - elev)*numpy.pi/180 - azim = azim*numpy.pi/180 - cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) - cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2) - - signX = numpy.sign(numpy.cos(azim)) - signY = numpy.sign(numpy.sin(azim)) - - cosDirX = numpy.copysign(cosDirX, signX) - cosDirY = numpy.copysign(cosDirY, signY) - return cosDirX, cosDirY - - def __calculateAngles(self, theta_x, theta_y, azimuth): - - dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2) - zenith_arr = numpy.arccos(dir_cosw) - azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180 - - dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr) - dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr) - - return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw - - def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly): - - if horOnly: - A = numpy.c_[dir_cosu,dir_cosv] - else: - A = numpy.c_[dir_cosu,dir_cosv,dir_cosw] - A = numpy.asmatrix(A) - A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose() - - return A1 - def __correctValues(self, heiRang, phi, velRadial, SNR): listPhi = phi.tolist() maxid = listPhi.index(max(listPhi)) minid = listPhi.index(min(listPhi)) rango = list(range(len(phi))) - heiRang1 = heiRang*math.cos(phi[maxid]) heiRangAux = heiRang*math.cos(phi[minid]) indOut = (heiRang1 < heiRangAux[0]).nonzero() @@ -5151,13 +5284,18 @@ class WindProfiler(Operation): for i in rango: x = heiRang*math.cos(phi[i]) y1 = velRadial[i,:] - f1 = interpolate.interp1d(x,y1,kind = 'cubic') - + vali= (numpy.isfinite(y1)==True).nonzero() + y1=y1[vali] + x = x[vali] + f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False) x1 = heiRang1 y11 = f1(x1) - y2 = SNR[i,:] - f2 = interpolate.interp1d(x,y2,kind = 'cubic') + x = heiRang*math.cos(phi[i]) + vali= (y2 != -1).nonzero() + y2 = y2[vali] + x = x[vali] + f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False) y21 = f2(x1) velRadial1[i,:] = y11 @@ -5165,713 +5303,124 @@ class WindProfiler(Operation): return heiRang1, velRadial1, SNR1 - def __calculateVelUVW(self, A, velRadial): - - #Operacion Matricial - velUVW = numpy.zeros((A.shape[0],velRadial.shape[1])) - velUVW[:,:] = numpy.dot(A,velRadial) - - - return velUVW + def run(self, dataOut, zenith, zenithCorrection,fileDrifts,beam_pos=None): + dataOut.lat = -11.95 + dataOut.lon = -76.87 + dataOut.spcst = 0.00666 + dataOut.pl = 0.0003 + dataOut.cbadn = 3 + dataOut.inttms = 300 + if numpy.any(beam_pos) : + dataOut.azw = beam_pos[0] # -115.687 + dataOut.elw = beam_pos[1] # 86.1095 + dataOut.aze = beam_pos[2] # 130.052 + dataOut.ele = beam_pos[3] # 87.6558 + dataOut.jro14 = numpy.log10(dataOut.spc_noise[0]/dataOut.normFactor) + dataOut.jro15 = numpy.log10(dataOut.spc_noise[1]/dataOut.normFactor) + dataOut.jro16 = numpy.log10(dataOut.spc_noise[2]/dataOut.normFactor) + dataOut.nwlos = numpy.log10(dataOut.spc_noise[3]/dataOut.normFactor) - def techniqueDBS(self, kwargs): - """ - Function that implements Doppler Beam Swinging (DBS) technique. + heiRang = dataOut.heightList + velRadial = dataOut.data_param[:,3,:] + velRadialm = dataOut.data_param[:,2:4,:]*-1 - Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, - Direction correction (if necessary), Ranges and SNR + rbufc=dataOut.data_paramC[:,:,0] + ebufc=dataOut.data_paramC[:,:,1] + SNR = dataOut.data_snr1_i + rbufi = dataOut.data_snr1_i + velRerr = dataOut.data_error[:,4,:] + range1 = dataOut.heightList + nhei = len(range1) + + sat_fits = dataOut.sat_fits + + channels = dataOut.channelList + nChan = len(channels) + my_nbeams = nChan/2 + if my_nbeams == 2: + moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]])) + else : + moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]])) + dataOut.moments=moments + #Incoherent + smooth_w = dataOut.clean_num_aver[0,:] + chisq_w = dataOut.data_error[0,0,:] + p_w0 = rbufi[0,:] + p_w1 = rbufi[1,:] - Output: Winds estimation (Zonal, Meridional and Vertical) + # Coherent + smooth_wC = ebufc[0,:] + p_w0C = rbufc[0,:] + p_w1C = rbufc[1,:] + w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1) + t_wC = rbufc[3,:] + val = (numpy.isfinite(p_w0)==False).nonzero() + p_w0[val]=0 + val = (numpy.isfinite(p_w1)==False).nonzero() + p_w1[val]=0 + val = (numpy.isfinite(p_w0C)==False).nonzero() + p_w0C[val]=0 + val = (numpy.isfinite(p_w1C)==False).nonzero() + p_w1C[val]=0 + val = (numpy.isfinite(smooth_w)==False).nonzero() + smooth_w[val]=0 + val = (numpy.isfinite(smooth_wC)==False).nonzero() + smooth_wC[val]=0 + + p_w0_a = (p_w0*smooth_w+p_w0C*smooth_wC)/(smooth_w+smooth_wC) + p_w1_a = (p_w1*smooth_w+p_w1C*smooth_wC)/(smooth_w+smooth_wC) - Parameters affected: Winds, height range, SNR - """ - velRadial0 = kwargs['velRadial'] - heiRang = kwargs['heightList'] - SNR0 = kwargs['SNR'] + if len(sat_fits) >0 : + p_w0C = p_w0C + sat_fits[0,:] + p_w1C = p_w1C + sat_fits[1,:] + + if my_nbeams == 1: + w = velRadial[0,:] + winds = velRadial.copy() + w_err = velRerr[0,:] + werrtmp = numpy.where(numpy.isfinite(w)==True,w_err,numpy.nan) + w_err = werrtmp.copy() + u = w*numpy.nan + u_err = w_err*numpy.nan + p_e0 = p_w0*numpy.nan + p_e1 = p_w1*numpy.nan + #snr1 = 10*numpy.log10(SNR[0]) + if my_nbeams == 2: - if 'dirCosx' in kwargs and 'dirCosy' in kwargs: - theta_x = numpy.array(kwargs['dirCosx']) - theta_y = numpy.array(kwargs['dirCosy']) - else: - elev = numpy.array(kwargs['elevation']) - azim = numpy.array(kwargs['azimuth']) - theta_x, theta_y = self.__calculateCosDir(elev, azim) - azimuth = kwargs['correctAzimuth'] - if 'horizontalOnly' in kwargs: - horizontalOnly = kwargs['horizontalOnly'] - else: horizontalOnly = False - if 'correctFactor' in kwargs: - correctFactor = kwargs['correctFactor'] - else: correctFactor = 1 - if 'channelList' in kwargs: - channelList = kwargs['channelList'] - if len(channelList) == 2: - horizontalOnly = True - arrayChannel = numpy.array(channelList) - param = param[arrayChannel,:,:] - theta_x = theta_x[arrayChannel] - theta_y = theta_y[arrayChannel] + zenith = numpy.array(zenith) + zenith -= zenithCorrection + zenith *= numpy.pi/180 + if zenithCorrection != 0 : + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) + else : + heiRang1 = heiRang + velRadial1 = velRadial + SNR1 = SNR + + alp = zenith[0] + bet = zenith[1] - azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) - heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0) - A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly) + w_w = velRadial1[0,:] + w_e = velRadial1[1,:] + w_w_err = velRerr[0,:] + w_e_err = velRerr[1,:] + smooth_e = dataOut.clean_num_aver[2,:] + chisq_e = dataOut.data_error[1,0,:] + p_e0 = rbufi[2,:] + p_e1 = rbufi[3,:] - #Calculo de Componentes de la velocidad con DBS - winds = self.__calculateVelUVW(A,velRadial1) + tini=time.localtime(dataOut.utctime) - return winds, heiRang1, SNR1 - - def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None): - - nPairs = len(pairs_ccf) - posx = numpy.asarray(posx) - posy = numpy.asarray(posy) - - #Rotacion Inversa para alinear con el azimuth - if azimuth!= None: - azimuth = azimuth*math.pi/180 - posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth) - posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth) - else: - posx1 = posx - posy1 = posy - - #Calculo de Distancias - distx = numpy.zeros(nPairs) - disty = numpy.zeros(nPairs) - dist = numpy.zeros(nPairs) - ang = numpy.zeros(nPairs) - - for i in range(nPairs): - distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]] - disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]] - dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) - ang[i] = numpy.arctan2(disty[i],distx[i]) - - return distx, disty, dist, ang - #Calculo de Matrices - - def __calculateVelVer(self, phase, lagTRange, _lambda): - - Ts = lagTRange[1] - lagTRange[0] - velW = -_lambda*phase/(4*math.pi*Ts) - - return velW - - def __calculateVelHorDir(self, dist, tau1, tau2, ang): - nPairs = tau1.shape[0] - nHeights = tau1.shape[1] - vel = numpy.zeros((nPairs,3,nHeights)) - dist1 = numpy.reshape(dist, (dist.size,1)) - - angCos = numpy.cos(ang) - angSin = numpy.sin(ang) - - vel0 = dist1*tau1/(2*tau2**2) - vel[:,0,:] = (vel0*angCos).sum(axis = 1) - vel[:,1,:] = (vel0*angSin).sum(axis = 1) - - ind = numpy.where(numpy.isinf(vel)) - vel[ind] = numpy.nan - - return vel - - def techniqueSA(self, kwargs): - - """ - Function that implements Spaced Antenna (SA) technique. - - Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, - Direction correction (if necessary), Ranges and SNR - - Output: Winds estimation (Zonal, Meridional and Vertical) - - Parameters affected: Winds - """ - position_x = kwargs['positionX'] - position_y = kwargs['positionY'] - azimuth = kwargs['azimuth'] - - if 'correctFactor' in kwargs: - correctFactor = kwargs['correctFactor'] - else: - correctFactor = 1 - - groupList = kwargs['groupList'] - pairs_ccf = groupList[1] - tau = kwargs['tau'] - _lambda = kwargs['_lambda'] - - #Cross Correlation pairs obtained - - indtau = tau.shape[0]/2 - tau1 = tau[:indtau,:] - tau2 = tau[indtau:-1,:] - phase1 = tau[-1,:] - - #--------------------------------------------------------------------- - #Metodo Directo - distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth) - winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) - winds = stats.nanmean(winds, axis=0) - #--------------------------------------------------------------------- - #Metodo General - - #--------------------------------------------------------------------- - winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) - winds = correctFactor*winds - return winds - - def __checkTime(self, currentTime, paramInterval, outputInterval): - - dataTime = currentTime + paramInterval - deltaTime = dataTime - self.__initime - - if deltaTime >= outputInterval or deltaTime < 0: - self.__dataReady = True - return - - def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax): - ''' - Function that implements winds estimation technique with detected meteors. - - Input: Detected meteors, Minimum meteor quantity to wind estimation - - Output: Winds estimation (Zonal and Meridional) - - Parameters affected: Winds - ''' - #Settings - nInt = (heightMax - heightMin)/2 - nInt = int(nInt) - winds = numpy.zeros((2,nInt))*numpy.nan - - #Filter errors - error = numpy.where(arrayMeteor[:,-1] == 0)[0] - finalMeteor = arrayMeteor[error,:] - - #Meteor Histogram - finalHeights = finalMeteor[:,2] - hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax)) - nMeteorsPerI = hist[0] - heightPerI = hist[1] - - #Sort of meteors - indSort = finalHeights.argsort() - finalMeteor2 = finalMeteor[indSort,:] - - # Calculating winds - ind1 = 0 - ind2 = 0 - - for i in range(nInt): - nMet = nMeteorsPerI[i] - ind1 = ind2 - ind2 = ind1 + nMet - - meteorAux = finalMeteor2[ind1:ind2,:] - - if meteorAux.shape[0] >= meteorThresh: - vel = meteorAux[:, 6] - zen = meteorAux[:, 4]*numpy.pi/180 - azim = meteorAux[:, 3]*numpy.pi/180 - - n = numpy.cos(zen) - l = numpy.sin(zen)*numpy.sin(azim) - m = numpy.sin(zen)*numpy.cos(azim) - - A = numpy.vstack((l, m)).transpose() - A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose()) - windsAux = numpy.dot(A1, vel) - - winds[0,i] = windsAux[0] - winds[1,i] = windsAux[1] - - return winds, heightPerI[:-1] - - def techniqueNSM_SA(self, **kwargs): - metArray = kwargs['metArray'] - heightList = kwargs['heightList'] - timeList = kwargs['timeList'] - - rx_location = kwargs['rx_location'] - groupList = kwargs['groupList'] - azimuth = kwargs['azimuth'] - dfactor = kwargs['dfactor'] - k = kwargs['k'] - - azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth) - d = dist*dfactor - #Phase calculation - metArray1 = self.__getPhaseSlope(metArray, heightList, timeList) - - metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities - - velEst = numpy.zeros((heightList.size,2))*numpy.nan - azimuth1 = azimuth1*numpy.pi/180 - - for i in range(heightList.size): - h = heightList[i] - indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0] - metHeight = metArray1[indH,:] - if metHeight.shape[0] >= 2: - velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities - iazim = metHeight[:,1].astype(int) - azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths - A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux))) - A = numpy.asmatrix(A) - A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose() - velHor = numpy.dot(A1,velAux) - - velEst[i,:] = numpy.squeeze(velHor) - return velEst - - def __getPhaseSlope(self, metArray, heightList, timeList): - meteorList = [] - #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2 - #Putting back together the meteor matrix - utctime = metArray[:,0] - uniqueTime = numpy.unique(utctime) - - phaseDerThresh = 0.5 - ippSeconds = timeList[1] - timeList[0] - sec = numpy.where(timeList>1)[0][0] - nPairs = metArray.shape[1] - 6 - nHeights = len(heightList) - - for t in uniqueTime: - metArray1 = metArray[utctime==t,:] - tmet = metArray1[:,1].astype(int) - hmet = metArray1[:,2].astype(int) - - metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1)) - metPhase[:,:] = numpy.nan - metPhase[:,hmet,tmet] = metArray1[:,6:].T - - #Delete short trails - metBool = ~numpy.isnan(metPhase[0,:,:]) - heightVect = numpy.sum(metBool, axis = 1) - metBool[heightVect phaseDerThresh)) - metPhase[phDerAux] = numpy.nan - - #--------------------------METEOR DETECTION ----------------------------------------- - indMet = numpy.where(numpy.any(metBool,axis=1))[0] - - for p in numpy.arange(nPairs): - phase = metPhase[p,:,:] - phDer = metDer[p,:,:] - - for h in indMet: - height = heightList[h] - phase1 = phase[h,:] #82 - phDer1 = phDer[h,:] - - phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap - - indValid = numpy.where(~numpy.isnan(phase1))[0] - initMet = indValid[0] - endMet = 0 - - for i in range(len(indValid)-1): - - #Time difference - inow = indValid[i] - inext = indValid[i+1] - idiff = inext - inow - #Phase difference - phDiff = numpy.abs(phase1[inext] - phase1[inow]) - - if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor - sizeTrail = inow - initMet + 1 - if sizeTrail>3*sec: #Too short meteors - x = numpy.arange(initMet,inow+1)*ippSeconds - y = phase1[initMet:inow+1] - ynnan = ~numpy.isnan(y) - x = x[ynnan] - y = y[ynnan] - slope, intercept, r_value, p_value, std_err = stats.linregress(x,y) - ylin = x*slope + intercept - rsq = r_value**2 - if rsq > 0.5: - vel = slope#*height*1000/(k*d) - estAux = numpy.array([utctime,p,height, vel, rsq]) - meteorList.append(estAux) - initMet = inext - metArray2 = numpy.array(meteorList) - - return metArray2 - - def __calculateAzimuth1(self, rx_location, pairslist, azimuth0): - - azimuth1 = numpy.zeros(len(pairslist)) - dist = numpy.zeros(len(pairslist)) - - for i in range(len(rx_location)): - ch0 = pairslist[i][0] - ch1 = pairslist[i][1] - - diffX = rx_location[ch0][0] - rx_location[ch1][0] - diffY = rx_location[ch0][1] - rx_location[ch1][1] - azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi - dist[i] = numpy.sqrt(diffX**2 + diffY**2) - - azimuth1 -= azimuth0 - return azimuth1, dist - - def techniqueNSM_DBS(self, **kwargs): - metArray = kwargs['metArray'] - heightList = kwargs['heightList'] - timeList = kwargs['timeList'] - azimuth = kwargs['azimuth'] - theta_x = numpy.array(kwargs['theta_x']) - theta_y = numpy.array(kwargs['theta_y']) - - utctime = metArray[:,0] - cmet = metArray[:,1].astype(int) - hmet = metArray[:,3].astype(int) - SNRmet = metArray[:,4] - vmet = metArray[:,5] - spcmet = metArray[:,6] - - nChan = numpy.max(cmet) + 1 - nHeights = len(heightList) - - azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) - hmet = heightList[hmet] - h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights - - velEst = numpy.zeros((heightList.size,2))*numpy.nan - - for i in range(nHeights - 1): - hmin = heightList[i] - hmax = heightList[i + 1] - - thisH = (h1met>=hmin) & (h1met8) & (vmet<50) & (spcmet<10) - indthisH = numpy.where(thisH) - - if numpy.size(indthisH) > 3: - - vel_aux = vmet[thisH] - chan_aux = cmet[thisH] - cosu_aux = dir_cosu[chan_aux] - cosv_aux = dir_cosv[chan_aux] - cosw_aux = dir_cosw[chan_aux] - - nch = numpy.size(numpy.unique(chan_aux)) - if nch > 1: - A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True) - velEst[i,:] = numpy.dot(A,vel_aux) - - return velEst - - def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs): - - param = dataOut.moments - if numpy.any(dataOut.abscissaList): - absc = dataOut.abscissaList[:-1] - # noise = dataOut.noise - heightList = dataOut.heightList - SNR = dataOut.data_snr - - if technique == 'DBS': - - kwargs['velRadial'] = param[:,1,:] #Radial velocity - kwargs['heightList'] = heightList - kwargs['SNR'] = SNR - - dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function - dataOut.utctimeInit = dataOut.utctime - dataOut.outputInterval = dataOut.paramInterval - - elif technique == 'SA': - - #Parameters - kwargs['groupList'] = dataOut.groupList - kwargs['tau'] = dataOut.data_param - kwargs['_lambda'] = dataOut.C/dataOut.frequency - dataOut.data_output = self.techniqueSA(kwargs) - dataOut.utctimeInit = dataOut.utctime - dataOut.outputInterval = dataOut.timeInterval - - elif technique == 'Meteors': - dataOut.flagNoData = True - self.__dataReady = False - - if 'nHours' in kwargs: - nHours = kwargs['nHours'] - else: - nHours = 1 - - if 'meteorsPerBin' in kwargs: - meteorThresh = kwargs['meteorsPerBin'] - else: - meteorThresh = 6 - - if 'hmin' in kwargs: - hmin = kwargs['hmin'] - else: hmin = 70 - if 'hmax' in kwargs: - hmax = kwargs['hmax'] - else: hmax = 110 - - dataOut.outputInterval = nHours*3600 - - if self.__isConfig == False: - #Get Initial LTC time - self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) - self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() - - self.__isConfig = True - - if self.__buffer is None: - self.__buffer = dataOut.data_param - self.__firstdata = copy.copy(dataOut) - - else: - self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - - self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready - - if self.__dataReady: - dataOut.utctimeInit = self.__initime - - self.__initime += dataOut.outputInterval #to erase time offset - - dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) - dataOut.flagNoData = False - self.__buffer = None - - elif technique == 'Meteors1': - dataOut.flagNoData = True - self.__dataReady = False - - if 'nMins' in kwargs: - nMins = kwargs['nMins'] - else: nMins = 20 - if 'rx_location' in kwargs: - rx_location = kwargs['rx_location'] - else: rx_location = [(0,1),(1,1),(1,0)] - if 'azimuth' in kwargs: - azimuth = kwargs['azimuth'] - else: azimuth = 51.06 - if 'dfactor' in kwargs: - dfactor = kwargs['dfactor'] - if 'mode' in kwargs: - mode = kwargs['mode'] - if 'theta_x' in kwargs: - theta_x = kwargs['theta_x'] - if 'theta_y' in kwargs: - theta_y = kwargs['theta_y'] - else: mode = 'SA' - - #Borrar luego esto - if dataOut.groupList is None: - dataOut.groupList = [(0,1),(0,2),(1,2)] - groupList = dataOut.groupList - C = 3e8 - freq = 50e6 - lamb = C/freq - k = 2*numpy.pi/lamb - - timeList = dataOut.abscissaList - heightList = dataOut.heightList - - if self.__isConfig == False: - dataOut.outputInterval = nMins*60 - #Get Initial LTC time - initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) - minuteAux = initime.minute - minuteNew = int(numpy.floor(minuteAux/nMins)*nMins) - self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() - - self.__isConfig = True - - if self.__buffer is None: - self.__buffer = dataOut.data_param - self.__firstdata = copy.copy(dataOut) - - else: - self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - - self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready - - if self.__dataReady: - dataOut.utctimeInit = self.__initime - self.__initime += dataOut.outputInterval #to erase time offset - - metArray = self.__buffer - if mode == 'SA': - dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList) - elif mode == 'DBS': - dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y) - dataOut.data_output = dataOut.data_output.T - dataOut.flagNoData = False - self.__buffer = None - - return dataOut - -class EWDriftsEstimation(Operation): - - def __init__(self): - Operation.__init__(self) - - def __correctValues(self, heiRang, phi, velRadial, SNR): - listPhi = phi.tolist() - maxid = listPhi.index(max(listPhi)) - minid = listPhi.index(min(listPhi)) - - rango = list(range(len(phi))) - heiRang1 = heiRang*math.cos(phi[maxid]) - heiRangAux = heiRang*math.cos(phi[minid]) - indOut = (heiRang1 < heiRangAux[0]).nonzero() - heiRang1 = numpy.delete(heiRang1,indOut) - - velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) - SNR1 = numpy.zeros([len(phi),len(heiRang1)]) - - for i in rango: - x = heiRang*math.cos(phi[i]) - y1 = velRadial[i,:] - vali= (numpy.isfinite(y1)==True).nonzero() - y1=y1[vali] - x = x[vali] - f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False) - x1 = heiRang1 - y11 = f1(x1) - y2 = SNR[i,:] - x = heiRang*math.cos(phi[i]) - vali= (y2 != -1).nonzero() - y2 = y2[vali] - x = x[vali] - f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False) - y21 = f2(x1) - - velRadial1[i,:] = y11 - SNR1[i,:] = y21 - - return heiRang1, velRadial1, SNR1 - - def run(self, dataOut, zenith, zenithCorrection,fileDrifts): - - dataOut.lat = -11.95 - dataOut.lon = -76.87 - dataOut.spcst = 0.00666 - dataOut.pl = 0.0003 - dataOut.cbadn = 3 - dataOut.inttms = 300 - dataOut.azw = -115.687 - dataOut.elw = 86.1095 - dataOut.aze = 130.052 - dataOut.ele = 87.6558 - dataOut.jro14 = numpy.log10(dataOut.spc_noise[0]/dataOut.normFactor) - dataOut.jro15 = numpy.log10(dataOut.spc_noise[1]/dataOut.normFactor) - dataOut.jro16 = numpy.log10(dataOut.spc_noise[2]/dataOut.normFactor) - dataOut.nwlos = numpy.log10(dataOut.spc_noise[3]/dataOut.normFactor) - - heiRang = dataOut.heightList - velRadial = dataOut.data_param[:,3,:] - velRadialm = dataOut.data_param[:,2:4,:]*-1 - - rbufc=dataOut.data_paramC[:,:,0] - ebufc=dataOut.data_paramC[:,:,1] - SNR = dataOut.data_snr1_i - rbufi = dataOut.data_snr1_i - velRerr = dataOut.data_error[:,4,:] - range1 = dataOut.heightList - nhei = len(range1) - - sat_fits = dataOut.sat_fits - - channels = dataOut.channelList - nChan = len(channels) - my_nbeams = nChan/2 - if my_nbeams == 2: - moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]])) - else : - moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]])) - dataOut.moments=moments - #Incoherent - smooth_w = dataOut.clean_num_aver[0,:] - chisq_w = dataOut.data_error[0,0,:] - p_w0 = rbufi[0,:] - p_w1 = rbufi[1,:] - - # Coherent - smooth_wC = ebufc[0,:] - p_w0C = rbufc[0,:] - p_w1C = rbufc[1,:] - w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1) - t_wC = rbufc[3,:] - val = (numpy.isfinite(p_w0)==False).nonzero() - p_w0[val]=0 - val = (numpy.isfinite(p_w1)==False).nonzero() - p_w1[val]=0 - val = (numpy.isfinite(p_w0C)==False).nonzero() - p_w0C[val]=0 - val = (numpy.isfinite(p_w1C)==False).nonzero() - p_w1C[val]=0 - val = (numpy.isfinite(smooth_w)==False).nonzero() - smooth_w[val]=0 - val = (numpy.isfinite(smooth_wC)==False).nonzero() - smooth_wC[val]=0 - - #p_w0 = (p_w0*smooth_w+p_w0C*smooth_wC)/(smooth_w+smooth_wC) - #p_w1 = (p_w1*smooth_w+p_w1C*smooth_wC)/(smooth_w+smooth_wC) - - if len(sat_fits) >0 : - p_w0C = p_w0C + sat_fits[0,:] - p_w1C = p_w1C + sat_fits[1,:] - - if my_nbeams == 1: - w = velRadial[0,:] - winds = velRadial.copy() - w_err = velRerr[0,:] - u = w*numpy.nan - u_err = w_err*numpy.nan - p_e0 = p_w0*numpy.nan - p_e1 = p_w1*numpy.nan - #snr1 = 10*numpy.log10(SNR[0]) - if my_nbeams == 2: - - zenith = numpy.array(zenith) - zenith -= zenithCorrection - zenith *= numpy.pi/180 - if zenithCorrection != 0 : - heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) - else : - heiRang1 = heiRang - velRadial1 = velRadial - SNR1 = SNR - - alp = zenith[0] - bet = zenith[1] - - w_w = velRadial1[0,:] - w_e = velRadial1[1,:] - w_w_err = velRerr[0,:] - w_e_err = velRerr[1,:] - smooth_e = dataOut.clean_num_aver[2,:] - chisq_e = dataOut.data_error[1,0,:] - p_e0 = rbufi[2,:] - p_e1 = rbufi[3,:] - - tini=time.localtime(dataOut.utctime) - - if tini[3] >= 6 and tini[3] < 18 : - w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w) - w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err) - else: - w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w) - w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp) - w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err) - w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp) - w_w = w_wtmp - w_w_err = w_w_errtmp + if tini[3] >= 6 and tini[3] < 18 : + w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w) + w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err) + else: + w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w) + w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp) + w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err) + w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp) + w_w = w_wtmp + w_w_err = w_w_errtmp #if my_nbeams == 2: smooth_eC=ebufc[4,:] @@ -5891,9 +5440,9 @@ class EWDriftsEstimation(Operation): smooth_e[val]=0 val = (numpy.isfinite(smooth_eC)==False).nonzero() smooth_eC[val]=0 - #p_e0 = (p_e0*smooth_e+p_e0C*smooth_eC)/(smooth_e+smooth_eC) - #p_e1 = (p_e1*smooth_e+p_e1C*smooth_eC)/(smooth_e+smooth_eC) - + p_e0_a = (p_e0*smooth_e+p_e0C*smooth_eC)/(smooth_e+smooth_eC) + p_e1_a = (p_e1*smooth_e+p_e1C*smooth_eC)/(smooth_e+smooth_eC) + #print(w_e,w_eC,p_e0C,sat_fits[2,:]) if len(sat_fits) >0 : p_e0C = p_e0C + sat_fits[2,:] p_e1C = p_e1C + sat_fits[3,:] @@ -5914,6 +5463,22 @@ class EWDriftsEstimation(Operation): w_err = numpy.sqrt((w_w_err*numpy.sin(bet))**2.+(w_e_err*numpy.sin(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp)) u_err = numpy.sqrt((w_w_err*numpy.cos(bet))**2.+(w_e_err*numpy.cos(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp)) + + #wtmp = numpy.where(w < -100 ,numpy.nan,w) + wtmp = numpy.where(abs(w) > 100 ,numpy.nan,w) + werrtmp = numpy.where(abs(w_err) > 100 ,numpy.nan,w_err) + werrtmp = numpy.where(numpy.isfinite(wtmp)==True,werrtmp,numpy.nan) + #wtmp = numpy.where(numpy.isfinite(werrtmp)==True,wtmp,numpy.nan) + + #utmp = numpy.where(u < -500,numpy.nan,u) + utmp = numpy.where(abs(u) > 500 ,numpy.nan,u) + uerrtmp = numpy.where(abs(u_err) > 500 ,numpy.nan,u_err) + uerrtmp = numpy.where(numpy.isfinite(utmp)==True,uerrtmp,numpy.nan) + #utmp = numpy.where(numpy.isfinite(uerrtmp)==True,utmp,numpy.nan) + w= wtmp.copy() + u= utmp.copy() + w_err= werrtmp.copy() + u_err= uerrtmp.copy() winds = numpy.vstack((w,u)) dataOut.heightList = heiRang1 @@ -5924,7 +5489,7 @@ class EWDriftsEstimation(Operation): #print('alt ',range1*numpy.sin(86.1*numpy.pi/180)) #print(numpy.min([dataOut.eldir7,dataOut.eldir8])) galt = range1*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.) - dataOut.params = numpy.vstack((range1,galt,w,w_err,u,u_err,w_w,w_w_err,w_e,w_e_err,numpy.log10(p_w0),numpy.log10(p_w0C),numpy.log10(p_w1),numpy.log10(p_w1C),numpy.log10(p_e0),numpy.log10(p_e0C),numpy.log10(p_e1),numpy.log10(p_e1C),chisq_w,chisq_e)) + dataOut.params = numpy.vstack((range1,galt,w,w_err,u,u_err,w_w,w_w_err,w_e,w_e_err,numpy.log10(p_w0),numpy.log10(p_w0C),numpy.log10(p_w1),numpy.log10(p_w1C),numpy.log10(p_e0),numpy.log10(p_e0C),numpy.log10(p_e1),numpy.log10(p_e1C),chisq_w,chisq_e,numpy.log10(p_w0_a),numpy.log10(p_w1_a),numpy.log10(p_e0_a),numpy.log10(p_e1_a))) #snr1 = 10*numpy.log10(SNR1[0]) #print(min(snr1), max(snr1)) snr1 = numpy.vstack((p_w0,p_w1,p_e0,p_e1)) @@ -6001,39 +5566,43 @@ class EWDriftsEstimation(Operation): uA_err[nhei_avg-1] = sigma uA = uA[6*h_avgs:nhei_avg-0] uA_err = uA_err[6*h_avgs:nhei_avg-0] - dataOut.drifts_avg = numpy.vstack((wA,uA)) + dataOut.drifts_avg = numpy.vstack((wA,uA)) if my_nbeams == 1: dataOut.drifts_avg = wA #deltahavg= wA*0.0+deltah dataOut.range = range1 galtavg = range_aver*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.) dataOut.params_avg = numpy.vstack((wA,wA_err,uA,uA_err,range_aver,galtavg,delta_h)) - - #print('comparando dim de avg ',wA.shape,deltahavg.shape,range_aver.shape) + + ''' tini=time.localtime(dataOut.utctime) datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2) nfile = fileDrifts+'/jro'+datefile+'drifts_sch3.txt' - + #nfile = '/home/pcondor/Database/ewdriftsschain2019/jro'+datefile+'drifts_sch3.txt' + #print(len(dataOut.drifts_avg),dataOut.drifts_avg.shape) f1 = open(nfile,'a') datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4]) driftavgstr=str(dataOut.drifts_avg) numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i') numpy.savetxt(f1,numpy.reshape(range_aver,(1,len(range_aver))) ,fmt='%10.2f') + #numpy.savetxt(f1,numpy.reshape(dataOut.drifts_avg,(7,len(dataOut.drifts_avg))) ,fmt='%10.2f') numpy.savetxt(f1,dataOut.drifts_avg[:,:],fmt='%10.2f') f1.close() - + ''' + ''' swfile = fileDrifts+'/jro'+datefile+'drifts_sw.txt' f1 = open(swfile,'a') numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i') numpy.savetxt(f1,numpy.reshape(heiRang,(1,len(heiRang))),fmt='%10.2f') numpy.savetxt(f1,dataOut.data_param[:,0,:],fmt='%10.2f') f1.close() - dataOut.heightListtmp = dataOut.heightList ''' + dataOut.heightListtmp = dataOut.heightList + #Envio data de drifts a mysql fechad = str(tini[0]).zfill(4)+'-'+str(tini[1]).zfill(2)+'-'+str(tini[2]).zfill(2)+' '+str(tini[3]).zfill(2)+':'+str(tini[4]).zfill(2)+':'+str(0).zfill(2) mydb = mysql.connector.connect( - host="10.10.110.213", + host="10.10.110.205", user="user_clima", password="5D.bh(B2)Y_wRNz9", database="clima_espacial" @@ -6056,7 +5625,7 @@ class EWDriftsEstimation(Operation): mydb.commit() print(mycursor.rowcount, "record inserted.") - ''' + return dataOut class setHeightDrifts(Operation): @@ -7317,7 +6886,6 @@ class SMOperations(): chan1Y = chan1 chan2Y = chan2 distances[2:4] = numpy.array([d1,d2]) - pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)] return pairslist, distances diff --git a/schainpy/scripts/comparar2.py b/schainpy/scripts/comparar2.py new file mode 100644 index 0000000..048b698 --- /dev/null +++ b/schainpy/scripts/comparar2.py @@ -0,0 +1,663 @@ +#primer windprofiler +class WindProfiler(Operation): + + __isConfig = False + + __initime = None + __lastdatatime = None + __integrationtime = None + + __buffer = None + + __dataReady = False + + __firstdata = None + + n = None + + def __init__(self): + Operation.__init__(self) + + def __calculateCosDir(self, elev, azim): + zen = (90 - elev)*numpy.pi/180 + azim = azim*numpy.pi/180 + cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) + cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2) + + signX = numpy.sign(numpy.cos(azim)) + signY = numpy.sign(numpy.sin(azim)) + + cosDirX = numpy.copysign(cosDirX, signX) + cosDirY = numpy.copysign(cosDirY, signY) + return cosDirX, cosDirY + + def __calculateAngles(self, theta_x, theta_y, azimuth): + + dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2) + zenith_arr = numpy.arccos(dir_cosw) + azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180 + + dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr) + dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr) + + return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw + + def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly): + + if horOnly: + A = numpy.c_[dir_cosu,dir_cosv] + else: + A = numpy.c_[dir_cosu,dir_cosv,dir_cosw] + A = numpy.asmatrix(A) + A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose() + + return A1 + + def __correctValues(self, heiRang, phi, velRadial, SNR): + listPhi = phi.tolist() + maxid = listPhi.index(max(listPhi)) + minid = listPhi.index(min(listPhi)) + + rango = list(range(len(phi))) + + heiRang1 = heiRang*math.cos(phi[maxid]) + heiRangAux = heiRang*math.cos(phi[minid]) + indOut = (heiRang1 < heiRangAux[0]).nonzero() + heiRang1 = numpy.delete(heiRang1,indOut) + + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) + SNR1 = numpy.zeros([len(phi),len(heiRang1)]) + + for i in rango: + x = heiRang*math.cos(phi[i]) + y1 = velRadial[i,:] + f1 = interpolate.interp1d(x,y1,kind = 'cubic') + + x1 = heiRang1 + y11 = f1(x1) + + y2 = SNR[i,:] + f2 = interpolate.interp1d(x,y2,kind = 'cubic') + y21 = f2(x1) + + velRadial1[i,:] = y11 + SNR1[i,:] = y21 + + return heiRang1, velRadial1, SNR1 + + def __calculateVelUVW(self, A, velRadial): + + #Operacion Matricial + velUVW = numpy.zeros((A.shape[0],velRadial.shape[1])) + velUVW[:,:] = numpy.dot(A,velRadial) + + + return velUVW + + def techniqueDBS(self, kwargs): + """ + Function that implements Doppler Beam Swinging (DBS) technique. + + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, + Direction correction (if necessary), Ranges and SNR + + Output: Winds estimation (Zonal, Meridional and Vertical) + + Parameters affected: Winds, height range, SNR + """ + velRadial0 = kwargs['velRadial'] + heiRang = kwargs['heightList'] + SNR0 = kwargs['SNR'] + + if 'dirCosx' in kwargs and 'dirCosy' in kwargs: + theta_x = numpy.array(kwargs['dirCosx']) + theta_y = numpy.array(kwargs['dirCosy']) + else: + elev = numpy.array(kwargs['elevation']) + azim = numpy.array(kwargs['azimuth']) + theta_x, theta_y = self.__calculateCosDir(elev, azim) + azimuth = kwargs['correctAzimuth'] + if 'horizontalOnly' in kwargs: + horizontalOnly = kwargs['horizontalOnly'] + else: horizontalOnly = False + if 'correctFactor' in kwargs: + correctFactor = kwargs['correctFactor'] + else: correctFactor = 1 + if 'channelList' in kwargs: + channelList = kwargs['channelList'] + if len(channelList) == 2: + horizontalOnly = True + arrayChannel = numpy.array(channelList) + param = param[arrayChannel,:,:] + theta_x = theta_x[arrayChannel] + theta_y = theta_y[arrayChannel] + + azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0) + A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly) + + #Calculo de Componentes de la velocidad con DBS + winds = self.__calculateVelUVW(A,velRadial1) + + return winds, heiRang1, SNR1 + + def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None): + + nPairs = len(pairs_ccf) + posx = numpy.asarray(posx) + posy = numpy.asarray(posy) + + #Rotacion Inversa para alinear con el azimuth + if azimuth!= None: + azimuth = azimuth*math.pi/180 + posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth) + posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth) + else: + posx1 = posx + posy1 = posy + + #Calculo de Distancias + distx = numpy.zeros(nPairs) + disty = numpy.zeros(nPairs) + dist = numpy.zeros(nPairs) + ang = numpy.zeros(nPairs) + + for i in range(nPairs): + distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]] + disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]] + dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) + ang[i] = numpy.arctan2(disty[i],distx[i]) + + return distx, disty, dist, ang + #Calculo de Matrices + + + def __calculateVelVer(self, phase, lagTRange, _lambda): + + Ts = lagTRange[1] - lagTRange[0] + velW = -_lambda*phase/(4*math.pi*Ts) + + return velW + + def __calculateVelHorDir(self, dist, tau1, tau2, ang): + nPairs = tau1.shape[0] + nHeights = tau1.shape[1] + vel = numpy.zeros((nPairs,3,nHeights)) + dist1 = numpy.reshape(dist, (dist.size,1)) + + angCos = numpy.cos(ang) + angSin = numpy.sin(ang) + + vel0 = dist1*tau1/(2*tau2**2) + vel[:,0,:] = (vel0*angCos).sum(axis = 1) + vel[:,1,:] = (vel0*angSin).sum(axis = 1) + + ind = numpy.where(numpy.isinf(vel)) + vel[ind] = numpy.nan + + return vel + + def techniqueSA(self, kwargs): + + """ + Function that implements Spaced Antenna (SA) technique. + + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, + Direction correction (if necessary), Ranges and SNR + + Output: Winds estimation (Zonal, Meridional and Vertical) + + Parameters affected: Winds + """ + position_x = kwargs['positionX'] + position_y = kwargs['positionY'] + azimuth = kwargs['azimuth'] + + if 'correctFactor' in kwargs: + correctFactor = kwargs['correctFactor'] + else: + correctFactor = 1 + + groupList = kwargs['groupList'] + pairs_ccf = groupList[1] + tau = kwargs['tau'] + _lambda = kwargs['_lambda'] + + #Cross Correlation pairs obtained + + indtau = tau.shape[0]/2 + tau1 = tau[:indtau,:] + tau2 = tau[indtau:-1,:] + phase1 = tau[-1,:] + + #--------------------------------------------------------------------- + #Metodo Directo + distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth) + winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) + winds = stats.nanmean(winds, axis=0) + #--------------------------------------------------------------------- + #Metodo General + + #--------------------------------------------------------------------- + winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) + winds = correctFactor*winds + return winds + + def __checkTime(self, currentTime, paramInterval, outputInterval): + + dataTime = currentTime + paramInterval + deltaTime = dataTime - self.__initime + + if deltaTime >= outputInterval or deltaTime < 0: + self.__dataReady = True + return + + def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax): + ''' + Function that implements winds estimation technique with detected meteors. + + Input: Detected meteors, Minimum meteor quantity to wind estimation + + Output: Winds estimation (Zonal and Meridional) + + Parameters affected: Winds + ''' + #Settings + nInt = (heightMax - heightMin)/2 + nInt = int(nInt) + winds = numpy.zeros((2,nInt))*numpy.nan + + #Filter errors + error = numpy.where(arrayMeteor[:,-1] == 0)[0] + finalMeteor = arrayMeteor[error,:] + + #Meteor Histogram + finalHeights = finalMeteor[:,2] + hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax)) + nMeteorsPerI = hist[0] + heightPerI = hist[1] + + #Sort of meteors + indSort = finalHeights.argsort() + finalMeteor2 = finalMeteor[indSort,:] + + # Calculating winds + ind1 = 0 + ind2 = 0 + + for i in range(nInt): + nMet = nMeteorsPerI[i] + ind1 = ind2 + ind2 = ind1 + nMet + + meteorAux = finalMeteor2[ind1:ind2,:] + + if meteorAux.shape[0] >= meteorThresh: + vel = meteorAux[:, 6] + zen = meteorAux[:, 4]*numpy.pi/180 + azim = meteorAux[:, 3]*numpy.pi/180 + + n = numpy.cos(zen) + # m = (1 - n**2)/(1 - numpy.tan(azim)**2) + # l = m*numpy.tan(azim) + l = numpy.sin(zen)*numpy.sin(azim) + m = numpy.sin(zen)*numpy.cos(azim) + + A = numpy.vstack((l, m)).transpose() + A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose()) + windsAux = numpy.dot(A1, vel) + + winds[0,i] = windsAux[0] + winds[1,i] = windsAux[1] + + return winds, heightPerI[:-1] + + def techniqueNSM_SA(self, **kwargs): + metArray = kwargs['metArray'] + heightList = kwargs['heightList'] + timeList = kwargs['timeList'] + + rx_location = kwargs['rx_location'] + groupList = kwargs['groupList'] + azimuth = kwargs['azimuth'] + dfactor = kwargs['dfactor'] + k = kwargs['k'] + + azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth) + d = dist*dfactor + #Phase calculation + metArray1 = self.__getPhaseSlope(metArray, heightList, timeList) + + metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities + + velEst = numpy.zeros((heightList.size,2))*numpy.nan + azimuth1 = azimuth1*numpy.pi/180 + + for i in range(heightList.size): + h = heightList[i] + indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0] + metHeight = metArray1[indH,:] + if metHeight.shape[0] >= 2: + velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities + iazim = metHeight[:,1].astype(int) + azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths + A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux))) + A = numpy.asmatrix(A) + A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose() + velHor = numpy.dot(A1,velAux) + + velEst[i,:] = numpy.squeeze(velHor) + return velEst + + def __getPhaseSlope(self, metArray, heightList, timeList): + meteorList = [] + #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2 + #Putting back together the meteor matrix + utctime = metArray[:,0] + uniqueTime = numpy.unique(utctime) + + phaseDerThresh = 0.5 + ippSeconds = timeList[1] - timeList[0] + sec = numpy.where(timeList>1)[0][0] + nPairs = metArray.shape[1] - 6 + nHeights = len(heightList) + + for t in uniqueTime: + metArray1 = metArray[utctime==t,:] +# phaseDerThresh = numpy.pi/4 #reducir Phase thresh + tmet = metArray1[:,1].astype(int) + hmet = metArray1[:,2].astype(int) + + metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1)) + metPhase[:,:] = numpy.nan + metPhase[:,hmet,tmet] = metArray1[:,6:].T + + #Delete short trails + metBool = ~numpy.isnan(metPhase[0,:,:]) + heightVect = numpy.sum(metBool, axis = 1) + metBool[heightVect phaseDerThresh)) + metPhase[phDerAux] = numpy.nan + + #--------------------------METEOR DETECTION ----------------------------------------- + indMet = numpy.where(numpy.any(metBool,axis=1))[0] + + for p in numpy.arange(nPairs): + phase = metPhase[p,:,:] + phDer = metDer[p,:,:] + + for h in indMet: + height = heightList[h] + phase1 = phase[h,:] #82 + phDer1 = phDer[h,:] + + phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap + + indValid = numpy.where(~numpy.isnan(phase1))[0] + initMet = indValid[0] + endMet = 0 + + for i in range(len(indValid)-1): + + #Time difference + inow = indValid[i] + inext = indValid[i+1] + idiff = inext - inow + #Phase difference + phDiff = numpy.abs(phase1[inext] - phase1[inow]) + + if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor + sizeTrail = inow - initMet + 1 + if sizeTrail>3*sec: #Too short meteors + x = numpy.arange(initMet,inow+1)*ippSeconds + y = phase1[initMet:inow+1] + ynnan = ~numpy.isnan(y) + x = x[ynnan] + y = y[ynnan] + slope, intercept, r_value, p_value, std_err = stats.linregress(x,y) + ylin = x*slope + intercept + rsq = r_value**2 + if rsq > 0.5: + vel = slope#*height*1000/(k*d) + estAux = numpy.array([utctime,p,height, vel, rsq]) + meteorList.append(estAux) + initMet = inext + metArray2 = numpy.array(meteorList) + + return metArray2 + + def __calculateAzimuth1(self, rx_location, pairslist, azimuth0): + + azimuth1 = numpy.zeros(len(pairslist)) + dist = numpy.zeros(len(pairslist)) + + for i in range(len(rx_location)): + ch0 = pairslist[i][0] + ch1 = pairslist[i][1] + + diffX = rx_location[ch0][0] - rx_location[ch1][0] + diffY = rx_location[ch0][1] - rx_location[ch1][1] + azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi + dist[i] = numpy.sqrt(diffX**2 + diffY**2) + + azimuth1 -= azimuth0 + return azimuth1, dist + + def techniqueNSM_DBS(self, **kwargs): + metArray = kwargs['metArray'] + heightList = kwargs['heightList'] + timeList = kwargs['timeList'] + azimuth = kwargs['azimuth'] + theta_x = numpy.array(kwargs['theta_x']) + theta_y = numpy.array(kwargs['theta_y']) + + utctime = metArray[:,0] + cmet = metArray[:,1].astype(int) + hmet = metArray[:,3].astype(int) + SNRmet = metArray[:,4] + vmet = metArray[:,5] + spcmet = metArray[:,6] + + nChan = numpy.max(cmet) + 1 + nHeights = len(heightList) + + azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) + hmet = heightList[hmet] + h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights + + velEst = numpy.zeros((heightList.size,2))*numpy.nan + + for i in range(nHeights - 1): + hmin = heightList[i] + hmax = heightList[i + 1] + + thisH = (h1met>=hmin) & (h1met8) & (vmet<50) & (spcmet<10) + indthisH = numpy.where(thisH) + + if numpy.size(indthisH) > 3: + + vel_aux = vmet[thisH] + chan_aux = cmet[thisH] + cosu_aux = dir_cosu[chan_aux] + cosv_aux = dir_cosv[chan_aux] + cosw_aux = dir_cosw[chan_aux] + + nch = numpy.size(numpy.unique(chan_aux)) + if nch > 1: + A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True) + velEst[i,:] = numpy.dot(A,vel_aux) + + return velEst + + def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs): + + param = dataOut.data_param + #if dataOut.abscissaList != None: + if numpy.any(dataOut.abscissaList): + absc = dataOut.abscissaList[:-1] + # noise = dataOut.noise + heightList = dataOut.heightList + SNR = dataOut.data_snr + + if technique == 'DBS': + + kwargs['velRadial'] = param[:,1,:] #Radial velocity + kwargs['heightList'] = heightList + kwargs['SNR'] = SNR + + dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function + dataOut.utctimeInit = dataOut.utctime + dataOut.outputInterval = dataOut.paramInterval + + elif technique == 'SA': + + #Parameters +# position_x = kwargs['positionX'] +# position_y = kwargs['positionY'] +# azimuth = kwargs['azimuth'] +# +# if kwargs.has_key('crosspairsList'): +# pairs = kwargs['crosspairsList'] +# else: +# pairs = None +# +# if kwargs.has_key('correctFactor'): +# correctFactor = kwargs['correctFactor'] +# else: +# correctFactor = 1 + +# tau = dataOut.data_param +# _lambda = dataOut.C/dataOut.frequency +# pairsList = dataOut.groupList +# nChannels = dataOut.nChannels + + kwargs['groupList'] = dataOut.groupList + kwargs['tau'] = dataOut.data_param + kwargs['_lambda'] = dataOut.C/dataOut.frequency +# dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) + dataOut.data_output = self.techniqueSA(kwargs) + dataOut.utctimeInit = dataOut.utctime + dataOut.outputInterval = dataOut.timeInterval + + elif technique == 'Meteors': + dataOut.flagNoData = True + self.__dataReady = False + + if 'nHours' in kwargs: + nHours = kwargs['nHours'] + else: + nHours = 1 + + if 'meteorsPerBin' in kwargs: + meteorThresh = kwargs['meteorsPerBin'] + else: + meteorThresh = 6 + + if 'hmin' in kwargs: + hmin = kwargs['hmin'] + else: hmin = 70 + if 'hmax' in kwargs: + hmax = kwargs['hmax'] + else: hmax = 110 + + dataOut.outputInterval = nHours*3600 + + if self.__isConfig == False: +# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) + #Get Initial LTC time + self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) + self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() + + self.__isConfig = True + + if self.__buffer is None: + self.__buffer = dataOut.data_param + self.__firstdata = copy.copy(dataOut) + + else: + self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) + + self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready + + if self.__dataReady: + dataOut.utctimeInit = self.__initime + + self.__initime += dataOut.outputInterval #to erase time offset + + dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) + dataOut.flagNoData = False + self.__buffer = None + + elif technique == 'Meteors1': + dataOut.flagNoData = True + self.__dataReady = False + + if 'nMins' in kwargs: + nMins = kwargs['nMins'] + else: nMins = 20 + if 'rx_location' in kwargs: + rx_location = kwargs['rx_location'] + else: rx_location = [(0,1),(1,1),(1,0)] + if 'azimuth' in kwargs: + azimuth = kwargs['azimuth'] + else: azimuth = 51.06 + if 'dfactor' in kwargs: + dfactor = kwargs['dfactor'] + if 'mode' in kwargs: + mode = kwargs['mode'] + if 'theta_x' in kwargs: + theta_x = kwargs['theta_x'] + if 'theta_y' in kwargs: + theta_y = kwargs['theta_y'] + else: mode = 'SA' + + #Borrar luego esto + if dataOut.groupList is None: + dataOut.groupList = [(0,1),(0,2),(1,2)] + groupList = dataOut.groupList + C = 3e8 + freq = 50e6 + lamb = C/freq + k = 2*numpy.pi/lamb + + timeList = dataOut.abscissaList + heightList = dataOut.heightList + + if self.__isConfig == False: + dataOut.outputInterval = nMins*60 +# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) + #Get Initial LTC time + initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) + minuteAux = initime.minute + minuteNew = int(numpy.floor(minuteAux/nMins)*nMins) + self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() + + self.__isConfig = True + + if self.__buffer is None: + self.__buffer = dataOut.data_param + self.__firstdata = copy.copy(dataOut) + + else: + self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) + + self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready + + if self.__dataReady: + dataOut.utctimeInit = self.__initime + self.__initime += dataOut.outputInterval #to erase time offset + + metArray = self.__buffer + if mode == 'SA': + dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList) + elif mode == 'DBS': + dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y) + dataOut.data_output = dataOut.data_output.T + dataOut.flagNoData = False + self.__buffer = None + + return