##// END OF EJS Templates
schain-jc
RSS
Add changes made by J. Oscanoa (jasmet)
Juan C. Espinoza -
r989:57782215c6e2
parent child
Show More
Show file before | Show file after | Hide comments
@@ -174,5 +174,5 class CorrelationProc(ProcessingUnit):
174 self.dataOut.lagRange = numpy.array(lags)*delta
174 self.dataOut.lagRange = numpy.array(lags)*delta
175 # self.dataOut.nCohInt = self.dataIn.nCohInt*nAvg
175 # self.dataOut.nCohInt = self.dataIn.nCohInt*nAvg
176 self.dataOut.flagNoData = False
176 self.dataOut.flagNoData = False
177 a = self.dataOut.normFactor
177 # a = self.dataOut.normFactor
178 return
178 return
Show file before | Show file after | Hide comments
@@ -1017,33 +1017,55 class WindProfiler(Operation):
1017 def techniqueNSM_DBS(self, **kwargs):
1017 def techniqueNSM_DBS(self, **kwargs):
1018 metArray = kwargs['metArray']
1018 metArray = kwargs['metArray']
1019 heightList = kwargs['heightList']
1019 heightList = kwargs['heightList']
1020 timeList = kwargs['timeList']
1020 timeList = kwargs['timeList']
1021 zenithList = kwargs['zenithList']
1021 azimuth = kwargs['azimuth']
1022 theta_x = numpy.array(kwargs['theta_x'])
1023 theta_y = numpy.array(kwargs['theta_y'])
1024
1025 utctime = metArray[:,0]
1026 cmet = metArray[:,1].astype(int)
1027 hmet = metArray[:,3].astype(int)
1028 SNRmet = metArray[:,4]
1029 vmet = metArray[:,5]
1030 spcmet = metArray[:,6]
1031
1022 nChan = numpy.max(cmet) + 1
1032 nChan = numpy.max(cmet) + 1
1023 nHeights = len(heightList)
1033 nHeights = len(heightList)
1024
1034
1025 utctime = metArray[:,0]
1035 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
1026 cmet = metArray[:,1]
1036 hmet = heightList[hmet]
1027 hmet = metArray1[:,3].astype(int)
1037 h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
1028 h1met = heightList[hmet]*zenithList[cmet]
1038
1029 vmet = metArray1[:,5]
1039 velEst = numpy.zeros((heightList.size,2))*numpy.nan
1030
1040
1031 for i in range(nHeights - 1):
1041 for i in range(nHeights - 1):
1032 hmin = heightList[i]
1042 hmin = heightList[i]
1033 hmax = heightList[i + 1]
1043 hmax = heightList[i + 1]
1034
1044
1035 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
1045 thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
1036
1046 indthisH = numpy.where(thisH)
1037
1047
1038
1048 if numpy.size(indthisH) > 3:
1039 return data_output
1049
1050 vel_aux = vmet[thisH]
1051 chan_aux = cmet[thisH]
1052 cosu_aux = dir_cosu[chan_aux]
1053 cosv_aux = dir_cosv[chan_aux]
1054 cosw_aux = dir_cosw[chan_aux]
1055
1056 nch = numpy.size(numpy.unique(chan_aux))
1057 if nch > 1:
1058 A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
1059 velEst[i,:] = numpy.dot(A,vel_aux)
1060
1061 return velEst
1040
1062
1041 def run(self, dataOut, technique, hmin=70, hmax=110, nHours=1, **kwargs):
1063 def run(self, dataOut, technique, **kwargs):
1042
1064
1043 param = dataOut.data_param
1065 param = dataOut.data_param
1044 if dataOut.abscissaList != None:
1066 if dataOut.abscissaList != None:
1045 absc = dataOut.abscissaList[:-1]
1067 absc = dataOut.abscissaList[:-1]
1046 #noise = dataOut.noise
1068 noise = dataOut.noise
1047 heightList = dataOut.heightList
1069 heightList = dataOut.heightList
1048 SNR = dataOut.data_SNR
1070 SNR = dataOut.data_SNR
1049
1071
@@ -1153,11 +1175,15 class WindProfiler(Operation):
1153 else: rx_location = [(0,1),(1,1),(1,0)]
1175 else: rx_location = [(0,1),(1,1),(1,0)]
1154 if kwargs.has_key('azimuth'):
1176 if kwargs.has_key('azimuth'):
1155 azimuth = kwargs['azimuth']
1177 azimuth = kwargs['azimuth']
1156 else: azimuth = 51
1178 else: azimuth = 51.06
1157 if kwargs.has_key('dfactor'):
1179 if kwargs.has_key('dfactor'):
1158 dfactor = kwargs['dfactor']
1180 dfactor = kwargs['dfactor']
1159 if kwargs.has_key('mode'):
1181 if kwargs.has_key('mode'):
1160 mode = kwargs['mode']
1182 mode = kwargs['mode']
1183 if kwargs.has_key('theta_x'):
1184 theta_x = kwargs['theta_x']
1185 if kwargs.has_key('theta_y'):
1186 theta_y = kwargs['theta_y']
1161 else: mode = 'SA'
1187 else: mode = 'SA'
1162
1188
1163 #Borrar luego esto
1189 #Borrar luego esto
@@ -1200,7 +1226,7 class WindProfiler(Operation):
1200 if mode == 'SA':
1226 if mode == 'SA':
1201 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1227 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1202 elif mode == 'DBS':
1228 elif mode == 'DBS':
1203 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
1229 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
1204 dataOut.data_output = dataOut.data_output.T
1230 dataOut.data_output = dataOut.data_output.T
1205 dataOut.flagNoData = False
1231 dataOut.flagNoData = False
1206 self.__buffer = None
1232 self.__buffer = None
@@ -1277,25 +1303,26 class EWDriftsEstimation(Operation):
1277
1303
1278 class NonSpecularMeteorDetection(Operation):
1304 class NonSpecularMeteorDetection(Operation):
1279
1305
1280 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1306 def run(self, dataOut, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1281 data_acf = self.dataOut.data_pre[0]
1307 data_acf = dataOut.data_pre[0]
1282 data_ccf = self.dataOut.data_pre[1]
1308 data_ccf = dataOut.data_pre[1]
1283
1309 pairsList = dataOut.groupList[1]
1284 lamb = self.dataOut.C/self.dataOut.frequency
1310
1285 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
1311 lamb = dataOut.C/dataOut.frequency
1286 paramInterval = self.dataOut.paramInterval
1312 tSamp = dataOut.ippSeconds*dataOut.nCohInt
1287
1313 paramInterval = dataOut.paramInterval
1314
1288 nChannels = data_acf.shape[0]
1315 nChannels = data_acf.shape[0]
1289 nLags = data_acf.shape[1]
1316 nLags = data_acf.shape[1]
1290 nProfiles = data_acf.shape[2]
1317 nProfiles = data_acf.shape[2]
1291 nHeights = self.dataOut.nHeights
1318 nHeights = dataOut.nHeights
1292 nCohInt = self.dataOut.nCohInt
1319 nCohInt = dataOut.nCohInt
1293 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
1320 sec = numpy.round(nProfiles/dataOut.paramInterval)
1294 heightList = self.dataOut.heightList
1321 heightList = dataOut.heightList
1295 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
1322 ippSeconds = dataOut.ippSeconds*dataOut.nCohInt*dataOut.nAvg
1296 utctime = self.dataOut.utctime
1323 utctime = dataOut.utctime
1297
1324
1298 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1325 dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1299
1326
1300 #------------------------ SNR --------------------------------------
1327 #------------------------ SNR --------------------------------------
1301 power = data_acf[:,0,:,:].real
1328 power = data_acf[:,0,:,:].real
@@ -1308,6 +1335,7 class NonSpecularMeteorDetection(Operation):
1308 SNRdB = 10*numpy.log10(SNR)
1335 SNRdB = 10*numpy.log10(SNR)
1309
1336
1310 if mode == 'SA':
1337 if mode == 'SA':
1338 dataOut.groupList = dataOut.groupList[1]
1311 nPairs = data_ccf.shape[0]
1339 nPairs = data_ccf.shape[0]
1312 #---------------------- Coherence and Phase --------------------------
1340 #---------------------- Coherence and Phase --------------------------
1313 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
1341 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
@@ -1315,8 +1343,8 class NonSpecularMeteorDetection(Operation):
1315 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1343 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1316
1344
1317 for p in range(nPairs):
1345 for p in range(nPairs):
1318 ch0 = self.dataOut.groupList[p][0]
1346 ch0 = pairsList[p][0]
1319 ch1 = self.dataOut.groupList[p][1]
1347 ch1 = pairsList[p][1]
1320 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1348 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1321 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1349 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1322 # phase1[p,:,:] = numpy.angle(ccf) #median filter
1350 # phase1[p,:,:] = numpy.angle(ccf) #median filter
@@ -1388,16 +1416,18 class NonSpecularMeteorDetection(Operation):
1388 data_param[:,6:] = phase[:,tmet,hmet].T
1416 data_param[:,6:] = phase[:,tmet,hmet].T
1389
1417
1390 elif mode == 'DBS':
1418 elif mode == 'DBS':
1391 self.dataOut.groupList = numpy.arange(nChannels)
1419 dataOut.groupList = numpy.arange(nChannels)
1392
1420
1393 #Radial Velocities
1421 #Radial Velocities
1394 # phase = numpy.angle(data_acf[:,1,:,:])
1422 phase = numpy.angle(data_acf[:,1,:,:])
1395 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1423 # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1396 velRad = phase*lamb/(4*numpy.pi*tSamp)
1424 velRad = phase*lamb/(4*numpy.pi*tSamp)
1397
1425
1398 #Spectral width
1426 #Spectral width
1399 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1427 # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1400 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1428 # acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1429 acf1 = data_acf[:,1,:,:]
1430 acf2 = data_acf[:,2,:,:]
1401
1431
1402 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1432 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1403 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
1433 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
@@ -1409,7 +1439,7 class NonSpecularMeteorDetection(Operation):
1409 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1439 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1410
1440
1411 #Radial velocity
1441 #Radial velocity
1412 boolMet2 = numpy.abs(velRad) < 30
1442 boolMet2 = numpy.abs(velRad) < 20
1413 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1443 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1414
1444
1415 #Spectral Width
1445 #Spectral Width
@@ -1436,9 +1466,9 class NonSpecularMeteorDetection(Operation):
1436
1466
1437 # self.dataOut.data_param = data_int
1467 # self.dataOut.data_param = data_int
1438 if len(data_param) == 0:
1468 if len(data_param) == 0:
1439 self.dataOut.flagNoData = True
1469 dataOut.flagNoData = True
1440 else:
1470 else:
1441 self.dataOut.data_param = data_param
1471 dataOut.data_param = data_param
1442
1472
1443 def __erase_small(self, binArray, threshX, threshY):
1473 def __erase_small(self, binArray, threshX, threshY):
1444 labarray, numfeat = ndimage.measurements.label(binArray)
1474 labarray, numfeat = ndimage.measurements.label(binArray)
@@ -2245,10 +2275,10 class SMPhaseCalibration(Operation):
2245
2275
2246 pairi = pairs[i]
2276 pairi = pairs[i]
2247
2277
2248 phip3 = phases[:,pairi[1]]
2278 phip3 = phases[:,pairi[0]]
2249 d3 = d[pairi[1]]
2279 d3 = d[pairi[0]]
2250 phip2 = phases[:,pairi[0]]
2280 phip2 = phases[:,pairi[1]]
2251 d2 = d[pairi[0]]
2281 d2 = d[pairi[1]]
2252 #Calculating gamma
2282 #Calculating gamma
2253 # jdcos = alp1/(k*d1)
2283 # jdcos = alp1/(k*d1)
2254 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2284 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
@@ -2303,8 +2333,8 class SMPhaseCalibration(Operation):
2303 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2333 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2304 meteorOps = SMOperations()
2334 meteorOps = SMOperations()
2305 nchan = 4
2335 nchan = 4
2306 pairx = pairsList[0]
2336 pairx = pairsList[0] #x es 0
2307 pairy = pairsList[1]
2337 pairy = pairsList[1] #y es 1
2308 center_xangle = 0
2338 center_xangle = 0
2309 center_yangle = 0
2339 center_yangle = 0
2310 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
2340 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
@@ -2331,14 +2361,28 class SMPhaseCalibration(Operation):
2331 # Iterations looking for the offset
2361 # Iterations looking for the offset
2332 for iy in range(int(nstepsy)):
2362 for iy in range(int(nstepsy)):
2333 for ix in range(int(nstepsx)):
2363 for ix in range(int(nstepsx)):
2334 jph[pairy[1]] = alpha_y[iy]
2364 d3 = d[pairsList[1][0]]
2335 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2365 d2 = d[pairsList[1][1]]
2336
2366 d5 = d[pairsList[0][0]]
2337 jph[pairx[1]] = alpha_x[ix]
2367 d4 = d[pairsList[0][1]]
2338 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2368
2339
2369 alp2 = alpha_y[iy] #gamma 1
2370 alp4 = alpha_x[ix] #gamma 0
2371
2372 alp3 = -alp2*d3/d2 - gammas[1]
2373 alp5 = -alp4*d5/d4 - gammas[0]
2374 # jph[pairy[1]] = alpha_y[iy]
2375 # jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2376
2377 # jph[pairx[1]] = alpha_x[ix]
2378 # jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2379 jph[pairsList[0][1]] = alp4
2380 jph[pairsList[0][0]] = alp5
2381 jph[pairsList[1][0]] = alp3
2382 jph[pairsList[1][1]] = alp2
2340 jph_array[:,ix,iy] = jph
2383 jph_array[:,ix,iy] = jph
2341
2384 # d = [2.0,2.5,2.5,2.0]
2385 #falta chequear si va a leer bien los meteoros
2342 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2386 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2343 error = meteorsArray1[:,-1]
2387 error = meteorsArray1[:,-1]
2344 ind1 = numpy.where(error==0)[0]
2388 ind1 = numpy.where(error==0)[0]
@@ -2387,7 +2431,8 class SMPhaseCalibration(Operation):
2387 k = 2*numpy.pi/lamb
2431 k = 2*numpy.pi/lamb
2388 azimuth = 0
2432 azimuth = 0
2389 h = (hmin, hmax)
2433 h = (hmin, hmax)
2390 pairs = ((0,1),(2,3))
2434 # pairs = ((0,1),(2,3)) #Estrella
2435 # pairs = ((1,0),(2,3)) #T
2391
2436
2392 if channelPositions is None:
2437 if channelPositions is None:
2393 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2438 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
@@ -2395,6 +2440,17 class SMPhaseCalibration(Operation):
2395 meteorOps = SMOperations()
2440 meteorOps = SMOperations()
2396 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2441 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2397
2442
2443 #Checking correct order of pairs
2444 pairs = []
2445 if distances[1] > distances[0]:
2446 pairs.append((1,0))
2447 else:
2448 pairs.append((0,1))
2449
2450 if distances[3] > distances[2]:
2451 pairs.append((3,2))
2452 else:
2453 pairs.append((2,3))
2398 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2454 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2399
2455
2400 meteorsArray = self.__buffer
2456 meteorsArray = self.__buffer
General Comments 0
You need to be logged in to leave comments. Login now