schain Commit - r842:ae31287a1940 · Jicamarca Repository

112

self.dataOut.utctimeInit = self.dataIn.utctime

112

self.dataOut.utctimeInit = self.dataIn.utctime

113

self.dataOut.paramInterval = self.dataIn.timeInterval

113

self.dataOut.paramInterval = self.dataIn.timeInterval

114

115

#------------------- Get Moments ----------------------------------

115

return

116

def GetMoments(self, channelList = None):

116

117

class SpectralMoments(Operation):

118

117

'''

119

'''

118

Function ~~Get~~Moments()

120

Function SpectralMoments()

121

122

Calculates moments (power, mean, standard deviation) and SNR of the signal

123

124

Type of dataIn: Spectra

119

125

120

Input:

126

Input:

121

channelList : simple channel list to select e.g. [2,3,7]

127

channelList : simple channel list to select e.g. [2,3,7]

122

self.dataOut.data_pre

128

self.dataOut.data_pre : Spectral data

123

self.dataOut.abscissaList

129

self.dataOut.abscissaList : List of frequencies

124

self.dataOut.noise

130

self.dataOut.noise : Noise level per channel

125

131

126

Affected:

132

Affected:

127

self.dataOut.data_param

133

self.dataOut.data_param : Parameters per channel

128

self.dataOut.data_SNR

134

self.dataOut.data_SNR : SNR per channel

129

135

130

'''

136

'''

131

self.dataOut.data_pre = self.dataOut.data_pre[0]

137

132

data = self.dataOut.data_pre

138

def run(self, dataOut, channelList=None):

133

absc = self.dataOut.abscissaList[:-1]

139

134

noise = self.dataOut.noise

140

dataOut.data_pre = dataOut.data_pre[0]

141

data = dataOut.data_pre

142

absc = dataOut.abscissaList[:-1]

143

noise = dataOut.noise

135

144

136

data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))

145

data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))

137

146

138

if channelList== None:

147

if channelList== None:

139

channelList = ~~self~~.~~dataIn~~.channelList

148

channelList = dataOut.channelList

140

~~self~~.dataOut.channelList = channelList

149

dataOut.channelList = channelList

141

150

142

for ind in channelList:

151

for ind in channelList:

143

data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])

152

data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])

144

153

145

~~self~~.dataOut.data_param = data_param[:,1:,:]

154

dataOut.data_param = data_param[:,1:,:]

146

~~self~~.dataOut.data_SNR = data_param[:,0]

155

dataOut.data_SNR = data_param[:,0]

147

return

156

return

148

157

149

def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):

158

def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):

221

#------------------ Get SA Parameters --------------------------

230

#------------------ Get SA Parameters --------------------------

222

231

223

def GetSAParameters(self):

232

def GetSAParameters(self):

233

#SA en frecuencia

224

pairslist = self.dataOut.groupList

234

pairslist = self.dataOut.groupList

225

num_pairs = len(pairslist)

235

num_pairs = len(pairslist)

226

236

249

cspc_par[i,:,:] = __calculateMoments(norm_cspectra)

259

cspc_par[i,:,:] = __calculateMoments(norm_cspectra)

250

#------------------- Get Lags ----------------------------------

260

#------------------- Get Lags ----------------------------------

251

261

252

def GetLags(self):

262

class SALags(Operation):

253

'''

263

'''

254

Function GetMoments()

264

Function GetMoments()

255

265

266

self.dataOut.data_param

276

self.dataOut.data_param

267

277

268

'''

278

'''

279

def run(self, dataOut):

280

data = dataOut.data_pre

281

data_acf = dataOut.data_pre[0]

282

data_ccf = dataOut.data_pre[1]

269

283

270

~~data~~ = ~~self~~.dataOut.~~data_pre~~

284

normFactor = dataOut.normFactor

271

normFactor = self.dataOut.normFactor

285

nHeights = dataOut.nHeights

272

nHeights = self.dataOut.nHeights

286

absc = dataOut.abscissaList[:-1]

273

absc = self.dataOut.abscissaList[:-1]

287

noise = dataOut.noise

274

noise = self.dataOut.noise

288

SNR = dataOut.data_SNR

275

SNR = self.dataOut.data_SNR

289

# pairsList = dataOut.groupList

276

pairsList = self.dataOut.groupList

290

nChannels = dataOut.nChannels

277

nChannels = self.dataOut.nChannels

291

# pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

278

pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

292

dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))

279

self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))

280

293

281

dataNorm = numpy.abs(data)

294

dataNorm = numpy.abs(data)

282

for l in range(len(pairsList)):

295

for l in range(len(pairsList)):

283

dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]

296

dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]

284

297

285

self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)

298

self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)

286

self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, ~~pairsAutoCorr~~, absc)

299

self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)

287

return

300

return

288

301

289

def __getPairsAutoCorr(self, pairsList, nChannels):

302

# def __getPairsAutoCorr(self, pairsList, nChannels):

290

303

#

291

pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

304

# pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

292

305

#

293

for l in range(len(pairsList)):

306

# for l in range(len(pairsList)):

294

firstChannel = pairsList[l][0]

307

# firstChannel = pairsList[l][0]

295

secondChannel = pairsList[l][1]

308

# secondChannel = pairsList[l][1]

296

309

#

297

#Obteniendo pares de Autocorrelacion

310

# #Obteniendo pares de Autocorrelacion

298

if firstChannel == secondChannel:

311

# if firstChannel == secondChannel:

299

pairsAutoCorr[firstChannel] = int(l)

312

# pairsAutoCorr[firstChannel] = int(l)

300

313

#

301

pairsAutoCorr = pairsAutoCorr.astype(int)

314

# pairsAutoCorr = pairsAutoCorr.astype(int)

302

315

#

303

pairsCrossCorr = range(len(pairsList))

316

# pairsCrossCorr = range(len(pairsList))

304

pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

317

# pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

305

318

#

306

return pairsAutoCorr, pairsCrossCorr

319

# return pairsAutoCorr, pairsCrossCorr

307

320

308

def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):

321

def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):

309

322

331

344

332

return tau

345

return tau

333

346

334

def __calculateLag1Phase(self, data, ~~pairs~~, lagTRange):

347

def __calculateLag1Phase(self, data, lagTRange):

335

data1 = stats.nanmean(data[~~pairs,:,:]~~, axis = 0)

348

data1 = stats.nanmean(data, axis = 0)

336

lag1 = numpy.where(lagTRange == 0)[0][0] + 1

349

lag1 = numpy.where(lagTRange == 0)[0][0] + 1

337

350

338

phase = numpy.angle(data1[lag1,:])

351

phase = numpy.angle(data1[lag1,:])

339

352

340

return phase

353

return phase

341

#------------------- Detect Meteors ------------------------------

342

343

def MeteorDetection(self, hei_ref = None, tauindex = 0,

344

phaseOffsets = None,

345

cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,

346

noise_timeStep = 4, noise_multiple = 4,

347

multDet_timeLimit = 1, multDet_rangeLimit = 3,

348

phaseThresh = 20, SNRThresh = 5,

349

hmin = 50, hmax=150, azimuth = 0,

350

# channel25X = 0, channel20X = 4, channelCentX = 2,

351

# channel25Y = 1, channel20Y = 3, channelCentY = 2,

352

channelPositions = None) :

353

354

355

class SpectralFitting(Operation):

354

'''

356

'''

355

Function ~~DetectMeteor~~s()

357

Function GetMoments()

356

Project developed with paper:

357

HOLDSWORTH ET AL. 2004

358

359

Input:

359

Input:

360

self.dataOut.data_pre

360

Output:

361

Variables modified:

362

'''

361

363

362

centerReceiverIndex: From the channels, which is the center receiver

364

def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):

363

365

364

hei_ref: Height reference for the Beacon signal extraction

365

tauindex:

366

predefinedPhaseShifts: Predefined phase offset for the voltge signals

367

366

368

cohDetection: Whether to user Coherent detection or not

367

if path != None:

369

cohDet_timeStep: Coherent Detection calculation time step

368

sys.path.append(path)

370

cohDet_thresh: Coherent Detection phase threshold to correct phases

369

self.dataOut.library = importlib.import_module(file)

371

370

372

noise_timeStep: Noise calculation time step

371

#To be inserted as a parameter

373

noise_multiple: Noise multiple to define signal threshold

372

groupArray = numpy.array(groupList)

373

# groupArray = numpy.array([[0,1],[2,3]])

374

self.dataOut.groupList = groupArray

374

375

multDet_timeLimit: Multiple Detection Removal time limit in seconds

376

nGroups = groupArray.shape[0]

376

multDet_rangeLimit: Multiple Detection Removal range limit in km

377

nChannels = self.dataIn.nChannels

378

nHeights=self.dataIn.heightList.size

377

379

378

phaseThresh: Maximum phase difference between receiver to be consider a meteor

380

#Parameters Array

379

SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor

381

self.dataOut.data_param = None

380

382

381

hmin: Minimum Height of the meteor to use it in the further wind estimations

383

#Set constants

382

hmax: Maximum Height of the meteor to use it in the further wind estimations

384

constants = self.dataOut.library.setConstants(self.dataIn)

383

azimuth: Azimuth angle correction

385

self.dataOut.constants = constants

386

M = self.dataIn.normFactor

387

N = self.dataIn.nFFTPoints

388

ippSeconds = self.dataIn.ippSeconds

389

K = self.dataIn.nIncohInt

390

pairsArray = numpy.array(self.dataIn.pairsList)

384

391

385

Affected:

392

#List of possible combinations

386

self.dataOut.data_param

393

listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)

394

indCross = numpy.zeros(len(list(listComb)), dtype = 'int')

387

395

388

Rejection Criteria (Errors):

396

if getSNR:

389

0: No error; analysis OK

397

listChannels = groupArray.reshape((groupArray.size))

390

1: SNR < SNR threshold

398

listChannels.sort()

391

2: angle of arrival (AOA) ambiguously determined

399

noise = self.dataIn.getNoise()

392

3: AOA estimate not feasible

400

self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])

393

4: Large difference in AOAs obtained from different antenna baselines

394

5: echo at start or end of time series

395

6: echo less than 5 examples long; too short for analysis

396

7: echo rise exceeds 0.3s

397

8: echo decay time less than twice rise time

398

9: large power level before echo

399

10: large power level after echo

400

11: poor fit to amplitude for estimation of decay time

401

12: poor fit to CCF phase variation for estimation of radial drift velocity

402

13: height unresolvable echo: not valid height within 70 to 110 km

403

14: height ambiguous echo: more then one possible height within 70 to 110 km

404

15: radial drift velocity or projected horizontal velocity exceeds 200 m/s

405

16: oscilatory echo, indicating event most likely not an underdense echo

406

401

407

17: phase difference in meteor Reestimation

402

for i in range(nGroups):

403

coord = groupArray[i,:]

408

404

409

Data Storage:

405

#Input data array

410

Meteors for Wind Estimation (8):

406

data = self.dataIn.data_spc[coord,:,:]/(M*N)

411

Utc Time | Range Height

407

data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))

412

Azimuth Zenith errorCosDir

413

VelRad errorVelRad

414

Phase0 Phase1 Phase2 Phase3

415

TypeError

416

408

417

'''

409

#Cross Spectra data array for Covariance Matrixes

418

#Getting Pairslist

410

ind = 0

419

if channelPositions == None:

411

for pairs in listComb:

420

# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T

412

pairsSel = numpy.array([coord[x],coord[y]])

421

channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella

413

indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])

422

meteorOps = MeteorOperations()

414

ind += 1

423

pairslist0, distances = meteorOps.getPhasePairs(channelPositions)

415

dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)

416

dataCross = dataCross**2/K

424

417

425

#Get Beacon signal

418

for h in range(nHeights):

426

newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])

419

# print self.dataOut.heightList[h]

427

420

428

if hei_ref != None:

421

#Input

429

newheis = numpy.where(self.dataOut.heightList>hei_ref)

422

d = data[:,h]

430

423

431

heiRang = self.dataOut.getHeiRange()

424

#Covariance Matrix

425

D = numpy.diag(d**2/K)

426

ind = 0

427

for pairs in listComb:

428

#Coordinates in Covariance Matrix

429

x = pairs[0]

430

y = pairs[1]

431

#Channel Index

432

S12 = dataCross[ind,:,h]

433

D12 = numpy.diag(S12)

434

#Completing Covariance Matrix with Cross Spectras

435

D[x*N:(x+1)*N,y*N:(y+1)*N] = D12

436

D[y*N:(y+1)*N,x*N:(x+1)*N] = D12

437

ind += 1

438

Dinv=numpy.linalg.inv(D)

439

L=numpy.linalg.cholesky(Dinv)

440

LT=L.T

432

441

433

# nChannel = self.dataOut.nChannels

442

dp = numpy.dot(LT,d)

434

# for i in range(nChannel):

435

# if i != centerReceiverIndex:

436

# pairslist.append((centerReceiverIndex,i))

437

443

438

#****************REMOVING HARDWARE PHASE DIFFERENCES***************

444

#Initial values

439

# see if the user put in pre defined phase shifts

445

data_spc = self.dataIn.data_spc[coord,:,h]

440

voltsPShift = self.dataOut.data_pre.copy()

441

446

442

# if predefinedPhaseShifts != None:

447

if (h>0)and(error1[3]<5):

443

# hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180

448

p0 = self.dataOut.data_param[i,:,h-1]

444

#

449

else:

445

# # elif beaconPhaseShifts:

450

p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))

446

# # #get hardware phase shifts using beacon signal

447

# # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)

448

# # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)

449

#

450

# else:

451

# hardwarePhaseShifts = numpy.zeros(5)

452

#

453

# voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')

454

# for i in range(self.dataOut.data_pre.shape[0]):

455

# voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])

456

451

457

#******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********

452

try:

453

#Least Squares

454

minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)

455

# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))

456

#Chi square error

457

error0 = numpy.sum(infodict['fvec']**2)/(2*N)

458

#Error with Jacobian

459

error1 = self.dataOut.library.errorFunction(minp,constants,LT)

460

except:

461

minp = p0*numpy.nan

462

error0 = numpy.nan

463

error1 = p0*numpy.nan

458

464

459

#Remove DC

465

#Save

460

voltsDC = numpy.mean(voltsPShift,1)

466

if self.dataOut.data_param == None:

461

voltsDC = numpy.mean(voltsDC,1)

467

self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan

462

for i in range(voltsDC.shape[0]):

468

self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan

463

voltsPShift[i] = voltsPShift[i] - voltsDC[i]

464

469

465

#Don't considerate last heights, theyre used to calculate Hardware Phase Shift

470

self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))

466

voltsPShift = voltsPShift[:,:,:newheis[0][0]]

471

self.dataOut.data_param[i,:,h] = minp

472

return

467

473

468

#************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********

474

def __residFunction(self, p, dp, LT, constants):

469

#Coherent Detection

470

if cohDetection:

471

#use coherent detection to get the net power

472

cohDet_thresh = cohDet_thresh*numpy.pi/180

473

voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh)

474

475

#Non-coherent detection!

476

fm = self.dataOut.library.modelFunction(p, constants)

476

powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)

477

fmp=numpy.dot(LT,fm)

477

#********** END OF COH/NON-COH POWER CALCULATION**********************

478

479

#********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************

479

return dp-fmp

480

#Get noise

481

noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)

482

# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)

483

#Get signal threshold

484

signalThresh = noise_multiple*noise

485

#Meteor echoes detection

486

listMeteors = self.__findMeteors(powerNet, signalThresh)

487

#******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********

488

480

489

#************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************

481

def __getSNR(self, z, noise):

490

#Parameters

491

heiRange = self.dataOut.getHeiRange()

492

rangeInterval = heiRange[1] - heiRange[0]

493

rangeLimit = multDet_rangeLimit/rangeInterval

494

timeLimit = multDet_timeLimit/self.dataOut.timeInterval

495

#Multiple detection removals

496

listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)

497

#************ END OF REMOVE MULTIPLE DETECTIONS **********************

498

482

499

#********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************

483

avg = numpy.average(z, axis=1)

500

#Parameters

484

SNR = (avg.T-noise)/noise

501

phaseThresh = phaseThresh*numpy.pi/180

485

SNR = SNR.T

502

thresh = [phaseThresh, noise_multiple, SNRThresh]

486

return SNR

503

#Meteor reestimation (Errors N 1, 6, 12, 17)

504

listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)

505

# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)

506

#Estimation of decay times (Errors N 7, 8, 11)

507

listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)

508

#******************* END OF METEOR REESTIMATION *******************

509

487

510

#********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************

488

def __chisq(p,chindex,hindex):

511

#Calculating Radial Velocity (Error N 15)

489

#similar to Resid but calculates CHI**2

512

radialStdThresh = 10

490

[LT,d,fm]=setupLTdfm(p,chindex,hindex)

513

listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval)

491

dp=numpy.dot(LT,d)

492

fmp=numpy.dot(LT,fm)

493

chisq=numpy.dot((dp-fmp).T,(dp-fmp))

494

return chisq

514

495

515

if len(listMeteors4) > 0:

496

class WindProfiler(Operation):

516

#Setting New Array

517

date = self.dataOut.utctime

518

arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)

519

497

520

#Correcting phase offset

498

__isConfig = False

521

if phaseOffsets != None:

522

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

523

arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

524

499

525

#Second Pairslist

500

__initime = None

526

pairsList = []

501

__lastdatatime = None

527

pairx = (0,1)

502

__integrationtime = None

528

pairy = (2,3)

529

pairsList.append(pairx)

530

pairsList.append(pairy)

531

503

532

jph = numpy.array([0,0,0,0])

504

__buffer = None

533

h = (hmin,hmax)

534

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

535

505

536

# #Calculate AOA (Error N 3, 4)

506

__dataReady = False

537

# #JONES ET AL. 1998

538

# error = arrayParameters[:,-1]

539

# AOAthresh = numpy.pi/8

540

# phases = -arrayParameters[:,9:13]

541

# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)

542

#

543

# #Calculate Heights (Error N 13 and 14)

544

# error = arrayParameters[:,-1]

545

# Ranges = arrayParameters[:,2]

546

# zenith = arrayParameters[:,5]

547

# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)

548

# error = arrayParameters[:,-1]

549

#********************* END OF PARAMETERS CALCULATION **************************

550

507

551

#***************************+ PASS DATA TO NEXT STEP **********************

508

__firstdata = None

552

# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))

553

self.dataOut.data_param = arrayParameters

554

509

555

if arrayParameters == None:

510

n = None

556

self.dataOut.flagNoData = True

557

else:

558

self.dataOut.flagNoData = True

559

511

560

return

512

def __init__(self):

513

Operation.__init__(self)

561

514

562

def CorrectMeteorPhases(self, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):

515

def __calculateCosDir(self, elev, azim):

516

zen = (90 - elev)*numpy.pi/180

517

azim = azim*numpy.pi/180

518

cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))

519

cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)

563

520

564

arrayParameters = self.dataOut.data_param

521

signX = numpy.sign(numpy.cos(azim))

565

pairsList = []

522

signY = numpy.sign(numpy.sin(azim))

566

pairx = (0,1)

567

pairy = (2,3)

568

pairsList.append(pairx)

569

pairsList.append(pairy)

570

jph = numpy.zeros(4)

571

523

572

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

524

cosDirX = numpy.copysign(cosDirX, signX)

573

# arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

525

cosDirY = numpy.copysign(cosDirY, signY)

574

arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))

526

return cosDirX, cosDirY

575

527

576

meteorOps = MeteorOperations()

528

def __calculateAngles(self, theta_x, theta_y, azimuth):

577

if channelPositions == None:

578

# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T

579

channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella

580

529

581

pairslist0, distances = meteorOps.getPhasePairs(channelPositions)

530

dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)

582

h = (hmin,hmax)

531

zenith_arr = numpy.arccos(dir_cosw)

532

azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180

583

533

584

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

534

dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)

535

dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)

585

536

586

self.dataOut.data_param = arrayParameters

537

return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw

587

return

588

538

589

def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):

539

def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):

590

540

591

minIndex = min(newheis[0])

541

#

592

maxIndex = max(newheis[0])

542

if horOnly:

543

A = numpy.c_[dir_cosu,dir_cosv]

544

else:

545

A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]

546

A = numpy.asmatrix(A)

547

A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()

593

548

594

voltage = voltage0[:,:,minIndex:maxIndex+1]

549

return A1

595

nLength = voltage.shape[1]/n

596

nMin = 0

597

nMax = 0

598

phaseOffset = numpy.zeros((len(pairslist),n))

599

550

600

for i in range(n):

551

def __correctValues(self, heiRang, phi, velRadial, SNR):

601

nMax += nLength

552

listPhi = phi.tolist()

602

phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))

553

maxid = listPhi.index(max(listPhi))

603

phaseCCF = numpy.mean(phaseCCF, axis = 2)

554

minid = listPhi.index(min(listPhi))

604

phaseOffset[:,i] = phaseCCF.transpose()

605

nMin = nMax

606

# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)

607

555

608

#Remove Outliers

556

rango = range(len(phi))

609

factor = 2

557

# rango = numpy.delete(rango,maxid)

610

wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))

611

dw = numpy.std(wt,axis = 1)

612

dw = dw.reshape((dw.size,1))

613

ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))

614

phaseOffset[ind] = numpy.nan

615

phaseOffset = stats.nanmean(phaseOffset, axis=1)

616

558

617

return phaseOffset

559

heiRang1 = heiRang*math.cos(phi[maxid])

560

heiRangAux = heiRang*math.cos(phi[minid])

561

indOut = (heiRang1 < heiRangAux[0]).nonzero()

562

heiRang1 = numpy.delete(heiRang1,indOut)

618

563

619

def __shiftPhase(self, data, phaseShift):

564

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

620

#this will shift the phase of a complex number

565

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

621

dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)

622

return dataShifted

623

566

624

def __estimatePhaseDifference(self, array, pairslist):

567

for i in rango:

625

nChannel = array.shape[0]

568

x = heiRang*math.cos(phi[i])

626

nHeights = array.shape[2]

569

y1 = velRadial[i,:]

627

numPairs = len(pairslist)

570

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

628

# phaseCCF = numpy.zeros((nChannel, 5, nHeights))

629

phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))

630

571

631

#Correct phases

572

x1 = heiRang1

632

derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]

573

y11 = f1(x1)

633

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

634

574

635

if indDer[0].shape[0] > 0:

575

y2 = SNR[i,:]

636

for i in range(indDer[0].shape[0]):

576

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

637

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])

577

y21 = f2(x1)

638

phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi

639

578

640

# for j in range(numSides):

579

velRadial1[i,:] = y11

641

# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])

580

SNR1[i,:] = y21

642

# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)

643

#

644

#Linear

645

phaseInt = numpy.zeros((numPairs,1))

646

angAllCCF = phaseCCF[:,[0,1,3,4],0]

647

for j in range(numPairs):

648

fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])

649

phaseInt[j] = fit[1]

650

#Phase Differences

651

phaseDiff = phaseInt - phaseCCF[:,2,:]

652

phaseArrival = phaseInt.reshape(phaseInt.size)

653

581

654

#Dealias

582

return heiRang1, velRadial1, SNR1

655

phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))

656

# indAlias = numpy.where(phaseArrival > numpy.pi)

657

# phaseArrival[indAlias] -= 2*numpy.pi

658

# indAlias = numpy.where(phaseArrival < -numpy.pi)

659

# phaseArrival[indAlias] += 2*numpy.pi

660

583

661

return phaseDiff, phaseArrival

584

def __calculateVelUVW(self, A, velRadial):

662

585

663

def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):

586

#Operacion Matricial

664

#this function will run the coherent detection used in Holdworth et al. 2004 and return the net power

587

# velUVW = numpy.zeros((velRadial.shape[1],3))

665

#find the phase shifts of each channel over 1 second intervals

588

# for ind in range(velRadial.shape[1]):

666

#only look at ranges below the beacon signal

589

# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])

667

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

590

# velUVW = velUVW.transpose()

668

numBlocks = int(volts.shape[1]/numProfPerBlock)

591

velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))

669

numHeights = volts.shape[2]

592

velUVW[:,:] = numpy.dot(A,velRadial)

670

nChannel = volts.shape[0]

671

voltsCohDet = volts.copy()

672

593

673

pairsarray = numpy.array(pairslist)

674

indSides = pairsarray[:,1]

675

# indSides = numpy.array(range(nChannel))

676

# indSides = numpy.delete(indSides, indCenter)

677

#

678

# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)

679

listBlocks = numpy.array_split(volts, numBlocks, 1)

680

594

681

startInd = 0

595

return velUVW

682

endInd = 0

683

596

684

for i in range(numBlocks):

597

# def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):

685

startInd = endInd

598

def techniqueDBS(self, velRadial0, heiRang, SNR0, kwargs):

686

endInd = endInd + listBlocks[i].shape[1]

599

"""

600

Function that implements Doppler Beam Swinging (DBS) technique.

687

601

688

arrayBlock = listBlocks[i]

602

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

689

# arrayBlockCenter = listCenter[i]

603

Direction correction (if necessary), Ranges and SNR

690

604

691

#Estimate the Phase Difference

605

Output: Winds estimation (Zonal, Meridional and Vertical)

692

phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)

693

#Phase Difference RMS

694

arrayPhaseRMS = numpy.abs(phaseDiff)

695

phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)

696

indPhase = numpy.where(phaseRMSaux==4)

697

#Shifting

698

if indPhase[0].shape[0] > 0:

699

for j in range(indSides.size):

700

arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())

701

voltsCohDet[:,startInd:endInd,:] = arrayBlock

702

606

703

return voltsCohDet

607

Parameters affected: Winds, height range, SNR

608

"""

704

609

705

def __calculateCCF(self, volts, pairslist ,laglist):

610

if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):

611

theta_x = numpy.array(kwargs['dirCosx'])

612

theta_y = numpy.array(kwargs['dirCosy'])

613

else:

614

elev = numpy.array(kwargs['elevation'])

615

azim = numpy.array(kwargs['azimuth'])

616

theta_x, theta_y = self.__calculateCosDir(elev, azim)

617

azimuth = kwargs['correctAzimuth']

618

if kwargs.has_key('horizontalOnly'):

619

horizontalOnly = kwargs['horizontalOnly']

620

else: horizontalOnly = False

621

if kwargs.has_key('correctFactor'):

622

correctFactor = kwargs['correctFactor']

623

else: correctFactor = 1

624

if kwargs.has_key('channelList'):

625

channelList = kwargs['channelList']

626

if len(channelList) == 2:

627

horizontalOnly = True

628

arrayChannel = numpy.array(channelList)

629

param = param[arrayChannel,:,:]

630

theta_x = theta_x[arrayChannel]

631

theta_y = theta_y[arrayChannel]

706

632

707

nHeights = volts.shape[2]

633

azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)

708

nPoints = volts.shape[1]

634

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)

709

voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')

635

A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)

710

636

711

for i in range(len(pairslist)):

637

#Calculo de Componentes de la velocidad con DBS

712

volts1 = volts[pairslist[i][0]]

638

winds = self.__calculateVelUVW(A,velRadial1)

713

volts2 = volts[pairslist[i][1]]

714

639

715

for t in range(len(laglist)):

640

return winds, heiRang1, SNR1

716

idxT = laglist[t]

717

if idxT >= 0:

718

vStacked = numpy.vstack((volts2[idxT:,:],

719

numpy.zeros((idxT, nHeights),dtype='complex')))

720

else:

721

vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),

722

volts2[:(nPoints + idxT),:]))

723

voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)

724

641

725

vStacked = None

642

def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):

726

return voltsCCF

727

643

728

def __getNoise(self, power, timeSegment, timeInterval):

644

posx = numpy.asarray(posx)

729

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

645

posy = numpy.asarray(posy)

730

numBlocks = int(power.shape[0]/numProfPerBlock)

731

numHeights = power.shape[1]

732

646

733

listPower = numpy.array_split(power, numBlocks, 0)

647

#Rotacion Inversa para alinear con el azimuth

734

noise = numpy.zeros((power.shape[0], power.shape[1]))

648

if azimuth!= None:

735

noise1 = numpy.zeros((power.shape[0], power.shape[1]))

649

azimuth = azimuth*math.pi/180

650

posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)

651

posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)

652

else:

653

posx1 = posx

654

posy1 = posy

736

655

737

startInd = 0

656

#Calculo de Distancias

738

endInd = 0

657

distx = numpy.zeros(pairsCrossCorr.size)

658

disty = numpy.zeros(pairsCrossCorr.size)

659

dist = numpy.zeros(pairsCrossCorr.size)

660

ang = numpy.zeros(pairsCrossCorr.size)

739

661

740

for i in range(numBlocks): #split por canal

662

for i in range(pairsCrossCorr.size):

741

startInd = endInd

663

distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]

742

endInd = endInd + listPower[i].shape[0]

664

disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]

665

dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)

666

ang[i] = numpy.arctan2(disty[i],distx[i])

667

#Calculo de Matrices

668

nPairs = len(pairs)

669

ang1 = numpy.zeros((nPairs, 2, 1))

670

dist1 = numpy.zeros((nPairs, 2, 1))

743

671

744

arrayBlock = listPower[i]

672

for j in range(nPairs):

745

noiseAux = numpy.mean(arrayBlock, 0)

673

dist1[j,0,0] = dist[pairs[j][0]]

746

# noiseAux = numpy.median(noiseAux)

674

dist1[j,1,0] = dist[pairs[j][1]]

747

# noiseAux = numpy.mean(arrayBlock)

675

ang1[j,0,0] = ang[pairs[j][0]]

748

noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux

676

ang1[j,1,0] = ang[pairs[j][1]]

749

677

750

noiseAux1 = numpy.mean(arrayBlock)

678

return distx,disty, dist1,ang1

751

noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1

752

679

753

return noise, noise1

680

def __calculateVelVer(self, phase, lagTRange, _lambda):

754

681

755

def __findMeteors(self, power, thresh):

682

Ts = lagTRange[1] - lagTRange[0]

756

nProf = power.shape[0]

683

velW = -_lambda*phase/(4*math.pi*Ts)

757

nHeights = power.shape[1]

758

listMeteors = []

759

684

760

for i in range(nHeights):

685

return velW

761

powerAux = power[:,i]

762

threshAux = thresh[:,i]

763

686

764

indUPthresh = numpy.where(powerAux > threshAux)[0]

687

def __calculateVelHorDir(self, dist, tau1, tau2, ang):

765

indDNthresh = numpy.where(powerAux <= threshAux)[0]

688

nPairs = tau1.shape[0]

689

vel = numpy.zeros((nPairs,3,tau1.shape[2]))

766

690

767

j = 0

691

angCos = numpy.cos(ang)

692

angSin = numpy.sin(ang)

768

693

769

while (j < indUPthresh.size - 2):

694

vel0 = dist*tau1/(2*tau2**2)

770

if (indUPthresh[j + 2] == indUPthresh[j] + 2):

695

vel[:,0,:] = (vel0*angCos).sum(axis = 1)

771

indDNAux = numpy.where(indDNthresh > indUPthresh[j])

696

vel[:,1,:] = (vel0*angSin).sum(axis = 1)

772

indDNthresh = indDNthresh[indDNAux]

773

697

774

if (indDNthresh.size > 0):

698

ind = numpy.where(numpy.isinf(vel))

775

indEnd = indDNthresh[0] - 1

699

vel[ind] = numpy.nan

776

indInit = indUPthresh[j]

777

700

778

meteor = powerAux[indInit:indEnd + 1]

701

return vel

779

indPeak = meteor.argmax() + indInit

780

FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))

781

702

782

listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!

703

def __getPairsAutoCorr(self, pairsList, nChannels):

783

j = numpy.where(indUPthresh == indEnd)[0] + 1

784

else: j+=1

785

else: j+=1

786

704

787

return listMeteors

705

pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

788

706

789

def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):

707

for l in range(len(pairsList)):

708

firstChannel = pairsList[l][0]

709

secondChannel = pairsList[l][1]

790

710

791

arrayMeteors = numpy.asarray(listMeteors)

711

#Obteniendo pares de Autocorrelacion

792

listMeteors1 = []

712

if firstChannel == secondChannel:

713

pairsAutoCorr[firstChannel] = int(l)

793

714

794

while arrayMeteors.shape[0] > 0:

715

pairsAutoCorr = pairsAutoCorr.astype(int)

795

FLAs = arrayMeteors[:,4]

796

maxFLA = FLAs.argmax()

797

listMeteors1.append(arrayMeteors[maxFLA,:])

798

716

799

MeteorInitTime = arrayMeteors[maxFLA,1]

717

pairsCrossCorr = range(len(pairsList))

800

MeteorEndTime = arrayMeteors[maxFLA,3]

718

pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

801

MeteorHeight = arrayMeteors[maxFLA,0]

802

719

803

#Check neighborhood

720

return pairsAutoCorr, pairsCrossCorr

804

maxHeightIndex = MeteorHeight + rangeLimit

805

minHeightIndex = MeteorHeight - rangeLimit

806

minTimeIndex = MeteorInitTime - timeLimit

807

maxTimeIndex = MeteorEndTime + timeLimit

808

721

809

#Check Heights

722

def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):

810

indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)

723

"""

811

indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)

724

Function that implements Spaced Antenna (SA) technique.

812

indBoth = numpy.where(numpy.logical_and(indTime,indHeight))

813

725

814

arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)

726

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

727

Direction correction (if necessary), Ranges and SNR

815

728

816

return listMeteors1

729

Output: Winds estimation (Zonal, Meridional and Vertical)

817

730

818

def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):

731

Parameters affected: Winds

819

numHeights = volts.shape[2]

732

"""

820

nChannel = volts.shape[0]

733

#Cross Correlation pairs obtained

734

pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

735

pairsArray = numpy.array(pairsList)[pairsCrossCorr]

736

pairsSelArray = numpy.array(pairsSelected)

737

pairs = []

821

738

822

thresholdPhase = thresh[0]

739

#Wind estimation pairs obtained

823

thresholdNoise = thresh[1]

740

for i in range(pairsSelArray.shape[0]/2):

824

thresholdDB = float(thresh[2])

741

ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]

742

ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]

743

pairs.append((ind1,ind2))

825

744

826

thresholdDB1 = 10**(thresholdDB/10)

745

indtau = tau.shape[0]/2

827

pairsarray = numpy.array(pairslist)

746

tau1 = tau[:indtau,:]

828

indSides = pairsarray[:,1]

747

tau2 = tau[indtau:-1,:]

748

tau1 = tau1[pairs,:]

749

tau2 = tau2[pairs,:]

750

phase1 = tau[-1,:]

829

751

830

pairslist1 = list(pairslist)

752

#---------------------------------------------------------------------

831

pairslist1.append((0,1))

753

#Metodo Directo

832

pairslist1.append((3,4))

754

distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)

755

winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)

756

winds = stats.nanmean(winds, axis=0)

757

#---------------------------------------------------------------------

758

#Metodo General

759

# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)

760

# #Calculo Coeficientes de Funcion de Correlacion

761

# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)

762

# #Calculo de Velocidades

763

# winds = self.calculateVelUV(F,G,A,B,H)

833

764

834

listMeteors1 = []

765

#---------------------------------------------------------------------

835

listPowerSeries = []

766

winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)

836

listVoltageSeries = []

767

winds = correctFactor*winds

837

#volts has the war data

768

return winds

838

769

839

if frequency == 30e6:

770

def __checkTime(self, currentTime, paramInterval, outputInterval):

840

timeLag = 45*10**-3

841

else:

842

timeLag = 15*10**-3

843

lag = numpy.ceil(timeLag/timeInterval)

844

771

845

for i in range(len(listMeteors)):

772

dataTime = currentTime + paramInterval

773

deltaTime = dataTime - self.__initime

846

774

847

###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################

775

if deltaTime >= outputInterval or deltaTime < 0:

848

meteorAux = numpy.zeros(16)

776

self.__dataReady = True

777

return

849

778

850

#Loading meteor Data (mHeight, mStart, mPeak, mEnd)

779

def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):

851

mHeight = listMeteors[i][0]

780

'''

852

mStart = listMeteors[i][1]

781

Function that implements winds estimation technique with detected meteors.

853

mPeak = listMeteors[i][2]

854

mEnd = listMeteors[i][3]

855

782

856

#get the volt data between the start and end times of the meteor

783

Input: Detected meteors, Minimum meteor quantity to wind estimation

857

meteorVolts = volts[:,mStart:mEnd+1,mHeight]

858

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

859

784

860

#3.6. Phase Difference estimation

785

Output: Winds estimation (Zonal and Meridional)

861

phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)

862

863

#3.7. Phase difference removal & meteor start, peak and end times reestimated

864

#meteorVolts0.- all Channels, all Profiles

865

meteorVolts0 = volts[:,:,mHeight]

866

meteorThresh = noise[:,mHeight]*thresholdNoise

867

meteorNoise = noise[:,mHeight]

868

meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting

869

powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power

870

871

#Times reestimation

872

mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]

873

if mStart1.size > 0:

874

mStart1 = mStart1[-1] + 1

875

876

else:

877

mStart1 = mPeak

878

786

879

mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1

787

Parameters affected: Winds

880

mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]

788

'''

881

if mEndDecayTime1.size == 0:

789

# print arrayMeteor.shape

882

mEndDecayTime1 = powerNet0.size

790

#Settings

883

else:

791

nInt = (heightMax - heightMin)/2

884

mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1

792

# print nInt

885

# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()

793

nInt = int(nInt)

794

# print nInt

795

winds = numpy.zeros((2,nInt))*numpy.nan

886

796

887

#meteorVolts1.- all Channels, from start to end

797

#Filter errors

888

meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]

798

error = numpy.where(arrayMeteor[:,-1] == 0)[0]

889

meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]

799

finalMeteor = arrayMeteor[error,:]

890

if meteorVolts2.shape[1] == 0:

891

meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]

892

meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)

893

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)

894

##################### END PARAMETERS REESTIMATION #########################

895

800

896

##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################

801

#Meteor Histogram

897

# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis

802

finalHeights = finalMeteor[:,2]

898

if meteorVolts2.shape[1] > 0:

803

hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))

899

#Phase Difference re-estimation

804

nMeteorsPerI = hist[0]

900

phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation

805

heightPerI = hist[1]

901

# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)

902

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])

903

phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))

904

meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting

905

806

906

#Phase Difference RMS

807

#Sort of meteors

907

phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))

808

indSort = finalHeights.argsort()

908

powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)

809

finalMeteor2 = finalMeteor[indSort,:]

909

#Data from Meteor

910

mPeak1 = powerNet1.argmax() + mStart1

911

mPeakPower1 = powerNet1.max()

912

noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])

913

mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux

914

Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])

915

Meteor1 = numpy.hstack((Meteor1,phaseDiffint))

916

PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]

917

#Vectorize

918

meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]

919

meteorAux[7:11] = phaseDiffint[0:4]

920

810

921

#Rejection Criterions

811

# Calculating winds

922

if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation

812

ind1 = 0

923

meteorAux[-1] = 17

813

ind2 = 0

924

elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB

925

meteorAux[-1] = 1

926

814

815

for i in range(nInt):

816

nMet = nMeteorsPerI[i]

817

ind1 = ind2

818

ind2 = ind1 + nMet

927

819

928

else:

820

meteorAux = finalMeteor2[ind1:ind2,:]

929

meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]

930

meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis

931

PowerSeries = 0

932

821

933

listMeteors1.append(meteorAux)

822

if meteorAux.shape[0] >= meteorThresh:

934

listPowerSeries.append(PowerSeries)

823

vel = meteorAux[:, 6]

935

listVoltageSeries.append(meteorVolts1)

824

zen = meteorAux[:, 4]*numpy.pi/180

825

azim = meteorAux[:, 3]*numpy.pi/180

936

826

937

return listMeteors1, listPowerSeries, listVoltageSeries

827

n = numpy.cos(zen)

828

# m = (1 - n**2)/(1 - numpy.tan(azim)**2)

829

# l = m*numpy.tan(azim)

830

l = numpy.sin(zen)*numpy.sin(azim)

831

m = numpy.sin(zen)*numpy.cos(azim)

938

832

939

def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):

833

A = numpy.vstack((l, m)).transpose()

834

A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())

835

windsAux = numpy.dot(A1, vel)

940

836

941

threshError = 10

837

winds[0,i] = windsAux[0]

942

#Depending if it is 30 or 50 MHz

838

winds[1,i] = windsAux[1]

943

if frequency == 30e6:

944

timeLag = 45*10**-3

945

else:

946

timeLag = 15*10**-3

947

lag = numpy.ceil(timeLag/timeInterval)

948

839

949

listMeteors1 = []

840

return winds, heightPerI[:-1]

950

841

951

for i in range(len(listMeteors)):

842

def techniqueNSM_SA(self, **kwargs):

952

meteorPower = listPower[i]

843

metArray = kwargs['metArray']

953

meteorAux = listMeteors[i]

844

heightList = kwargs['heightList']

845

timeList = kwargs['timeList']

954

846

955

if meteorAux[-1] == 0:

847

rx_location = kwargs['rx_location']

848

groupList = kwargs['groupList']

849

azimuth = kwargs['azimuth']

850

dfactor = kwargs['dfactor']

851

k = kwargs['k']

956

852

957

try:

853

azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)

958

indmax = meteorPower.argmax()

854

d = dist*dfactor

959

indlag = indmax + lag

855

#Phase calculation

856

metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)

960

857

961

y = meteorPower[indlag:]

858

metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities

962

x = numpy.arange(0, y.size)*timeLag

963

859

964

#first guess

860

velEst = numpy.zeros((heightList.size,2))*numpy.nan

965

a = y[0]

861

azimuth1 = azimuth1*numpy.pi/180

966

tau = timeLag

967

#exponential fit

968

popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])

969

y1 = self.__exponential_function(x, *popt)

970

#error estimation

971

error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))

972

862

973

decayTime = popt[1]

863

for i in range(heightList.size):

974

riseTime = indmax*timeInterval

864

h = heightList[i]

975

meteorAux[11:13] = [decayTime, error]

865

indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]

866

metHeight = metArray1[indH,:]

867

if metHeight.shape[0] >= 2:

868

velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities

869

iazim = metHeight[:,1].astype(int)

870

azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths

871

A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))

872

A = numpy.asmatrix(A)

873

A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()

874

velHor = numpy.dot(A1,velAux)

976

875

977

#Table items 7, 8 and 11

876

velEst[i,:] = numpy.squeeze(velHor)

978

if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s

877

return velEst

979

meteorAux[-1] = 7

980

elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time

981

meteorAux[-1] = 8

982

if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time

983

meteorAux[-1] = 11

984

878

879

def __getPhaseSlope(self, metArray, heightList, timeList):

880

meteorList = []

881

#utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2

882

#Putting back together the meteor matrix

883

utctime = metArray[:,0]

884

uniqueTime = numpy.unique(utctime)

985

885

986

except:

886

phaseDerThresh = 0.5

987

meteorAux[-1] = 11

887

ippSeconds = timeList[1] - timeList[0]

888

sec = numpy.where(timeList>1)[0][0]

889

nPairs = metArray.shape[1] - 6

890

nHeights = len(heightList)

988

891

892

for t in uniqueTime:

893

metArray1 = metArray[utctime==t,:]

894

# phaseDerThresh = numpy.pi/4 #reducir Phase thresh

895

tmet = metArray1[:,1].astype(int)

896

hmet = metArray1[:,2].astype(int)

989

897

990

listMeteors1.append(meteorAux)

898

metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))

899

metPhase[:,:] = numpy.nan

900

metPhase[:,hmet,tmet] = metArray1[:,6:].T

991

901

992

return listMeteors1

902

#Delete short trails

903

metBool = ~numpy.isnan(metPhase[0,:,:])

904

heightVect = numpy.sum(metBool, axis = 1)

905

metBool[heightVect<sec,:] = False

906

metPhase[:,heightVect<sec,:] = numpy.nan

993

907

994

#Exponential Function

908

#Derivative

909

metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])

910

phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))

911

metPhase[phDerAux] = numpy.nan

995

912

996

def __exponential_function(self, x, a, tau):

913

#--------------------------METEOR DETECTION -----------------------------------------

997

y = a*numpy.exp(-x/tau)

914

indMet = numpy.where(numpy.any(metBool,axis=1))[0]

998

return y

999

915

1000

def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):

916

for p in numpy.arange(nPairs):

917

phase = metPhase[p,:,:]

918

phDer = metDer[p,:,:]

1001

919

1002

pairslist1 = list(pairslist)

920

for h in indMet:

1003

pairslist1.append((0,1))

921

height = heightList[h]

1004

pairslist1.append((3,4))

922

phase1 = phase[h,:] #82

1005

numPairs = len(pairslist1)

923

phDer1 = phDer[h,:]

1006

#Time Lag

1007

timeLag = 45*10**-3

1008

c = 3e8

1009

lag = numpy.ceil(timeLag/timeInterval)

1010

freq = 30e6

1011

924

1012

listMeteors1 = []

925

phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap

1013

926

1014

for i in range(len(listMeteors)):

927

indValid = numpy.where(~numpy.isnan(phase1))[0]

1015

meteorAux = listMeteors[i]

928

initMet = indValid[0]

1016

if meteorAux[-1] == 0:

929

endMet = 0

1017

mStart = listMeteors[i][1]

1018

mPeak = listMeteors[i][2]

1019

mLag = mPeak - mStart + lag

1020

930

1021

#get the volt data between the start and end times of the meteor

931

for i in range(len(indValid)-1):

1022

meteorVolts = listVolts[i]

1023

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

1024

932

1025

#Get CCF

933

#Time difference

1026

allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])

934

inow = indValid[i]

935

inext = indValid[i+1]

936

idiff = inext - inow

937

#Phase difference

938

phDiff = numpy.abs(phase1[inext] - phase1[inow])

1027

939

1028

#Method 2

940

if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor

1029

slopes = numpy.zeros(numPairs)

941

sizeTrail = inow - initMet + 1

1030

time = numpy.array([-2,-1,1,2])*timeInterval

942

if sizeTrail>3*sec: #Too short meteors

1031

angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])

943

x = numpy.arange(initMet,inow+1)*ippSeconds

944

y = phase1[initMet:inow+1]

945

ynnan = ~numpy.isnan(y)

946

x = x[ynnan]

947

y = y[ynnan]

948

slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)

949

ylin = x*slope + intercept

950

rsq = r_value**2

951

if rsq > 0.5:

952

vel = slope#*height*1000/(k*d)

953

estAux = numpy.array([utctime,p,height, vel, rsq])

954

meteorList.append(estAux)

955

initMet = inext

956

metArray2 = numpy.array(meteorList)

1032

957

1033

#Correct phases

958

return metArray2

1034

derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]

1035

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

1036

959

1037

if indDer[0].shape[0] > 0:

960

def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):

1038

for i in range(indDer[0].shape[0]):

1039

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])

1040

angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi

1041

961

1042

# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))

962

azimuth1 = numpy.zeros(len(pairslist))

1043

for j in range(numPairs):

963

dist = numpy.zeros(len(pairslist))

1044

fit = stats.linregress(time, angAllCCF[j,:])

1045

slopes[j] = fit[0]

1046

964

1047

#Remove Outlier

965

for i in range(len(rx_location)):

1048

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

966

ch0 = pairslist[i][0]

1049

# slopes = numpy.delete(slopes,indOut)

967

ch1 = pairslist[i][1]

1050

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

1051

# slopes = numpy.delete(slopes,indOut)

1052

968

1053

radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)

969

diffX = rx_location[ch0][0] - rx_location[ch1][0]

1054

radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)

970

diffY = rx_location[ch0][1] - rx_location[ch1][1]

1055

meteorAux[-2] = radialError

971

azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi

1056

meteorAux[-3] = radialVelocity

972

dist[i] = numpy.sqrt(diffX**2 + diffY**2)

1057

973

1058

#Setting Error

974

azimuth1 -= azimuth0

1059

#Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s

975

return azimuth1, dist

1060

if numpy.abs(radialVelocity) > 200:

1061

meteorAux[-1] = 15

1062

#Number 12: Poor fit to CCF variation for estimation of radial drift velocity

1063

elif radialError > radialStdThresh:

1064

meteorAux[-1] = 12

1065

976

1066

listMeteors1.append(meteorAux)

977

def techniqueNSM_DBS(self, **kwargs):

1067

return listMeteors1

978

metArray = kwargs['metArray']

979

heightList = kwargs['heightList']

980

timeList = kwargs['timeList']

981

zenithList = kwargs['zenithList']

982

nChan = numpy.max(cmet) + 1

983

nHeights = len(heightList)

1068

984

1069

def __setNewArrays(self, listMeteors, date, heiRang):

985

utctime = metArray[:,0]

986

cmet = metArray[:,1]

987

hmet = metArray1[:,3].astype(int)

988

h1met = heightList[hmet]*zenithList[cmet]

989

vmet = metArray1[:,5]

1070

990

1071

#New arrays

991

for i in range(nHeights - 1):

1072

arrayMeteors = numpy.array(listMeteors)

992

hmin = heightList[i]

1073

arrayParameters = numpy.zeros((len(listMeteors), 13))

993

hmax = heightList[i + 1]

1074

994

1075

#Date inclusion

995

vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]

1076

# date = re.findall(r'\((.*?)\)', date)

1077

# date = date[0].split(',')

1078

# date = map(int, date)

1079

#

1080

# if len(date)<6:

1081

# date.append(0)

1082

#

1083

# date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]

1084

# arrayDate = numpy.tile(date, (len(listMeteors), 1))

1085

arrayDate = numpy.tile(date, (len(listMeteors)))

1086

996

1087

#Meteor array

1088

# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]

1089

# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))

1090

997

1091

#Parameters Array

1092

arrayParameters[:,0] = arrayDate #Date

1093

arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range

1094

arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error

1095

arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases

1096

arrayParameters[:,-1] = arrayMeteors[:,-1] #Error

1097

998

999

return data_output

1098

1000

1099

return arrayParameters

1001

def run(self, dataOut, technique, **kwargs):

1100

1002

1101

# def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):

1003

param = dataOut.data_param

1102

#

1004

if dataOut.abscissaList != None:

1103

# arrayAOA = numpy.zeros((phases.shape[0],3))

1005

absc = dataOut.abscissaList[:-1]

1104

# cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)

1006

noise = dataOut.noise

1105

#

1007

heightList = dataOut.heightList

1106

# arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)

1008

SNR = dataOut.data_SNR

1107

# cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)

1108

# arrayAOA[:,2] = cosDirError

1109

#

1110

# azimuthAngle = arrayAOA[:,0]

1111

# zenithAngle = arrayAOA[:,1]

1112

#

1113

# #Setting Error

1114

# #Number 3: AOA not fesible

1115

# indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]

1116

# error[indInvalid] = 3

1117

# #Number 4: Large difference in AOAs obtained from different antenna baselines

1118

# indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]

1119

# error[indInvalid] = 4

1120

# return arrayAOA, error

1121

#

1122

# def __getDirectionCosines(self, arrayPhase, pairsList):

1123

#

1124

# #Initializing some variables

1125

# ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi

1126

# ang_aux = ang_aux.reshape(1,ang_aux.size)

1127

#

1128

# cosdir = numpy.zeros((arrayPhase.shape[0],2))

1129

# cosdir0 = numpy.zeros((arrayPhase.shape[0],2))

1130

#

1131

#

1132

# for i in range(2):

1133

# #First Estimation

1134

# phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]

1135

# #Dealias

1136

# indcsi = numpy.where(phi0_aux > numpy.pi)

1137

# phi0_aux[indcsi] -= 2*numpy.pi

1138

# indcsi = numpy.where(phi0_aux < -numpy.pi)

1139

# phi0_aux[indcsi] += 2*numpy.pi

1140

# #Direction Cosine 0

1141

# cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)

1142

#

1143

# #Most-Accurate Second Estimation

1144

# phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]

1145

# phi1_aux = phi1_aux.reshape(phi1_aux.size,1)

1146

# #Direction Cosine 1

1147

# cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)

1148

#

1149

# #Searching the correct Direction Cosine

1150

# cosdir0_aux = cosdir0[:,i]

1151

# cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)

1152

# #Minimum Distance

1153

# cosDiff = (cosdir1 - cosdir0_aux)**2

1154

# indcos = cosDiff.argmin(axis = 1)

1155

# #Saving Value obtained

1156

# cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]

1157

#

1158

# return cosdir0, cosdir

1159

#

1160

# def __calculateAOA(self, cosdir, azimuth):

1161

# cosdirX = cosdir[:,0]

1162

# cosdirY = cosdir[:,1]

1163

#

1164

# zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi

1165

# azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east

1166

# angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()

1167

#

1168

# return angles

1169

#

1170

# def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):

1171

#

1172

# Ramb = 375 #Ramb = c/(2*PRF)

1173

# Re = 6371 #Earth Radius

1174

# heights = numpy.zeros(Ranges.shape)

1175

#

1176

# R_aux = numpy.array([0,1,2])*Ramb

1177

# R_aux = R_aux.reshape(1,R_aux.size)

1178

#

1179

# Ranges = Ranges.reshape(Ranges.size,1)

1180

#

1181

# Ri = Ranges + R_aux

1182

# hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re

1183

#

1184

# #Check if there is a height between 70 and 110 km

1185

# h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)

1186

# ind_h = numpy.where(h_bool == 1)[0]

1187

#

1188

# hCorr = hi[ind_h, :]

1189

# ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))

1190

#

1191

# hCorr = hi[ind_hCorr]

1192

# heights[ind_h] = hCorr

1193

#

1194

# #Setting Error

1195

# #Number 13: Height unresolvable echo: not valid height within 70 to 110 km

1196

# #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km

1197

#

1198

# indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]

1199

# error[indInvalid2] = 14

1200

# indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]

1201

# error[indInvalid1] = 13

1202

#

1203

# return heights, error

1204

1009

1205

def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):

1010

if technique == 'DBS':

1011

# if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):

1012

# theta_x = numpy.array(kwargs['dirCosx'])

1013

# theta_y = numpy.array(kwargs['dirCosy'])

1014

# else:

1015

# elev = numpy.array(kwargs['elevation'])

1016

# azim = numpy.array(kwargs['azimuth'])

1017

# theta_x, theta_y = self.__calculateCosDir(elev, azim)

1018

# azimuth = kwargs['correctAzimuth']

1019

# if kwargs.has_key('horizontalOnly'):

1020

# horizontalOnly = kwargs['horizontalOnly']

1021

# else: horizontalOnly = False

1022

# if kwargs.has_key('correctFactor'):

1023

# correctFactor = kwargs['correctFactor']

1024

# else: correctFactor = 1

1025

# if kwargs.has_key('channelList'):

1026

# channelList = kwargs['channelList']

1027

# if len(channelList) == 2:

1028

# horizontalOnly = True

1029

# arrayChannel = numpy.array(channelList)

1030

# param = param[arrayChannel,:,:]

1031

# theta_x = theta_x[arrayChannel]

1032

# theta_y = theta_y[arrayChannel]

1206

1033

1207

'''

1034

velRadial0 = param[:,1,:] #Radial velocity

1208

Function GetMoments()

1035

# dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function

1036

dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, heightList, SNR, kwargs) #DBS Function

1037

dataOut.utctimeInit = dataOut.utctime

1038

dataOut.outputInterval = dataOut.paramInterval

1209

1039

1210

Input:

1040

elif technique == 'SA':

1211

Output:

1212

Variables modified:

1213

'''

1214

if path != None:

1215

sys.path.append(path)

1216

self.dataOut.library = importlib.import_module(file)

1217

1041

1218

#To be inserted as a parameter

1042

#Parameters

1219

groupArray = numpy.array(groupList)

1043

position_x = kwargs['positionX']

1220

# groupArray = numpy.array([[0,1],[2,3]])

1044

position_y = kwargs['positionY']

1221

self.dataOut.groupList = groupArray

1045

azimuth = kwargs['azimuth']

1222

1046

1223

nGroups = groupArray.shape[0]

1047

if kwargs.has_key('crosspairsList'):

1224

nChannels = self.dataIn.nChannels

1048

pairs = kwargs['crosspairsList']

1225

nHeights=self.dataIn.heightList.size

1049

else:

1050

pairs = None

1226

1051

1227

#Parameters Array

1052

if kwargs.has_key('correctFactor'):

1228

self.dataOut.data_param = None

1053

correctFactor = kwargs['correctFactor']

1054

else:

1055

correctFactor = 1

1229

1056

1230

#Set constants

1057

tau = dataOut.data_param

1231

constants = self.dataOut.library.setConstants(self.dataIn)

1058

_lambda = dataOut.C/dataOut.frequency

1232

self.dataOut.constants = constants

1059

pairsList = dataOut.groupList

1233

M = self.dataIn.normFactor

1060

nChannels = dataOut.nChannels

1234

N = self.dataIn.nFFTPoints

1235

ippSeconds = self.dataIn.ippSeconds

1236

K = self.dataIn.nIncohInt

1237

pairsArray = numpy.array(self.dataIn.pairsList)

1238

1061

1239

#List of possible combinations

1062

dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)

1240

listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)

1063

dataOut.utctimeInit = dataOut.utctime

1241

indCross = numpy.zeros(len(list(listComb)), dtype = 'int')

1064

dataOut.outputInterval = dataOut.timeInterval

1242

1065

1243

if getSNR:

1066

elif technique == 'Meteors':

1244

listChannels = groupArray.reshape((groupArray.size))

1067

dataOut.flagNoData = True

1245

listChannels.sort()

1068

self.__dataReady = False

1246

noise = self.dataIn.getNoise()

1247

self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])

1248

1069

1249

for i in range(nGroups):

1070

if kwargs.has_key('nHours'):

1250

coord = groupArray[i,:]

1071

nHours = kwargs['nHours']

1072

else:

1073

nHours = 1

1251

1074

1252

#Input data array

1075

if kwargs.has_key('meteorsPerBin'):

1253

data = self.dataIn.data_spc[coord,:,:]/(M*N)

1076

meteorThresh = kwargs['meteorsPerBin']

1254

data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))

1077

else:

1078

meteorThresh = 6

1255

1079

1256

#Cross Spectra data array for Covariance Matrixes

1080

if kwargs.has_key('hmin'):

1257

ind = 0

1081

hmin = kwargs['hmin']

1258

for pairs in listComb:

1082

else: hmin = 70

1259

pairsSel = numpy.array([coord[x],coord[y]])

1083

if kwargs.has_key('hmax'):

1260

indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])

1084

hmax = kwargs['hmax']

1261

~~ind~~ += 1

1085

else: hmax = 110

1262

dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)

1263

dataCross = dataCross**2/K

1264

1086

1265

for h in range(nHeights):

1087

dataOut.outputInterval = nHours*3600

1266

# print self.dataOut.heightList[h]

1267

1088

1268

#Input

1089

if self.__isConfig == False:

1269

d = data[:,h]

1090

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

1091

#Get Initial LTC time

1092

self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

1093

self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

1270

1094

1271

#Covariance Matrix

1095

self.__isConfig = True

1272

D = numpy.diag(d**2/K)

1273

ind = 0

1274

for pairs in listComb:

1275

#Coordinates in Covariance Matrix

1276

x = pairs[0]

1277

y = pairs[1]

1278

#Channel Index

1279

S12 = dataCross[ind,:,h]

1280

D12 = numpy.diag(S12)

1281

#Completing Covariance Matrix with Cross Spectras

1282

D[x*N:(x+1)*N,y*N:(y+1)*N] = D12

1283

D[y*N:(y+1)*N,x*N:(x+1)*N] = D12

1284

ind += 1

1285

Dinv=numpy.linalg.inv(D)

1286

L=numpy.linalg.cholesky(Dinv)

1287

LT=L.T

1288

1096

1289

dp = numpy.dot(LT,d)

1097

if self.__buffer == None:

1098

self.__buffer = dataOut.data_param

1099

self.__firstdata = copy.copy(dataOut)

1290

1100

1291

#Initial values

1101

else:

1292

data_spc = self.dataIn.data_spc[coord,:,h]

1102

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

1103

1104

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

1105

1106

if self.__dataReady:

1107

dataOut.utctimeInit = self.__initime

1108

1109

self.__initime += dataOut.outputInterval #to erase time offset

1110

1111

dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)

1112

dataOut.flagNoData = False

1113

self.__buffer = None

1114

1115

elif technique == 'Meteors1':

1116

dataOut.flagNoData = True

1117

self.__dataReady = False

1118

1119

if kwargs.has_key('nMins'):

1120

nMins = kwargs['nMins']

1121

else: nMins = 20

1122

if kwargs.has_key('rx_location'):

1123

rx_location = kwargs['rx_location']

1124

else: rx_location = [(0,1),(1,1),(1,0)]

1125

if kwargs.has_key('azimuth'):

1126

azimuth = kwargs['azimuth']

1127

else: azimuth = 51

1128

if kwargs.has_key('dfactor'):

1129

dfactor = kwargs['dfactor']

1130

if kwargs.has_key('mode'):

1131

mode = kwargs['mode']

1132

else: mode = 'SA'

1133

1134

#Borrar luego esto

1135

if dataOut.groupList == None:

1136

dataOut.groupList = [(0,1),(0,2),(1,2)]

1137

groupList = dataOut.groupList

1138

C = 3e8

1139

freq = 50e6

1140

lamb = C/freq

1141

k = 2*numpy.pi/lamb

1142

1143

timeList = dataOut.abscissaList

1144

heightList = dataOut.heightList

1145

1146

if self.__isConfig == False:

1147

dataOut.outputInterval = nMins*60

1148

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

1149

#Get Initial LTC time

1150

initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

1151

minuteAux = initime.minute

1152

minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)

1153

self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

1154

1155

self.__isConfig = True

1156

1157

if self.__buffer == None:

1158

self.__buffer = dataOut.data_param

1159

self.__firstdata = copy.copy(dataOut)

1293

1160

1294

if (h>0)and(error1[3]<5):

1295

p0 = self.dataOut.data_param[i,:,h-1]

1296

else:

1161

else:

1297

p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))

1162

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

1298

1163

1299

try:

1164

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

1300

#Least Squares

1301

minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)

1302

# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))

1303

#Chi square error

1304

error0 = numpy.sum(infodict['fvec']**2)/(2*N)

1305

#Error with Jacobian

1306

error1 = self.dataOut.library.errorFunction(minp,constants,LT)

1307

except:

1308

minp = p0*numpy.nan

1309

error0 = numpy.nan

1310

error1 = p0*numpy.nan

1311

1165

1312

#Save

1166

if self.__dataReady:

1313

if self.~~dataOut~~.~~data_param~~ == ~~None~~:

1167

dataOut.utctimeInit = self.__initime

1314

self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan

1168

self.__initime += dataOut.outputInterval #to erase time offset

1315

self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan

1169

1170

metArray = self.__buffer

1171

if mode == 'SA':

1172

dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)

1173

elif mode == 'DBS':

1174

dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)

1175

dataOut.data_output = dataOut.data_output.T

1176

dataOut.flagNoData = False

1177

self.__buffer = None

1316

1178

1317

self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))

1318

self.dataOut.data_param[i,:,h] = minp

1319

return

1179

return

1320

1180

1321

def __residFunction(self, p, dp, LT, constants):

1181

class EWDriftsEstimation(Operation):

1322

1182

1323

fm = self.dataOut.library.modelFunction(p, constants)

1183

def __init__(self):

1324

fmp=numpy.dot(LT,fm)

1184

Operation.__init__(self)

1325

1185

1326

return dp-fmp

1186

def __correctValues(self, heiRang, phi, velRadial, SNR):

1187

listPhi = phi.tolist()

1188

maxid = listPhi.index(max(listPhi))

1189

minid = listPhi.index(min(listPhi))

1190

1191

rango = range(len(phi))

1192

# rango = numpy.delete(rango,maxid)

1193

1194

heiRang1 = heiRang*math.cos(phi[maxid])

1195

heiRangAux = heiRang*math.cos(phi[minid])

1196

indOut = (heiRang1 < heiRangAux[0]).nonzero()

1197

heiRang1 = numpy.delete(heiRang1,indOut)

1198

1199

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

1200

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

1201

1202

for i in rango:

1203

x = heiRang*math.cos(phi[i])

1204

y1 = velRadial[i,:]

1205

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

1206

1207

x1 = heiRang1

1208

y11 = f1(x1)

1209

1210

y2 = SNR[i,:]

1211

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

1212

y21 = f2(x1)

1213

1214

velRadial1[i,:] = y11

1215

SNR1[i,:] = y21

1216

1217

return heiRang1, velRadial1, SNR1

1218

1219

def run(self, dataOut, zenith, zenithCorrection):

1220

heiRang = dataOut.heightList

1221

velRadial = dataOut.data_param[:,3,:]

1222

SNR = dataOut.data_SNR

1223

1224

zenith = numpy.array(zenith)

1225

zenith -= zenithCorrection

1226

zenith *= numpy.pi/180

1227

1228

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)

1229

1230

alp = zenith[0]

1231

bet = zenith[1]

1232

1233

w_w = velRadial1[0,:]

1234

w_e = velRadial1[1,:]

1235

1236

w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))

1237

u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))

1238

1239

winds = numpy.vstack((u,w))

1240

1241

dataOut.heightList = heiRang1

1242

dataOut.data_output = winds

1243

dataOut.data_SNR = SNR1

1327

1244

1328

def __getSNR(self, z, noise):

1245

dataOut.utctimeInit = dataOut.utctime

1246

dataOut.outputInterval = dataOut.timeInterval

1247

return

1329

1248

1330

avg = numpy.average(z, axis=1)

1249

#--------------- Non Specular Meteor ----------------

1331

SNR = (avg.T-noise)/noise

1332

SNR = SNR.T

1333

return SNR

1334

1250

1335

def __chisq(p,chindex,hindex):

1251

class NonSpecularMeteorDetection(Operation):

1336

#similar to Resid but calculates CHI**2

1337

[LT,d,fm]=setupLTdfm(p,chindex,hindex)

1338

dp=numpy.dot(LT,d)

1339

fmp=numpy.dot(LT,fm)

1340

chisq=numpy.dot((dp-fmp).T,(dp-fmp))

1341

return chisq

1342

1252

1343

def ~~NonSpecularMeteorDetectio~~n(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):

1253

def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):

1344

data_acf = self.dataOut.data_pre[0]

1254

data_acf = self.dataOut.data_pre[0]

1345

data_ccf = self.dataOut.data_pre[1]

1255

data_ccf = self.dataOut.data_pre[1]

1346

1256

1523

1433

1524

return binArray1

1434

return binArray1

1525

1435

1526

def WeirdEcho(self):

1436

#--------------- Specular Meteor ----------------

1527

# data_pre = self.dataOut.data_pre

1528

# nHeights = self.dataOut.nHeights

1529

# # nProfiles = self.dataOut.data_pre.shape[1]

1530

# # data_param = numpy.zeros((len(pairslist),nProfiles,nHeights))

1531

# data_param = numpy.zeros((len(pairslist),nHeights))

1532

#

1533

# for i in range(len(pairslist)):

1534

# chan0 = data_pre[pairslist[i][0],:]

1535

# chan1 = data_pre[pairslist[i][1],:]

1536

# #calcular correlacion cruzada

1537

# #magnitud es coherencia

1538

# #fase es dif fase

1539

# correl = chan0*numpy.conj(chan1)

1540

# coherence = numpy.abs(correl)/(numpy.abs(chan0)*numpy.abs(chan1))

1541

# phase = numpy.angle(correl)

1542

# # data_param[2*i,:,:] = coherence

1543

# data_param[i,:] = phase

1544

#

1545

# self.dataOut.data_param = data_param

1546

# self.dataOut.groupList = pairslist

1547

data_cspc = self.dataOut.data_pre[1]

1548

ccf = numpy.average(data_cspc,axis=1)

1549

phases = numpy.angle(ccf).T

1550

1551

meteorOps = MeteorOperations()

1552

pairsList = ((0,1),(2,3))

1553

jph = numpy.array([0,0,0,0])

1554

AOAthresh = numpy.pi/8

1555

azimuth = 45

1556

error = numpy.zeros((phases.shape[0],1))

1557

AOA,error = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)

1558

self.dataOut.data_param = AOA.T

1559

1560

class WindProfiler(Operation):

1561

1562

__isConfig = False

1563

1564

__initime = None

1565

__lastdatatime = None

1566

__integrationtime = None

1567

1568

__buffer = None

1569

1570

__dataReady = False

1571

1572

__firstdata = None

1573

1574

n = None

1575

1576

def __init__(self):

1577

Operation.__init__(self)

1578

1579

def __calculateCosDir(self, elev, azim):

1580

zen = (90 - elev)*numpy.pi/180

1581

azim = azim*numpy.pi/180

1582

cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))

1583

cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)

1584

1585

signX = numpy.sign(numpy.cos(azim))

1586

signY = numpy.sign(numpy.sin(azim))

1587

1588

cosDirX = numpy.copysign(cosDirX, signX)

1589

cosDirY = numpy.copysign(cosDirY, signY)

1590

return cosDirX, cosDirY

1591

1592

def __calculateAngles(self, theta_x, theta_y, azimuth):

1593

1594

dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)

1595

zenith_arr = numpy.arccos(dir_cosw)

1596

azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180

1597

1598

dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)

1599

dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)

1600

1601

return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw

1602

1603

def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):

1604

1605

#

1606

if horOnly:

1607

A = numpy.c_[dir_cosu,dir_cosv]

1608

else:

1609

A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]

1610

A = numpy.asmatrix(A)

1611

A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()

1612

1613

return A1

1614

1615

def __correctValues(self, heiRang, phi, velRadial, SNR):

1616

listPhi = phi.tolist()

1617

maxid = listPhi.index(max(listPhi))

1618

minid = listPhi.index(min(listPhi))

1619

1620

rango = range(len(phi))

1621

# rango = numpy.delete(rango,maxid)

1622

1437

1623

heiRang1 = heiRang*math.cos(phi[maxid])

1438

class MeteorDetection(Operation):

1624

heiRangAux = heiRang*math.cos(phi[minid])

1439

'''

1625

indOut = (heiRang1 < heiRangAux[0]).nonzero()

1440

Function DetectMeteors()

1626

heiRang1 = numpy.delete(heiRang1,indOut)

1441

Project developed with paper:

1442

HOLDSWORTH ET AL. 2004

1627

1443

1628

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

1444

Input:

1629

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

1445

self.dataOut.data_pre

1630

1446

1631

for i in rango:

1447

centerReceiverIndex: From the channels, which is the center receiver

1632

x = heiRang*math.cos(phi[i])

1633

y1 = velRadial[i,:]

1634

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

1635

1448

1636

x1 = heiRang1

1449

hei_ref: Height reference for the Beacon signal extraction

1637

y11 = f1(x1)

1450

tauindex:

1451

predefinedPhaseShifts: Predefined phase offset for the voltge signals

1638

1452

1639

y2 = SNR[i,:]

1453

cohDetection: Whether to user Coherent detection or not

1640

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

1454

cohDet_timeStep: Coherent Detection calculation time step

1641

y21 = f2(x1)

1455

cohDet_thresh: Coherent Detection phase threshold to correct phases

1642

1456

1643

velRadial1[i,:] = y11

1457

noise_timeStep: Noise calculation time step

1644

SNR1[i,:] = y21

1458

noise_multiple: Noise multiple to define signal threshold

1645

1459

1646

return heiRang1, velRadial1, SNR1

1460

multDet_timeLimit: Multiple Detection Removal time limit in seconds

1461

multDet_rangeLimit: Multiple Detection Removal range limit in km

1647

1462

1648

def __calculateVelUVW(self, A, velRadial):

1463

phaseThresh: Maximum phase difference between receiver to be consider a meteor

1464

SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor

1649

1465

1650

#Operacion Matricial

1466

hmin: Minimum Height of the meteor to use it in the further wind estimations

1651

# velUVW = numpy.zeros((velRadial.shape[1],3))

1467

hmax: Maximum Height of the meteor to use it in the further wind estimations

1652

# for ind in range(velRadial.shape[1]):

1468

azimuth: Azimuth angle correction

1653

# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])

1654

# velUVW = velUVW.transpose()

1655

velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))

1656

velUVW[:,:] = numpy.dot(A,velRadial)

1657

1469

1470

Affected:

1471

self.dataOut.data_param

1658

1472

1659

return velUVW

1473

Rejection Criteria (Errors):

1474

0: No error; analysis OK

1475

1: SNR < SNR threshold

1476

2: angle of arrival (AOA) ambiguously determined

1477

3: AOA estimate not feasible

1478

4: Large difference in AOAs obtained from different antenna baselines

1479

5: echo at start or end of time series

1480

6: echo less than 5 examples long; too short for analysis

1481

7: echo rise exceeds 0.3s

1482

8: echo decay time less than twice rise time

1483

9: large power level before echo

1484

10: large power level after echo

1485

11: poor fit to amplitude for estimation of decay time

1486

12: poor fit to CCF phase variation for estimation of radial drift velocity

1487

13: height unresolvable echo: not valid height within 70 to 110 km

1488

14: height ambiguous echo: more then one possible height within 70 to 110 km

1489

15: radial drift velocity or projected horizontal velocity exceeds 200 m/s

1490

16: oscilatory echo, indicating event most likely not an underdense echo

1660

1491

1661

def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):

1492

17: phase difference in meteor Reestimation

1662

"""

1663

Function that implements Doppler Beam Swinging (DBS) technique.

1664

1493

1665

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

1494

Data Storage:

1666

Direction correction (if necessary), Ranges and SNR

1495

Meteors for Wind Estimation (8):

1496

Utc Time | Range Height

1497

Azimuth Zenith errorCosDir

1498

VelRad errorVelRad

1499

Phase0 Phase1 Phase2 Phase3

1500

TypeError

1667

1501

1668

Output: Winds estimation (Zonal, Meridional and Vertical)

1502

'''

1669

1503

1670

Parameters affected: Winds, height range, SNR

1504

def run(self, dataOut, hei_ref = None, tauindex = 0,

1671

"""

1505

phaseOffsets = None,

1672

azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)

1506

cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,

1673

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)

1507

noise_timeStep = 4, noise_multiple = 4,

1674

A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)

1508

multDet_timeLimit = 1, multDet_rangeLimit = 3,

1509

phaseThresh = 20, SNRThresh = 5,

1510

hmin = 50, hmax=150, azimuth = 0,

1511

channelPositions = None) :

1675

1512

1676

#Calculo de Componentes de la velocidad con DBS

1677

winds = self.__calculateVelUVW(A,velRadial1)

1678

1513

1679

return winds, heiRang1, SNR1

1514

#Getting Pairslist

1515

if channelPositions == None:

1516

# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T

1517

channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella

1518

meteorOps = MeteorOperations()

1519

pairslist0, distances = meteorOps.getPhasePairs(channelPositions)

1680

1520

1681

def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):

1521

#Get Beacon signal

1522

newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])

1682

1523

1683

posx = numpy.asarray(posx)

1524

if hei_ref != None:

1684

posy = numpy.asarray(posy)

1525

newheis = numpy.where(self.dataOut.heightList>hei_ref)

1685

1526

1686

#Rotacion Inversa para alinear con el azimuth

1527

heiRang = dataOut.getHeiRange()

1687

if azimuth!= None:

1688

azimuth = azimuth*math.pi/180

1689

posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)

1690

posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)

1691

else:

1692

posx1 = posx

1693

posy1 = posy

1694

1528

1695

#Calculo de Distancias

1529

#****************REMOVING HARDWARE PHASE DIFFERENCES***************

1696

distx = numpy.zeros(pairsCrossCorr.size)

1530

# see if the user put in pre defined phase shifts

1697

disty = numpy.zeros(pairsCrossCorr.size)

1531

voltsPShift = self.dataOut.data_pre.copy()

1698

dist = numpy.zeros(pairsCrossCorr.size)

1699

ang = numpy.zeros(pairsCrossCorr.size)

1700

1532

1701

for i in range(pairsCrossCorr.size):

1533

# if predefinedPhaseShifts != None:

1702

distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]

1534

# hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180

1703

disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]

1535

#

1704

dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)

1536

# # elif beaconPhaseShifts:

1705

ang[i] = numpy.arctan2(disty[i],distx[i])

1537

# # #get hardware phase shifts using beacon signal

1706

#Calculo de Matrices

1538

# # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)

1707

nPairs = len(pairs)

1539

# # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)

1708

ang1 = numpy.zeros((nPairs, 2, 1))

1540

#

1709

dist1 = numpy.zeros((nPairs, 2, 1))

1541

# else:

1542

# hardwarePhaseShifts = numpy.zeros(5)

1543

#

1544

# voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')

1545

# for i in range(self.dataOut.data_pre.shape[0]):

1546

# voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])

1710

1547

1711

for j in range(nPairs):

1548

#******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********

1712

dist1[j,0,0] = dist[pairs[j][0]]

1713

dist1[j,1,0] = dist[pairs[j][1]]

1714

ang1[j,0,0] = ang[pairs[j][0]]

1715

ang1[j,1,0] = ang[pairs[j][1]]

1716

1549

1717

return distx,disty, dist1,ang1

1550

#Remove DC

1551

voltsDC = numpy.mean(voltsPShift,1)

1552

voltsDC = numpy.mean(voltsDC,1)

1553

for i in range(voltsDC.shape[0]):

1554

voltsPShift[i] = voltsPShift[i] - voltsDC[i]

1718

1555

1719

def __calculateVelVer(self, phase, lagTRange, _lambda):

1556

#Don't considerate last heights, theyre used to calculate Hardware Phase Shift

1557

voltsPShift = voltsPShift[:,:,:newheis[0][0]]

1720

1558

1721

Ts = lagTRange[1] - lagTRange[0]

1559

#************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********

1722

velW = -_lambda*phase/(4*math.pi*Ts)

1560

#Coherent Detection

1561

if cohDetection:

1562

#use coherent detection to get the net power

1563

cohDet_thresh = cohDet_thresh*numpy.pi/180

1564

voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh)

1723

1565

1724

return velW

1566

#Non-coherent detection!

1567

powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)

1568

#********** END OF COH/NON-COH POWER CALCULATION**********************

1725

1569

1726

def __calculateVelHorDir(self, dist, tau1, tau2, ang):

1570

#********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************

1727

nPairs = tau1.shape[0]

1571

#Get noise

1728

vel = numpy.zeros((nPairs,3,tau1.shape[2]))

1572

noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)

1573

# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)

1574

#Get signal threshold

1575

signalThresh = noise_multiple*noise

1576

#Meteor echoes detection

1577

listMeteors = self.__findMeteors(powerNet, signalThresh)

1578

#******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********

1729

1579

1730

angCos = numpy.cos(ang)

1580

#************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************

1731

angSin = numpy.sin(ang)

1581

#Parameters

1582

heiRange = self.dataOut.getHeiRange()

1583

rangeInterval = heiRange[1] - heiRange[0]

1584

rangeLimit = multDet_rangeLimit/rangeInterval

1585

timeLimit = multDet_timeLimit/self.dataOut.timeInterval

1586

#Multiple detection removals

1587

listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)

1588

#************ END OF REMOVE MULTIPLE DETECTIONS **********************

1732

1589

1733

vel0 = dist*tau1/(2*tau2**2)

1590

#********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************

1734

vel[:,0,:] = (vel0*angCos).sum(axis = 1)

1591

#Parameters

1735

vel[:,1,:] = (vel0*angSin).sum(axis = 1)

1592

phaseThresh = phaseThresh*numpy.pi/180

1593

thresh = [phaseThresh, noise_multiple, SNRThresh]

1594

#Meteor reestimation (Errors N 1, 6, 12, 17)

1595

listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)

1596

# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)

1597

#Estimation of decay times (Errors N 7, 8, 11)

1598

listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)

1599

#******************* END OF METEOR REESTIMATION *******************

1736

1600

1737

ind = numpy.where(numpy.isinf(vel))

1601

#********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************

1738

vel[ind] = numpy.nan

1602

#Calculating Radial Velocity (Error N 15)

1603

radialStdThresh = 10

1604

listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval)

1739

1605

1740

return vel

1606

if len(listMeteors4) > 0:

1607

#Setting New Array

1608

date = self.dataOut.utctime

1609

arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)

1741

1610

1742

def __getPairsAutoCorr(self, pairsList, nChannels):

1611

#Correcting phase offset

1612

if phaseOffsets != None:

1613

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

1614

arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

1743

1615

1744

pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan

1616

#Second Pairslist

1617

pairsList = []

1618

pairx = (0,1)

1619

pairy = (2,3)

1620

pairsList.append(pairx)

1621

pairsList.append(pairy)

1745

1622

1746

for l in range(len(pairsList)):

1623

jph = numpy.array([0,0,0,0])

1747

firstChannel = pairsList[l][0]

1624

h = (hmin,hmax)

1748

secondChannel = pairsList[l][1]

1625

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

1749

1626

1750

#Obteniendo pares de Autocorrelacion

1627

# #Calculate AOA (Error N 3, 4)

1751

if firstChannel == secondChannel:

1628

# #JONES ET AL. 1998

1752

pairsAutoCorr[firstChannel] = int(l)

1629

# error = arrayParameters[:,-1]

1630

# AOAthresh = numpy.pi/8

1631

# phases = -arrayParameters[:,9:13]

1632

# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)

1633

#

1634

# #Calculate Heights (Error N 13 and 14)

1635

# error = arrayParameters[:,-1]

1636

# Ranges = arrayParameters[:,2]

1637

# zenith = arrayParameters[:,5]

1638

# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)

1639

# error = arrayParameters[:,-1]

1640

#********************* END OF PARAMETERS CALCULATION **************************

1753

1641

1754

pairsAutoCorr = pairsAutoCorr.astype(int)

1642

#***************************+ PASS DATA TO NEXT STEP **********************

1643

# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))

1644

self.dataOut.data_param = arrayParameters

1755

1645

1756

pairsCrossCorr = range(len(pairsList))

1646

if arrayParameters == None:

1757

pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)

1647

self.dataOut.flagNoData = True

1648

else:

1649

self.dataOut.flagNoData = True

1758

1650

1759

return pairsAutoCorr, pairsCrossCorr

1651

return

1760

1652

1761

def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):

1653

def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):

1762

"""

1763

Function that implements Spaced Antenna (SA) technique.

1764

1654

1765

Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,

1655

minIndex = min(newheis[0])

1766

Direction correction (if necessary), Ranges and SNR

1656

maxIndex = max(newheis[0])

1767

1657

1768

Output: Winds estimation (Zonal, Meridional and Vertical)

1658

voltage = voltage0[:,:,minIndex:maxIndex+1]

1659

nLength = voltage.shape[1]/n

1660

nMin = 0

1661

nMax = 0

1662

phaseOffset = numpy.zeros((len(pairslist),n))

1769

1663

1770

Parameters affected: Winds

1664

for i in range(n):

1771

"""

1665

nMax += nLength

1772

#Cross Correlation pairs obtained

1666

phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))

1773

pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)

1667

phaseCCF = numpy.mean(phaseCCF, axis = 2)

1774

pairsArray = numpy.array(pairsList)[pairsCrossCorr]

1668

phaseOffset[:,i] = phaseCCF.transpose()

1775

pairsSelArray = numpy.array(pairsSelected)

1669

nMin = nMax

1776

pairs = []

1670

# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)

1777

1671

1778

#Wind estimation pairs obtained

1672

#Remove Outliers

1779

for i in range(pairsSelArray.shape[0]/2):

1673

factor = 2

1780

ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]

1674

wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))

1781

ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]

1675

dw = numpy.std(wt,axis = 1)

1782

pairs.append((ind1,ind2))

1676

dw = dw.reshape((dw.size,1))

1677

ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))

1678

phaseOffset[ind] = numpy.nan

1679

phaseOffset = stats.nanmean(phaseOffset, axis=1)

1783

1680

1784

indtau = tau.shape[0]/2

1681

return phaseOffset

1785

tau1 = tau[:indtau,:]

1786

tau2 = tau[indtau:-1,:]

1787

tau1 = tau1[pairs,:]

1788

tau2 = tau2[pairs,:]

1789

phase1 = tau[-1,:]

1790

1682

1791

#---------------------------------------------------------------------

1683

def __shiftPhase(self, data, phaseShift):

1792

#Metodo Directo

1684

#this will shift the phase of a complex number

1793

distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)

1685

dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)

1794

winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)

1686

return dataShifted

1795

winds = stats.nanmean(winds, axis=0)

1796

#---------------------------------------------------------------------

1797

#Metodo General

1798

# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)

1799

# #Calculo Coeficientes de Funcion de Correlacion

1800

# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)

1801

# #Calculo de Velocidades

1802

# winds = self.calculateVelUV(F,G,A,B,H)

1803

1687

1804

#---------------------------------------------------------------------

1688

def __estimatePhaseDifference(self, array, pairslist):

1805

winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)

1689

nChannel = array.shape[0]

1806

winds = correctFactor*winds

1690

nHeights = array.shape[2]

1807

return winds

1691

numPairs = len(pairslist)

1692

# phaseCCF = numpy.zeros((nChannel, 5, nHeights))

1693

phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))

1808

1694

1809

def __checkTime(self, currentTime, paramInterval, outputInterval):

1695

#Correct phases

1696

derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]

1697

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

1810

1698

1811

dataTime = currentTime + paramInterval

1699

if indDer[0].shape[0] > 0:

1812

deltaTime = dataTime - self.__initime

1700

for i in range(indDer[0].shape[0]):

1701

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])

1702

phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi

1813

1703

1814

if deltaTime >= outputInterval or deltaTime < 0:

1704

# for j in range(numSides):

1815

self.__dataReady = True

1705

# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])

1816

return

1706

# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)

1707

#

1708

#Linear

1709

phaseInt = numpy.zeros((numPairs,1))

1710

angAllCCF = phaseCCF[:,[0,1,3,4],0]

1711

for j in range(numPairs):

1712

fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])

1713

phaseInt[j] = fit[1]

1714

#Phase Differences

1715

phaseDiff = phaseInt - phaseCCF[:,2,:]

1716

phaseArrival = phaseInt.reshape(phaseInt.size)

1817

1717

1818

def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):

1718

#Dealias

1819

'''

1719

phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))

1820

Function that implements winds estimation technique with detected meteors.

1720

# indAlias = numpy.where(phaseArrival > numpy.pi)

1721

# phaseArrival[indAlias] -= 2*numpy.pi

1722

# indAlias = numpy.where(phaseArrival < -numpy.pi)

1723

# phaseArrival[indAlias] += 2*numpy.pi

1821

1724

1822

Input: Detected meteors, Minimum meteor quantity to wind estimation

1725

return phaseDiff, phaseArrival

1823

1726

1824

Output: Winds estimation (Zonal and Meridional)

1727

def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):

1728

#this function will run the coherent detection used in Holdworth et al. 2004 and return the net power

1729

#find the phase shifts of each channel over 1 second intervals

1730

#only look at ranges below the beacon signal

1731

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

1732

numBlocks = int(volts.shape[1]/numProfPerBlock)

1733

numHeights = volts.shape[2]

1734

nChannel = volts.shape[0]

1735

voltsCohDet = volts.copy()

1825

1736

1826

Parameters affected: Winds

1737

pairsarray = numpy.array(pairslist)

1827

'''

1738

indSides = pairsarray[:,1]

1828

# print arrayMeteor.shape

1739

# indSides = numpy.array(range(nChannel))

1829

#Settings

1740

# indSides = numpy.delete(indSides, indCenter)

1830

nInt = (heightMax - heightMin)/2

1741

#

1831

# print nInt

1742

# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)

1832

nInt = int(nInt)

1743

listBlocks = numpy.array_split(volts, numBlocks, 1)

1833

# print nInt

1834

winds = numpy.zeros((2,nInt))*numpy.nan

1835

1744

1836

#Filter errors

1745

startInd = 0

1837

error = numpy.where(arrayMeteor[:,-1] == 0)[0]

1746

endInd = 0

1838

finalMeteor = arrayMeteor[error,:]

1839

1747

1840

#Meteor Histogram

1748

for i in range(numBlocks):

1841

finalHeights = finalMeteor[:,2]

1749

startInd = endInd

1842

hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))

1750

endInd = endInd + listBlocks[i].shape[1]

1843

nMeteorsPerI = hist[0]

1844

heightPerI = hist[1]

1845

1751

1846

#Sort of meteors

1752

arrayBlock = listBlocks[i]

1847

indSort = finalHeights.argsort()

1753

# arrayBlockCenter = listCenter[i]

1848

finalMeteor2 = finalMeteor[indSort,:]

1849

1754

1850

# Calculating winds

1755

#Estimate the Phase Difference

1851

ind1 = 0

1756

phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)

1852

ind2 = 0

1757

#Phase Difference RMS

1758

arrayPhaseRMS = numpy.abs(phaseDiff)

1759

phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)

1760

indPhase = numpy.where(phaseRMSaux==4)

1761

#Shifting

1762

if indPhase[0].shape[0] > 0:

1763

for j in range(indSides.size):

1764

arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())

1765

voltsCohDet[:,startInd:endInd,:] = arrayBlock

1853

1766

1854

for i in range(nInt):

1767

return voltsCohDet

1855

nMet = nMeteorsPerI[i]

1856

ind1 = ind2

1857

ind2 = ind1 + nMet

1858

1768

1859

meteorAux = finalMeteor2[ind1:ind2,:]

1769

def __calculateCCF(self, volts, pairslist ,laglist):

1860

1770

1861

if meteorAux.shape[0] >= meteorThresh:

1771

nHeights = volts.shape[2]

1862

vel = meteorAux[:, 6]

1772

nPoints = volts.shape[1]

1863

zen = meteorAux[:, 4]*numpy.pi/180

1773

voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')

1864

azim = meteorAux[:, 3]*numpy.pi/180

1865

1774

1866

n = numpy.cos(zen)

1775

for i in range(len(pairslist)):

1867

# m = (1 - n**2)/(1 - numpy.tan(azim)**2)

1776

volts1 = volts[pairslist[i][0]]

1868

# l = m*numpy.tan(azim)

1777

volts2 = volts[pairslist[i][1]]

1869

l = numpy.sin(zen)*numpy.sin(azim)

1870

m = numpy.sin(zen)*numpy.cos(azim)

1871

1778

1872

A = numpy.vstack((l, m)).transpose()

1779

for t in range(len(laglist)):

1873

A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())

1780

idxT = laglist[t]

1874

windsAux = numpy.dot(A1, vel)

1781

if idxT >= 0:

1782

vStacked = numpy.vstack((volts2[idxT:,:],

1783

numpy.zeros((idxT, nHeights),dtype='complex')))

1784

else:

1785

vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),

1786

volts2[:(nPoints + idxT),:]))

1787

voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)

1875

1788

1876

winds[0,i] = windsAux[0]

1789

vStacked = None

1877

winds[1,i] = windsAux[1]

1790

return voltsCCF

1878

1791

1879

return winds, heightPerI[:-1]

1792

def __getNoise(self, power, timeSegment, timeInterval):

1793

numProfPerBlock = numpy.ceil(timeSegment/timeInterval)

1794

numBlocks = int(power.shape[0]/numProfPerBlock)

1795

numHeights = power.shape[1]

1880

1796

1881

def techniqueNSM_SA(self, **kwargs):

1797

listPower = numpy.array_split(power, numBlocks, 0)

1882

metArray = kwargs['metArray']

1798

noise = numpy.zeros((power.shape[0], power.shape[1]))

1883

heightList = kwargs['heightList']

1799

noise1 = numpy.zeros((power.shape[0], power.shape[1]))

1884

timeList = kwargs['timeList']

1885

1800

1886

rx_location = kwargs['rx_location']

1801

startInd = 0

1887

groupList = kwargs['groupList']

1802

endInd = 0

1888

azimuth = kwargs['azimuth']

1889

dfactor = kwargs['dfactor']

1890

k = kwargs['k']

1891

1803

1892

azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)

1804

for i in range(numBlocks): #split por canal

1893

d = dist*dfactor

1805

startInd = endInd

1894

#Phase calculation

1806

endInd = endInd + listPower[i].shape[0]

1895

metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)

1896

1807

1897

metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities

1808

arrayBlock = listPower[i]

1809

noiseAux = numpy.mean(arrayBlock, 0)

1810

# noiseAux = numpy.median(noiseAux)

1811

# noiseAux = numpy.mean(arrayBlock)

1812

noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux

1898

1813

1899

velEst = numpy.zeros((heightList.size,2))*numpy.nan

1814

noiseAux1 = numpy.mean(arrayBlock)

1900

azimuth1 = azimuth1*numpy.pi/180

1815

noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1

1901

1816

1902

for i in range(heightList.size):

1817

return noise, noise1

1903

h = heightList[i]

1904

indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]

1905

metHeight = metArray1[indH,:]

1906

if metHeight.shape[0] >= 2:

1907

velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities

1908

iazim = metHeight[:,1].astype(int)

1909

azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths

1910

A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))

1911

A = numpy.asmatrix(A)

1912

A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()

1913

velHor = numpy.dot(A1,velAux)

1914

1818

1915

velEst[i,:] = numpy.squeeze(velHor)

1819

def __findMeteors(self, power, thresh):

1916

return velEst

1820

nProf = power.shape[0]

1821

nHeights = power.shape[1]

1822

listMeteors = []

1917

1823

1918

def __getPhaseSlope(self, metArray, heightList, timeList):

1824

for i in range(nHeights):

1919

meteorList = []

1825

powerAux = power[:,i]

1920

#utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2

1826

threshAux = thresh[:,i]

1921

#Putting back together the meteor matrix

1922

utctime = metArray[:,0]

1923

uniqueTime = numpy.unique(utctime)

1924

1827

1925

phaseDerThresh = 0.5

1828

indUPthresh = numpy.where(powerAux > threshAux)[0]

1926

ippSeconds = timeList[1] - timeList[0]

1829

indDNthresh = numpy.where(powerAux <= threshAux)[0]

1927

sec = numpy.where(timeList>1)[0][0]

1928

nPairs = metArray.shape[1] - 6

1929

nHeights = len(heightList)

1930

1830

1931

for t in uniqueTime:

1831

j = 0

1932

metArray1 = metArray[utctime==t,:]

1933

# phaseDerThresh = numpy.pi/4 #reducir Phase thresh

1934

tmet = metArray1[:,1].astype(int)

1935

hmet = metArray1[:,2].astype(int)

1936

1832

1937

metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))

1833

while (j < indUPthresh.size - 2):

1938

metPhase[:,:] = numpy.nan

1834

if (indUPthresh[j + 2] == indUPthresh[j] + 2):

1939

metPhase[:,hmet,tmet] = metArray1[:,6:].T

1835

indDNAux = numpy.where(indDNthresh > indUPthresh[j])

1836

indDNthresh = indDNthresh[indDNAux]

1940

1837

1941

#Delete short trails

1838

if (indDNthresh.size > 0):

1942

metBool = ~numpy.isnan(metPhase[0,:,:])

1839

indEnd = indDNthresh[0] - 1

1943

heightVect = numpy.sum(metBool, axis = 1)

1840

indInit = indUPthresh[j]

1944

metBool[heightVect<sec,:] = False

1945

metPhase[:,heightVect<sec,:] = numpy.nan

1946

1841

1947

#Derivative

1842

meteor = powerAux[indInit:indEnd + 1]

1948

metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])

1843

indPeak = meteor.argmax() + indInit

1949

phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))

1844

FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))

1950

metPhase[phDerAux] = numpy.nan

1951

1845

1952

#--------------------------METEOR DETECTION -----------------------------------------

1846

listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!

1953

indMet = numpy.where(numpy.any(metBool,axis=1))[0]

1847

j = numpy.where(indUPthresh == indEnd)[0] + 1

1848

else: j+=1

1849

else: j+=1

1954

1850

1955

for p in numpy.arange(nPairs):

1851

return listMeteors

1956

phase = metPhase[p,:,:]

1957

phDer = metDer[p,:,:]

1958

1852

1959

for h in indMet:

1853

def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):

1960

height = heightList[h]

1961

phase1 = phase[h,:] #82

1962

phDer1 = phDer[h,:]

1963

1854

1964

phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap

1855

arrayMeteors = numpy.asarray(listMeteors)

1856

listMeteors1 = []

1965

1857

1966

indValid = numpy.where(~numpy.isnan(phase1))[0]

1858

while arrayMeteors.shape[0] > 0:

1967

initMet = indValid[0]

1859

FLAs = arrayMeteors[:,4]

1968

endMet = 0

1860

maxFLA = FLAs.argmax()

1861

listMeteors1.append(arrayMeteors[maxFLA,:])

1969

1862

1970

for i in range(len(indValid)-1):

1863

MeteorInitTime = arrayMeteors[maxFLA,1]

1864

MeteorEndTime = arrayMeteors[maxFLA,3]

1865

MeteorHeight = arrayMeteors[maxFLA,0]

1971

1866

1972

#Time difference

1867

#Check neighborhood

1973

inow = indValid[i]

1868

maxHeightIndex = MeteorHeight + rangeLimit

1974

inext = indValid[i+1]

1869

minHeightIndex = MeteorHeight - rangeLimit

1975

idiff = inext - inow

1870

minTimeIndex = MeteorInitTime - timeLimit

1976

#Phase difference

1871

maxTimeIndex = MeteorEndTime + timeLimit

1977

phDiff = numpy.abs(phase1[inext] - phase1[inow])

1978

1872

1979

if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor

1873

#Check Heights

1980

sizeTrail = inow - initMet + 1

1874

indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)

1981

if sizeTrail>3*sec: #Too short meteors

1875

indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)

1982

x = numpy.arange(initMet,inow+1)*ippSeconds

1876

indBoth = numpy.where(numpy.logical_and(indTime,indHeight))

1983

y = phase1[initMet:inow+1]

1984

ynnan = ~numpy.isnan(y)

1985

x = x[ynnan]

1986

y = y[ynnan]

1987

slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)

1988

ylin = x*slope + intercept

1989

rsq = r_value**2

1990

if rsq > 0.5:

1991

vel = slope#*height*1000/(k*d)

1992

estAux = numpy.array([utctime,p,height, vel, rsq])

1993

meteorList.append(estAux)

1994

initMet = inext

1995

metArray2 = numpy.array(meteorList)

1996

1877

1997

return metArray2

1878

arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)

1998

1879

1999

def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):

1880

return listMeteors1

2000

1881

2001

azimuth1 = numpy.zeros(len(pairslist))

1882

def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):

2002

dist = numpy.zeros(len(pairslist))

1883

numHeights = volts.shape[2]

1884

nChannel = volts.shape[0]

2003

1885

2004

for i in range(len(rx_location)):

1886

thresholdPhase = thresh[0]

2005

ch0 = pairslist[i][0]

1887

thresholdNoise = thresh[1]

2006

ch1 = pairslist[i][1]

1888

thresholdDB = float(thresh[2])

2007

1889

2008

diffX = rx_location[ch0][0] - rx_location[ch1][0]

1890

thresholdDB1 = 10**(thresholdDB/10)

2009

diffY = rx_location[ch0][1] - rx_location[ch1][1]

1891

pairsarray = numpy.array(pairslist)

2010

azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi

1892

indSides = pairsarray[:,1]

2011

dist[i] = numpy.sqrt(diffX**2 + diffY**2)

2012

1893

2013

azimuth1 -= azimuth0

1894

pairslist1 = list(pairslist)

2014

return azimuth1, dist

1895

pairslist1.append((0,1))

1896

pairslist1.append((3,4))

2015

1897

2016

def techniqueNSM_DBS(self, **kwargs):

1898

listMeteors1 = []

2017

metArray = kwargs['metArray']

1899

listPowerSeries = []

2018

heightList = kwargs['heightList']

1900

listVoltageSeries = []

2019

timeList = kwargs['timeList']

1901

#volts has the war data

2020

zenithList = kwargs['zenithList']

2021

nChan = numpy.max(cmet) + 1

2022

nHeights = len(heightList)

2023

1902

2024

utctime = metArray[:,0]

1903

if frequency == 30e6:

2025

cmet = metArray[:,1]

1904

timeLag = 45*10**-3

2026

hmet = metArray1[:,3].astype(int)

1905

else:

2027

h1met = heightList[hmet]*zenithList[cmet]

1906

timeLag = 15*10**-3

2028

vmet = metArray1[:,5]

1907

lag = numpy.ceil(timeLag/timeInterval)

2029

1908

2030

for i in range(~~nHeights~~ - 1):

1909

for i in range(len(listMeteors)):

2031

hmin = heightList[i]

2032

hmax = heightList[i + 1]

2033

1910

2034

vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]

1911

###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################

1912

meteorAux = numpy.zeros(16)

2035

1913

1914

#Loading meteor Data (mHeight, mStart, mPeak, mEnd)

1915

mHeight = listMeteors[i][0]

1916

mStart = listMeteors[i][1]

1917

mPeak = listMeteors[i][2]

1918

mEnd = listMeteors[i][3]

2036

1919

1920

#get the volt data between the start and end times of the meteor

1921

meteorVolts = volts[:,mStart:mEnd+1,mHeight]

1922

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

2037

1923

2038

return data_output

1924

#3.6. Phase Difference estimation

1925

phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)

2039

1926

2040

def run(self, dataOut, technique, **kwargs):

1927

#3.7. Phase difference removal & meteor start, peak and end times reestimated

1928

#meteorVolts0.- all Channels, all Profiles

1929

meteorVolts0 = volts[:,:,mHeight]

1930

meteorThresh = noise[:,mHeight]*thresholdNoise

1931

meteorNoise = noise[:,mHeight]

1932

meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting

1933

powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power

2041

1934

2042

param = dataOut.data_param

1935

#Times reestimation

2043

if dataOut.abscissaList != None:

1936

mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]

2044

absc = dataOut.abscissaList[:-1]

1937

if mStart1.size > 0:

2045

noise = dataOut.noise

1938

mStart1 = mStart1[-1] + 1

2046

heightList = dataOut.heightList

2047

SNR = dataOut.data_SNR

2048

1939

2049

if technique == 'DBS':

1940

else:

1941

mStart1 = mPeak

2050

1942

2051

if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):

1943

mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1

2052

theta_x = numpy.array(kwargs['dirCosx'])

1944

mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]

2053

theta_y = numpy.array(kwargs['dirCosy'])

1945

if mEndDecayTime1.size == 0:

1946

mEndDecayTime1 = powerNet0.size

2054

else:

1947

else:

2055

elev = numpy.array(kwargs['elevation'])

1948

mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1

2056

azim = numpy.array(kwargs['azimuth'])

1949

# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()

2057

theta_x, theta_y = self.__calculateCosDir(elev, azim)

2058

azimuth = kwargs['correctAzimuth']

2059

if kwargs.has_key('horizontalOnly'):

2060

horizontalOnly = kwargs['horizontalOnly']

2061

else: horizontalOnly = False

2062

if kwargs.has_key('correctFactor'):

2063

correctFactor = kwargs['correctFactor']

2064

else: correctFactor = 1

2065

if kwargs.has_key('channelList'):

2066

channelList = kwargs['channelList']

2067

if len(channelList) == 2:

2068

horizontalOnly = True

2069

arrayChannel = numpy.array(channelList)

2070

param = param[arrayChannel,:,:]

2071

theta_x = theta_x[arrayChannel]

2072

theta_y = theta_y[arrayChannel]

2073

1950

2074

velRadial0 = param[:,1,:] #Radial velocity

1951

#meteorVolts1.- all Channels, from start to end

2075

dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function

1952

meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]

2076

dataOut.utctimeInit = dataOut.utctime

1953

meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]

2077

dataOut.outputInterval = dataOut.paramInterval

1954

if meteorVolts2.shape[1] == 0:

1955

meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]

1956

meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)

1957

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)

1958

##################### END PARAMETERS REESTIMATION #########################

2078

1959

2079

elif technique == 'SA':

1960

##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################

1961

# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis

1962

if meteorVolts2.shape[1] > 0:

1963

#Phase Difference re-estimation

1964

phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation

1965

# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)

1966

meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])

1967

phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))

1968

meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting

2080

1969

2081

~~#Parameters~~

1970

#Phase Difference RMS

2082

position_x = kwargs['positionX']

1971

phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))

2083

position_y = kwargs['positionY']

1972

powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)

2084

azimuth = kwargs['azimuth']

1973

#Data from Meteor

1974

mPeak1 = powerNet1.argmax() + mStart1

1975

mPeakPower1 = powerNet1.max()

1976

noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])

1977

mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux

1978

Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])

1979

Meteor1 = numpy.hstack((Meteor1,phaseDiffint))

1980

PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]

1981

#Vectorize

1982

meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]

1983

meteorAux[7:11] = phaseDiffint[0:4]

1984

1985

#Rejection Criterions

1986

if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation

1987

meteorAux[-1] = 17

1988

elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB

1989

meteorAux[-1] = 1

2085

1990

2086

if kwargs.has_key('crosspairsList'):

2087

pairs = kwargs['crosspairsList']

2088

else:

2089

pairs = None

2090

1991

2091

if kwargs.has_key('correctFactor'):

2092

correctFactor = kwargs['correctFactor']

2093

else:

1992

else:

2094

correctFactor = 1

1993

meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]

1994

meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis

1995

PowerSeries = 0

2095

1996

2096

tau = dataOut.data_param

1997

listMeteors1.append(meteorAux)

2097

_lambda = dataOut.C/dataOut.frequency

1998

listPowerSeries.append(PowerSeries)

2098

pairsList = dataOut.groupList

1999

listVoltageSeries.append(meteorVolts1)

2099

nChannels = dataOut.nChannels

2100

2000

2101

dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)

2001

return listMeteors1, listPowerSeries, listVoltageSeries

2102

dataOut.utctimeInit = dataOut.utctime

2103

dataOut.outputInterval = dataOut.timeInterval

2104

2002

2105

elif technique == 'Meteors':

2003

def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):

2106

dataOut.flagNoData = True

2107

self.__dataReady = False

2108

2004

2109

if kwargs.has_key('nHours'):

2005

threshError = 10

2110

nHours = kwargs['nHours']

2006

#Depending if it is 30 or 50 MHz

2007

if frequency == 30e6:

2008

timeLag = 45*10**-3

2111

else:

2009

else:

2112

nHours = 1

2010

timeLag = 15*10**-3

2011

lag = numpy.ceil(timeLag/timeInterval)

2113

2012

2114

if kwargs.has_key('meteorsPerBin'):

2013

listMeteors1 = []

2115

meteorThresh = kwargs['meteorsPerBin']

2116

else:

2117

meteorThresh = 6

2118

2014

2119

if kwargs.has_key('hmin'):

2015

for i in range(len(listMeteors)):

2120

hmin = kwargs['hmin']

2016

meteorPower = listPower[i]

2121

else: hmin = 70

2017

meteorAux = listMeteors[i]

2122

if kwargs.has_key('hmax'):

2123

hmax = kwargs['hmax']

2124

else: hmax = 110

2125

2018

2126

dataOut.outputInterval = nHours*3600

2019

if meteorAux[-1] == 0:

2127

2020

2128

if self.__isConfig == False:

2021

try:

2129

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

2022

indmax = meteorPower.argmax()

2130

#Get Initial LTC time

2023

indlag = indmax + lag

2131

self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

2132

self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

2133

2024

2134

self.__isConfig = True

2025

y = meteorPower[indlag:]

2026

x = numpy.arange(0, y.size)*timeLag

2135

2027

2136

if self.__buffer == None:

2028

#first guess

2137

self.__buffer = dataOut.data_param

2029

a = y[0]

2138

self.__firstdata = copy.copy(dataOut)

2030

tau = timeLag

2031

#exponential fit

2032

popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])

2033

y1 = self.__exponential_function(x, *popt)

2034

#error estimation

2035

error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))

2139

2036

2140

else:

2037

decayTime = popt[1]

2141

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

2038

riseTime = indmax*timeInterval

2039

meteorAux[11:13] = [decayTime, error]

2142

2040

2143

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

2041

#Table items 7, 8 and 11

2042

if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s

2043

meteorAux[-1] = 7

2044

elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time

2045

meteorAux[-1] = 8

2046

if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time

2047

meteorAux[-1] = 11

2144

2048

2145

if self.__dataReady:

2146

dataOut.utctimeInit = self.__initime

2147

2049

2148

self.__initime += dataOut.outputInterval #to erase time offset

2050

except:

2051

meteorAux[-1] = 11

2149

2052

2150

dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)

2151

dataOut.flagNoData = False

2152

self.__buffer = None

2153

2053

2154

elif technique == 'Meteors1':

2054

listMeteors1.append(meteorAux)

2155

dataOut.flagNoData = True

2156

self.__dataReady = False

2157

2055

2158

if kwargs.has_key('nMins'):

2056

return listMeteors1

2159

nMins = kwargs['nMins']

2160

else: nMins = 20

2161

if kwargs.has_key('rx_location'):

2162

rx_location = kwargs['rx_location']

2163

else: rx_location = [(0,1),(1,1),(1,0)]

2164

if kwargs.has_key('azimuth'):

2165

azimuth = kwargs['azimuth']

2166

else: azimuth = 51

2167

if kwargs.has_key('dfactor'):

2168

dfactor = kwargs['dfactor']

2169

if kwargs.has_key('mode'):

2170

mode = kwargs['mode']

2171

else: mode = 'SA'

2172

2057

2173

#Borrar luego esto

2058

#Exponential Function

2174

if dataOut.groupList == None:

2175

dataOut.groupList = [(0,1),(0,2),(1,2)]

2176

groupList = dataOut.groupList

2177

C = 3e8

2178

freq = 50e6

2179

lamb = C/freq

2180

k = 2*numpy.pi/lamb

2181

2059

2182

timeList = dataOut.abscissaList

2060

def __exponential_function(self, x, a, tau):

2183

heightList = dataOut.heightList

2061

y = a*numpy.exp(-x/tau)

2062

return y

2184

2063

2185

if self.__isConfig == False:

2064

def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):

2186

dataOut.outputInterval = nMins*60

2065

2187

# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)

2066

pairslist1 = list(pairslist)

2188

#Get Initial LTC time

2067

pairslist1.append((0,1))

2189

initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)

2068

pairslist1.append((3,4))

2190

minuteAux = initime.minute

2069

numPairs = len(pairslist1)

2191

minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)

2070

#Time Lag

2192

self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()

2071

timeLag = 45*10**-3

2072

c = 3e8

2073

lag = numpy.ceil(timeLag/timeInterval)

2074

freq = 30e6

2193

2075

2194

self.__isConfig = True

2076

listMeteors1 = []

2195

2077

2196

if self.__buffer == None:

2078

for i in range(len(listMeteors)):

2197

self.__buffer = dataOut.data_param

2079

meteorAux = listMeteors[i]

2198

self.__firstdata = copy.copy(dataOut)

2080

if meteorAux[-1] == 0:

2081

mStart = listMeteors[i][1]

2082

mPeak = listMeteors[i][2]

2083

mLag = mPeak - mStart + lag

2199

2084

2200

else:

2085

#get the volt data between the start and end times of the meteor

2201

self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))

2086

meteorVolts = listVolts[i]

2087

meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)

2202

2088

2203

self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready

2089

#Get CCF

2090

allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])

2204

2091

2205

if self.__dataReady:

2092

#Method 2

2206

dataOut.utctimeInit = self.__initime

2093

slopes = numpy.zeros(numPairs)

2207

self.__initime += dataOut.outputInterval #to erase time offset

2094

time = numpy.array([-2,-1,1,2])*timeInterval

2095

angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])

2208

2096

2209

metArray = self.__buffer

2097

#Correct phases

2210

if mode == 'SA':

2098

derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]

2211

dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)

2099

indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)

2212

elif mode == 'DBS':

2213

dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)

2214

dataOut.data_output = dataOut.data_output.T

2215

dataOut.flagNoData = False

2216

self.__buffer = None

2217

2100

2218

return

2101

if indDer[0].shape[0] > 0:

2102

for i in range(indDer[0].shape[0]):

2103

signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])

2104

angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi

2219

2105

2220

class EWDriftsEstimation(Operation):

2106

# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))

2107

for j in range(numPairs):

2108

fit = stats.linregress(time, angAllCCF[j,:])

2109

slopes[j] = fit[0]

2221

2110

2222

def __init__(self):

2111

#Remove Outlier

2223

Operation.__init__(self)

2112

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

2113

# slopes = numpy.delete(slopes,indOut)

2114

# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))

2115

# slopes = numpy.delete(slopes,indOut)

2224

2116

2225

def __correctValues(self, heiRang, phi, velRadial, SNR):

2117

radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)

2226

listPhi = phi.tolist()

2118

radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)

2227

maxid = listPhi.index(max(listPhi))

2119

meteorAux[-2] = radialError

2228

minid = listPhi.index(min(listPhi))

2120

meteorAux[-3] = radialVelocity

2229

2121

2230

rango = range(len(phi))

2122

#Setting Error

2231

# rango = numpy.delete(rango,maxid)

2123

#Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s

2124

if numpy.abs(radialVelocity) > 200:

2125

meteorAux[-1] = 15

2126

#Number 12: Poor fit to CCF variation for estimation of radial drift velocity

2127

elif radialError > radialStdThresh:

2128

meteorAux[-1] = 12

2232

2129

2233

heiRang1 = heiRang*math.cos(phi[maxid])

2130

listMeteors1.append(meteorAux)

2234

heiRangAux = heiRang*math.cos(phi[minid])

2131

return listMeteors1

2235

indOut = (heiRang1 < heiRangAux[0]).nonzero()

2236

heiRang1 = numpy.delete(heiRang1,indOut)

2237

2132

2238

velRadial1 = numpy.zeros([len(phi),len(heiRang1)])

2133

def __setNewArrays(self, listMeteors, date, heiRang):

2239

SNR1 = numpy.zeros([len(phi),len(heiRang1)])

2240

2134

2241

for i in rango:

2135

#New arrays

2242

x = heiRang*math.cos(phi[i])

2136

arrayMeteors = numpy.array(listMeteors)

2243

y1 = velRadial[i,:]

2137

arrayParameters = numpy.zeros((len(listMeteors), 13))

2244

f1 = interpolate.interp1d(x,y1,kind = 'cubic')

2245

2138

2246

x1 = heiRang1

2139

#Date inclusion

2247

y11 = f1(x1)

2140

# date = re.findall(r'\((.*?)\)', date)

2141

# date = date[0].split(',')

2142

# date = map(int, date)

2143

#

2144

# if len(date)<6:

2145

# date.append(0)

2146

#

2147

# date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]

2148

# arrayDate = numpy.tile(date, (len(listMeteors), 1))

2149

arrayDate = numpy.tile(date, (len(listMeteors)))

2248

2150

2249

y2 = SNR[i,:]

2151

#Meteor array

2250

f2 = interpolate.interp1d(x,y2,kind = 'cubic')

2152

# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]

2251

y21 = f2(x1)

2153

# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))

2252

2154

2253

velRadial1[i,:] = y11

2155

#Parameters Array

2254

SNR1[i,:] = y21

2156

arrayParameters[:,0] = arrayDate #Date

2157

arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range

2158

arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error

2159

arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases

2160

arrayParameters[:,-1] = arrayMeteors[:,-1] #Error

2255

2161

2256

return heiRang1, velRadial1, SNR1

2257

2162

2258

def run(self, dataOut, zenith, zenithCorrection):

2163

return arrayParameters

2259

heiRang = dataOut.heightList

2260

velRadial = dataOut.data_param[:,3,:]

2261

SNR = dataOut.data_SNR

2262

2164

2263

zenith = numpy.array(zenith)

2165

class CorrectMeteorPhases(Operation):

2264

zenith -= zenithCorrection

2265

zenith *= numpy.pi/180

2266

2166

2267

heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)

2167

def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):

2268

2168

2269

alp = zenith[0]

2169

arrayParameters = dataOut.data_param

2270

bet = ~~zenith~~[1]

2170

pairsList = []

2171

pairx = (0,1)

2172

pairy = (2,3)

2173

pairsList.append(pairx)

2174

pairsList.append(pairy)

2175

jph = numpy.zeros(4)

2271

2176

2272

w_w = velRadial1[0,:]

2177

phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180

2273

w_e = velRadial1[1,:]

2178

# arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)

2179

arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))

2274

2180

2275

w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))

2181

meteorOps = MeteorOperations()

2276

u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))

2182

if channelPositions == None:

2183

# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T

2184

channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella

2277

2185

2278

winds = numpy.vstack((u,w))

2186

pairslist0, distances = meteorOps.getPhasePairs(channelPositions)

2187

h = (hmin,hmax)

2279

2188

2280

dataOut.heightList = heiRang1

2189

arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)

2281

dataOut.data_output = winds

2282

dataOut.data_SNR = SNR1

2283

2190

2284

dataOut.utctimeInit = dataOut.utctime

2191

dataOut.data_param = arrayParameters

2285

dataOut.outputInterval = dataOut.timeInterval

2286

return

2192

return

2287

2193

2288

class PhaseCalibration(Operation):

2194

class PhaseCalibration(Operation):

2710

pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]

2616

pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]

2711

2617

2712

return pairslist, distances

2618

return pairslist, distances

2619

# def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):

2620

#

2621

# arrayAOA = numpy.zeros((phases.shape[0],3))

2622

# cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)

2623

#

2624

# arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)

2625

# cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)

2626

# arrayAOA[:,2] = cosDirError

2627

#

2628

# azimuthAngle = arrayAOA[:,0]

2629

# zenithAngle = arrayAOA[:,1]

2630

#

2631

# #Setting Error

2632

# #Number 3: AOA not fesible

2633

# indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]

2634

# error[indInvalid] = 3

2635

# #Number 4: Large difference in AOAs obtained from different antenna baselines

2636

# indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]

2637

# error[indInvalid] = 4

2638

# return arrayAOA, error

2639

#

2640

# def __getDirectionCosines(self, arrayPhase, pairsList):

2641

#

2642

# #Initializing some variables

2643

# ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi

2644

# ang_aux = ang_aux.reshape(1,ang_aux.size)

2645

#

2646

# cosdir = numpy.zeros((arrayPhase.shape[0],2))

2647

# cosdir0 = numpy.zeros((arrayPhase.shape[0],2))

2648

#

2649

#

2650

# for i in range(2):

2651

# #First Estimation

2652

# phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]

2653

# #Dealias

2654

# indcsi = numpy.where(phi0_aux > numpy.pi)

2655

# phi0_aux[indcsi] -= 2*numpy.pi

2656

# indcsi = numpy.where(phi0_aux < -numpy.pi)

2657

# phi0_aux[indcsi] += 2*numpy.pi

2658

# #Direction Cosine 0

2659

# cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)

2660

#

2661

# #Most-Accurate Second Estimation

2662

# phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]

2663

# phi1_aux = phi1_aux.reshape(phi1_aux.size,1)

2664

# #Direction Cosine 1

2665

# cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)

2666

#

2667

# #Searching the correct Direction Cosine

2668

# cosdir0_aux = cosdir0[:,i]

2669

# cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)

2670

# #Minimum Distance

2671

# cosDiff = (cosdir1 - cosdir0_aux)**2

2672

# indcos = cosDiff.argmin(axis = 1)

2673

# #Saving Value obtained

2674

# cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]

2675

#

2676

# return cosdir0, cosdir

2677

#

2678

# def __calculateAOA(self, cosdir, azimuth):

2679

# cosdirX = cosdir[:,0]

2680

# cosdirY = cosdir[:,1]

2681

#

2682

# zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi

2683

# azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east

2684

# angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()

2685

#

2686

# return angles

2687

#

2688

# def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):

2689

#

2690

# Ramb = 375 #Ramb = c/(2*PRF)

2691

# Re = 6371 #Earth Radius

2692

# heights = numpy.zeros(Ranges.shape)

2693

#

2694

# R_aux = numpy.array([0,1,2])*Ramb

2695

# R_aux = R_aux.reshape(1,R_aux.size)

2696

#

2697

# Ranges = Ranges.reshape(Ranges.size,1)

2698

#

2699

# Ri = Ranges + R_aux

2700

# hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re

2701

#

2702

# #Check if there is a height between 70 and 110 km

2703

# h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)

2704

# ind_h = numpy.where(h_bool == 1)[0]

2705

#

2706

# hCorr = hi[ind_h, :]

2707

# ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))

2708

#

2709

# hCorr = hi[ind_hCorr]

2710

# heights[ind_h] = hCorr

2711

#

2712

# #Setting Error

2713

# #Number 13: Height unresolvable echo: not valid height within 70 to 110 km

2714

# #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km

2715

#

2716

# indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]

2717

# error[indInvalid2] = 14

2718

# indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]

2719

# error[indInvalid1] = 13

2720

#

2721

# return heights, error

2722

No newline at end of file

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages