| @@ -174,5 +174,5 class CorrelationProc(ProcessingUnit): | |||
|
|
174 | 174 | self.dataOut.lagRange = numpy.array(lags)*delta |
|
|
175 | 175 | # self.dataOut.nCohInt = self.dataIn.nCohInt*nAvg |
|
|
176 | 176 | self.dataOut.flagNoData = False |
|
|
177 | a = self.dataOut.normFactor | |
|
|
177 | # a = self.dataOut.normFactor | |
|
|
178 | 178 | return |
| @@ -1017,33 +1017,55 class WindProfiler(Operation): | |||
|
|
1017 | 1017 | def techniqueNSM_DBS(self, **kwargs): |
|
|
1018 | 1018 | metArray = kwargs['metArray'] |
|
|
1019 | 1019 | heightList = kwargs['heightList'] |
|
|
1020 | timeList = kwargs['timeList'] | |
|
|
1021 |
z |
|
|
|
1020 | timeList = kwargs['timeList'] | |
|
|
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 | 1032 | nChan = numpy.max(cmet) + 1 |
|
|
1023 | 1033 | nHeights = len(heightList) |
|
|
1024 | 1034 | |
|
|
1025 | utctime = metArray[:,0] | |
|
|
1026 |
|
|
|
|
1027 | hmet = metArray1[:,3].astype(int) | |
|
|
1028 | h1met = heightList[hmet]*zenithList[cmet] | |
|
|
1029 | vmet = metArray1[:,5] | |
|
|
1035 | azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) | |
|
|
1036 | hmet = heightList[hmet] | |
|
|
1037 | h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights | |
|
|
1038 | ||
|
|
1039 | velEst = numpy.zeros((heightList.size,2))*numpy.nan | |
|
|
1030 | 1040 | |
|
|
1031 | 1041 | for i in range(nHeights - 1): |
|
|
1032 | 1042 | hmin = heightList[i] |
|
|
1033 | 1043 | hmax = heightList[i + 1] |
|
|
1034 | 1044 | |
|
|
1035 |
|
|
|
|
1036 | ||
|
|
1037 | ||
|
|
1038 | ||
|
|
1039 | return data_output | |
|
|
1045 | thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10) | |
|
|
1046 | indthisH = numpy.where(thisH) | |
|
|
1047 | ||
|
|
1048 | if numpy.size(indthisH) > 3: | |
|
|
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, |
|
|
|
1063 | def run(self, dataOut, technique, **kwargs): | |
|
|
1042 | 1064 | |
|
|
1043 | 1065 | param = dataOut.data_param |
|
|
1044 | 1066 | if dataOut.abscissaList != None: |
|
|
1045 | 1067 | absc = dataOut.abscissaList[:-1] |
|
|
1046 |
|
|
|
|
1068 | noise = dataOut.noise | |
|
|
1047 | 1069 | heightList = dataOut.heightList |
|
|
1048 | 1070 | SNR = dataOut.data_SNR |
|
|
1049 | 1071 | |
| @@ -1153,11 +1175,15 class WindProfiler(Operation): | |||
|
|
1153 | 1175 | else: rx_location = [(0,1),(1,1),(1,0)] |
|
|
1154 | 1176 | if kwargs.has_key('azimuth'): |
|
|
1155 | 1177 | azimuth = kwargs['azimuth'] |
|
|
1156 | else: azimuth = 51 | |
|
|
1178 | else: azimuth = 51.06 | |
|
|
1157 | 1179 | if kwargs.has_key('dfactor'): |
|
|
1158 | 1180 | dfactor = kwargs['dfactor'] |
|
|
1159 | 1181 | if kwargs.has_key('mode'): |
|
|
1160 | 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 | 1187 | else: mode = 'SA' |
|
|
1162 | 1188 | |
|
|
1163 | 1189 | #Borrar luego esto |
| @@ -1200,7 +1226,7 class WindProfiler(Operation): | |||
|
|
1200 | 1226 | if mode == 'SA': |
|
|
1201 | 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 | 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 | 1230 | dataOut.data_output = dataOut.data_output.T |
|
|
1205 | 1231 | dataOut.flagNoData = False |
|
|
1206 | 1232 | self.__buffer = None |
| @@ -1277,25 +1303,26 class EWDriftsEstimation(Operation): | |||
|
|
1277 | 1303 | |
|
|
1278 | 1304 | class NonSpecularMeteorDetection(Operation): |
|
|
1279 | 1305 | |
|
|
1280 | def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False): | |
|
|
1281 |
data_acf = |
|
|
|
1282 |
data_ccf = |
|
|
|
1283 | ||
|
|
1284 | lamb = self.dataOut.C/self.dataOut.frequency | |
|
|
1285 | tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt | |
|
|
1286 | paramInterval = self.dataOut.paramInterval | |
|
|
1287 | ||
|
|
1306 | def run(self, dataOut, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False): | |
|
|
1307 | data_acf = dataOut.data_pre[0] | |
|
|
1308 | data_ccf = dataOut.data_pre[1] | |
|
|
1309 | pairsList = dataOut.groupList[1] | |
|
|
1310 | ||
|
|
1311 | lamb = dataOut.C/dataOut.frequency | |
|
|
1312 | tSamp = dataOut.ippSeconds*dataOut.nCohInt | |
|
|
1313 | paramInterval = dataOut.paramInterval | |
|
|
1314 | ||
|
|
1288 | 1315 | nChannels = data_acf.shape[0] |
|
|
1289 | 1316 | nLags = data_acf.shape[1] |
|
|
1290 | 1317 | nProfiles = data_acf.shape[2] |
|
|
1291 |
nHeights = |
|
|
|
1292 |
nCohInt = |
|
|
|
1293 |
sec = numpy.round(nProfiles/ |
|
|
|
1294 |
heightList = |
|
|
|
1295 |
ippSeconds = |
|
|
|
1296 |
utctime = |
|
|
|
1297 | ||
|
|
1298 |
|
|
|
|
1318 | nHeights = dataOut.nHeights | |
|
|
1319 | nCohInt = dataOut.nCohInt | |
|
|
1320 | sec = numpy.round(nProfiles/dataOut.paramInterval) | |
|
|
1321 | heightList = dataOut.heightList | |
|
|
1322 | ippSeconds = dataOut.ippSeconds*dataOut.nCohInt*dataOut.nAvg | |
|
|
1323 | utctime = dataOut.utctime | |
|
|
1324 | ||
|
|
1325 | dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds) | |
|
|
1299 | 1326 | |
|
|
1300 | 1327 | #------------------------ SNR -------------------------------------- |
|
|
1301 | 1328 | power = data_acf[:,0,:,:].real |
| @@ -1308,6 +1335,7 class NonSpecularMeteorDetection(Operation): | |||
|
|
1308 | 1335 | SNRdB = 10*numpy.log10(SNR) |
|
|
1309 | 1336 | |
|
|
1310 | 1337 | if mode == 'SA': |
|
|
1338 | dataOut.groupList = dataOut.groupList[1] | |
|
|
1311 | 1339 | nPairs = data_ccf.shape[0] |
|
|
1312 | 1340 | #---------------------- Coherence and Phase -------------------------- |
|
|
1313 | 1341 | phase = numpy.zeros(data_ccf[:,0,:,:].shape) |
| @@ -1315,8 +1343,8 class NonSpecularMeteorDetection(Operation): | |||
|
|
1315 | 1343 | coh1 = numpy.zeros(data_ccf[:,0,:,:].shape) |
|
|
1316 | 1344 | |
|
|
1317 | 1345 | for p in range(nPairs): |
|
|
1318 |
ch0 = |
|
|
|
1319 |
ch1 = |
|
|
|
1346 | ch0 = pairsList[p][0] | |
|
|
1347 | ch1 = pairsList[p][1] | |
|
|
1320 | 1348 | ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:]) |
|
|
1321 | 1349 | phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter |
|
|
1322 | 1350 | # phase1[p,:,:] = numpy.angle(ccf) #median filter |
| @@ -1388,16 +1416,18 class NonSpecularMeteorDetection(Operation): | |||
|
|
1388 | 1416 | data_param[:,6:] = phase[:,tmet,hmet].T |
|
|
1389 | 1417 | |
|
|
1390 | 1418 | elif mode == 'DBS': |
|
|
1391 |
|
|
|
|
1419 | dataOut.groupList = numpy.arange(nChannels) | |
|
|
1392 | 1420 | |
|
|
1393 | 1421 | #Radial Velocities |
|
|
1394 |
|
|
|
|
1395 | phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1)) | |
|
|
1422 | phase = numpy.angle(data_acf[:,1,:,:]) | |
|
|
1423 | # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1)) | |
|
|
1396 | 1424 | velRad = phase*lamb/(4*numpy.pi*tSamp) |
|
|
1397 | 1425 | |
|
|
1398 | 1426 | #Spectral width |
|
|
1399 | 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)) | |
|
|
1427 | # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), 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 | 1432 | spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2)) |
|
|
1403 | 1433 | # velRad = ndimage.median_filter(velRad, size = (1,5,1)) |
| @@ -1409,7 +1439,7 class NonSpecularMeteorDetection(Operation): | |||
|
|
1409 | 1439 | boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5)) |
|
|
1410 | 1440 | |
|
|
1411 | 1441 | #Radial velocity |
|
|
1412 |
boolMet2 = numpy.abs(velRad) < |
|
|
|
1442 | boolMet2 = numpy.abs(velRad) < 20 | |
|
|
1413 | 1443 | boolMet2 = ndimage.median_filter(boolMet2, (1,5,5)) |
|
|
1414 | 1444 | |
|
|
1415 | 1445 | #Spectral Width |
| @@ -1436,9 +1466,9 class NonSpecularMeteorDetection(Operation): | |||
|
|
1436 | 1466 | |
|
|
1437 | 1467 | # self.dataOut.data_param = data_int |
|
|
1438 | 1468 | if len(data_param) == 0: |
|
|
1439 |
|
|
|
|
1469 | dataOut.flagNoData = True | |
|
|
1440 | 1470 | else: |
|
|
1441 |
|
|
|
|
1471 | dataOut.data_param = data_param | |
|
|
1442 | 1472 | |
|
|
1443 | 1473 | def __erase_small(self, binArray, threshX, threshY): |
|
|
1444 | 1474 | labarray, numfeat = ndimage.measurements.label(binArray) |
| @@ -2245,10 +2275,10 class SMPhaseCalibration(Operation): | |||
|
|
2245 | 2275 | |
|
|
2246 | 2276 | pairi = pairs[i] |
|
|
2247 | 2277 | |
|
|
2248 |
phip3 = phases[:,pairi[ |
|
|
|
2249 |
d3 = d[pairi[ |
|
|
|
2250 |
phip2 = phases[:,pairi[ |
|
|
|
2251 |
d2 = d[pairi[ |
|
|
|
2278 | phip3 = phases[:,pairi[0]] | |
|
|
2279 | d3 = d[pairi[0]] | |
|
|
2280 | phip2 = phases[:,pairi[1]] | |
|
|
2281 | d2 = d[pairi[1]] | |
|
|
2252 | 2282 | #Calculating gamma |
|
|
2253 | 2283 | # jdcos = alp1/(k*d1) |
|
|
2254 | 2284 | # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0))) |
| @@ -2303,8 +2333,8 class SMPhaseCalibration(Operation): | |||
|
|
2303 | 2333 | def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray): |
|
|
2304 | 2334 | meteorOps = SMOperations() |
|
|
2305 | 2335 | nchan = 4 |
|
|
2306 | pairx = pairsList[0] | |
|
|
2307 | pairy = pairsList[1] | |
|
|
2336 | pairx = pairsList[0] #x es 0 | |
|
|
2337 | pairy = pairsList[1] #y es 1 | |
|
|
2308 | 2338 | center_xangle = 0 |
|
|
2309 | 2339 | center_yangle = 0 |
|
|
2310 | 2340 | range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4]) |
| @@ -2331,14 +2361,28 class SMPhaseCalibration(Operation): | |||
|
|
2331 | 2361 | # Iterations looking for the offset |
|
|
2332 | 2362 | for iy in range(int(nstepsy)): |
|
|
2333 | 2363 | for ix in range(int(nstepsx)): |
|
|
2334 |
|
|
|
|
2335 | jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]] | |
|
|
2336 | ||
|
|
2337 |
|
|
|
|
2338 | jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]] | |
|
|
2339 | ||
|
|
2364 | d3 = d[pairsList[1][0]] | |
|
|
2365 | d2 = d[pairsList[1][1]] | |
|
|
2366 | d5 = d[pairsList[0][0]] | |
|
|
2367 | d4 = d[pairsList[0][1]] | |
|
|
2368 | ||
|
|
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 | 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 | 2386 | meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph) |
|
|
2343 | 2387 | error = meteorsArray1[:,-1] |
|
|
2344 | 2388 | ind1 = numpy.where(error==0)[0] |
| @@ -2387,7 +2431,8 class SMPhaseCalibration(Operation): | |||
|
|
2387 | 2431 | k = 2*numpy.pi/lamb |
|
|
2388 | 2432 | azimuth = 0 |
|
|
2389 | 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 | 2437 | if channelPositions is None: |
|
|
2393 | 2438 | # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T |
| @@ -2395,6 +2440,17 class SMPhaseCalibration(Operation): | |||
|
|
2395 | 2440 | meteorOps = SMOperations() |
|
|
2396 | 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 | 2454 | # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb] |
|
|
2399 | 2455 | |
|
|
2400 | 2456 | meteorsArray = self.__buffer |
General Comments 0
You need to be logged in to leave comments.
Login now