#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