##// END OF EJS Templates
schain
RSS
Rem Yagi profiles from amisr data
joabAM -
r1590:353012d3677c
parent child
Show More
Side by Side Unified
Context file:
r1590:353012d3677c Loading diff...:
Show file before | Show file after | Hide comments
@@ -16,6 +16,8 import numpy
16 import h5py
16 import h5py
17 import re
17 import re
18 import time
18 import time
19 from schainpy.utils import log
20
19 ################################################################################
21 ################################################################################
20 ################################################################################
22 ################################################################################
21 ################################################################################
23 ################################################################################
@@ -41,30 +43,33 def folder_in_range(folder, start_date, end_date, pattern):
41 ################################################################################
43 ################################################################################
42 ################################################################################
44 ################################################################################
43 def getHei_index( minHei, maxHei, heightList):
45 def getHei_index( minHei, maxHei, heightList):
44 if (minHei < heightList[0]):
45 minHei = heightList[0]
46
47 if (maxHei > heightList[-1]):
48 maxHei = heightList[-1]
49
50 minIndex = 0
51 maxIndex = 0
52 heights = numpy.asarray(heightList)
53
54 inda = numpy.where(heights >= minHei)
55 indb = numpy.where(heights <= maxHei)
56
57 try:
46 try:
58 minIndex = inda[0][0]
47 if (minHei < heightList[0]):
59 except:
48 minHei = heightList[0]
60 minIndex = 0
61
49
62 try:
50 if (maxHei > heightList[-1]):
63 maxIndex = indb[0][-1]
51 maxHei = heightList[-1]
64 except:
65 maxIndex = len(heightList)
66 return minIndex,maxIndex
67
52
53 minIndex = 0
54 maxIndex = 0
55 heights = numpy.asarray(heightList)
56
57 inda = numpy.where(heights >= minHei)
58 indb = numpy.where(heights <= maxHei)
59
60 try:
61 minIndex = inda[0][0]
62 except:
63 minIndex = 0
64
65 try:
66 maxIndex = indb[0][-1]
67 except:
68 maxIndex = len(heightList)
69 return minIndex,maxIndex
70 except Exception as e:
71 log.error("In getHei_index: ", __name__)
72 log.error(e , __name__)
68 ################################################################################
73 ################################################################################
69 ################################################################################
74 ################################################################################
70 ################################################################################
75 ################################################################################
Show file before | Show file after | Hide comments
@@ -11,6 +11,8 import numpy
11 #import copy
11 #import copy
12 from schainpy.model.data import _noise
12 from schainpy.model.data import _noise
13
13
14 from matplotlib import pyplot as plt
15
14 class VoltageProc(ProcessingUnit):
16 class VoltageProc(ProcessingUnit):
15
17
16 def __init__(self):
18 def __init__(self):
@@ -2115,6 +2117,8 class removeProfileByFaradayHS(Operation):
2115 # plt.grid()
2117 # plt.grid()
2116 # plt.show()
2118 # plt.show()
2117
2119
2120
2121
2118
2122
2119 outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
2123 outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
2120 outliers_IDs = numpy.unique(outliers_IDs)
2124 outliers_IDs = numpy.unique(outliers_IDs)
@@ -3232,11 +3236,32 class RemoveProfileSats2(Operation):
3232
3236
3233 Omite los perfiles contaminados con seΓ±al de satΓ©lites, usando una altura de referencia
3237 Omite los perfiles contaminados con seΓ±al de satΓ©lites, usando una altura de referencia
3234 promedia todas las alturas para los cΓ‘lculos
3238 promedia todas las alturas para los cΓ‘lculos
3235 In: minHei = min_sat_range
3239 In:
3236 max_sat_range
3240 n = Cantidad de perfiles que se acumularan, usualmente 10 segundos
3237 min_hei_ref
3241 navg = Porcentaje de perfiles que puede considerarse como satΓ©lite, mΓ‘ximo 90%
3238 max_hei_ref
3242 minHei =
3239 th = diference between profiles mean, ref and sats
3243 minRef =
3244 maxRef =
3245 nBins =
3246 profile_margin =
3247 th_hist_outlier =
3248 nProfilesOut =
3249
3250 Pensado para remover interferencias de las YAGI, se puede adaptar a otras interferencias
3251
3252 remYagi = Activa la funcion de remociΓ³n de interferencias de la YAGI
3253 nProfYagi = Cantidad de perfiles que son afectados, acorde NTX de la YAGI
3254 offYagi =
3255 minHJULIA = Altura mΓ­nima donde aparece la seΓ±al referencia de JULIA (-50)
3256 maxHJULIA = Altura mΓ‘xima donde aparece la seΓ±al referencia de JULIA (-15)
3257
3258 debug = Activa los grΓ‘ficos, recomendable ejecutar para ajustar los parΓ‘metros
3259 para un experimento en especΓ­fico.
3260
3261 ** se modifica para remover interferencias puntuales, es decir, desde otros radares.
3262 Inicialmente se ha configurado para omitir tambiΓ©n los perfiles de la YAGI en los datos
3263 de AMISR-ISR.
3264
3240 Out:
3265 Out:
3241 profile clean
3266 profile clean
3242 '''
3267 '''
@@ -3252,14 +3277,16 class RemoveProfileSats2(Operation):
3252
3277
3253 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
3278 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
3254 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
3279 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
3255 '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'thdB','prev_pnoise')
3280 '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'debug','prev_pnoise')
3256 def __init__(self, **kwargs):
3281 def __init__(self, **kwargs):
3257
3282
3258 Operation.__init__(self, **kwargs)
3283 Operation.__init__(self, **kwargs)
3259 self.isConfig = False
3284 self.isConfig = False
3285 self.currentTime = None
3260
3286
3261 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=15,
3287 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=15,minHei=None, maxHei=None, nBins=10,
3262 minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
3288 minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None,
3289 idate=None,startH=None,endH=None ):
3263
3290
3264 if n == None and timeInterval == None:
3291 if n == None and timeInterval == None:
3265 raise ValueError("nprofiles or timeInterval should be specified ...")
3292 raise ValueError("nprofiles or timeInterval should be specified ...")
@@ -3281,6 +3308,7 class RemoveProfileSats2(Operation):
3281 self.__profIndex = 0
3308 self.__profIndex = 0
3282 self.first_utcBlock = None
3309 self.first_utcBlock = None
3283 self.prev_pnoise = None
3310 self.prev_pnoise = None
3311 self.nBins = nBins
3284 #self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
3312 #self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
3285 minHei = minHei
3313 minHei = minHei
3286 maxHei = maxHei
3314 maxHei = maxHei
@@ -3293,7 +3321,26 class RemoveProfileSats2(Operation):
3293 self.nChannels = dataOut.nChannels
3321 self.nChannels = dataOut.nChannels
3294 self.nHeights = dataOut.nHeights
3322 self.nHeights = dataOut.nHeights
3295 self.test_counter = 0
3323 self.test_counter = 0
3296 self.thdB = thdB
3324 self.debug = debug
3325 self.remYagi = remYagi
3326
3327 if self.remYagi :
3328 if minHJULIA==None or maxHJULIA==None:
3329 raise ValueError("Parameters minHYagi and minHYagi are necessary!")
3330 return
3331 if idate==None or startH==None or endH==None:
3332 raise ValueError("Date and hour parameters are necessary!")
3333 return
3334 self.minHJULIA_idx,self.maxHJULIA_idx = getHei_index(minHJULIA, maxHJULIA, dataOut.heightList)
3335 self.offYagi = offYagi
3336 self.nTxYagi = nProfYagi
3337
3338 self.startTime = datetime.datetime.combine(idate,startH)
3339 self.endTime = datetime.datetime.combine(idate,endH)
3340
3341 log.warning("Be careful with the selection of parameters for sat removal! Is recommended to \
3342 activate the debug parameter in this operation for calibration", self.name)
3343
3297
3344
3298 def filterSatsProfiles(self):
3345 def filterSatsProfiles(self):
3299
3346
@@ -3302,10 +3349,22 class RemoveProfileSats2(Operation):
3302 nChannels, profiles, heights = data.shape
3349 nChannels, profiles, heights = data.shape
3303 indexes=numpy.zeros([], dtype=int)
3350 indexes=numpy.zeros([], dtype=int)
3304 indexes = numpy.delete(indexes,0)
3351 indexes = numpy.delete(indexes,0)
3352
3353 indexesYagi=numpy.zeros([], dtype=int)
3354 indexesYagi = numpy.delete(indexesYagi,0)
3355
3356 indexesYagi_up=numpy.zeros([], dtype=int)
3357 indexesYagi_up = numpy.delete(indexesYagi_up,0)
3358 indexesYagi_down=numpy.zeros([], dtype=int)
3359 indexesYagi_down = numpy.delete(indexesYagi_down,0)
3360
3361
3362 indexesJULIA=numpy.zeros([], dtype=int)
3363 indexesJULIA = numpy.delete(indexesJULIA,0)
3364
3305 outliers_IDs=[]
3365 outliers_IDs=[]
3306 #nBins = 16 # min Th = 12
3366
3307 nBins = 16
3367 div = profiles//self.nBins
3308 div = profiles//nBins
3309
3368
3310 for c in range(nChannels):
3369 for c in range(nChannels):
3311 #print(self.min_ref,self.max_ref)
3370 #print(self.min_ref,self.max_ref)
@@ -3331,7 +3390,7 class RemoveProfileSats2(Operation):
3331
3390
3332
3391
3333
3392
3334 power = (data[c,:,self.minHei_idx:] * numpy.conjugate(data[c,:,self.minHei_idx:])).real
3393 power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx])).real
3335 heis = len(power[0,:])
3394 heis = len(power[0,:])
3336 power = power.mean(axis=1)
3395 power = power.mean(axis=1)
3337
3396
@@ -3347,8 +3406,23 class RemoveProfileSats2(Operation):
3347 if index[0].size < int(self.navg*profiles):
3406 if index[0].size < int(self.navg*profiles):
3348 indexes = numpy.append(indexes, index[0])
3407 indexes = numpy.append(indexes, index[0])
3349
3408
3409 #print("sdas ", noise_ref.mean())
3410
3411 if self.remYagi :
3412 #print(self.minHJULIA_idx, self.maxHJULIA_idx)
3413 powerJULIA = (data[c,:,self.minHJULIA_idx:self.maxHJULIA_idx] * numpy.conjugate(data[c,:,self.minHJULIA_idx:self.maxHJULIA_idx])).real
3414 powerJULIA = powerJULIA.mean(axis=1)
3415 th_JULIA = powerJULIA.mean()*0.85
3416 indexJULIA = numpy.where(powerJULIA >= th_JULIA )
3417
3418 indexesJULIA= numpy.append(indexesJULIA, indexJULIA[0])
3419
3420 # fig, ax = plt.subplots()
3421 # ax.plot(powerJULIA)
3422 # ax.axhline(th_JULIA, color='r')
3423 # plt.grid()
3424 # plt.show()
3350
3425
3351 # from matplotlib import pyplot as plt
3352 # fig, ax = plt.subplots()
3426 # fig, ax = plt.subplots()
3353 # ax.plot(fpower)
3427 # ax.plot(fpower)
3354 # ax.axhline(th, color='r')
3428 # ax.axhline(th, color='r')
@@ -3360,6 +3434,35 class RemoveProfileSats2(Operation):
3360
3434
3361 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
3435 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
3362 #outliers_IDs = numpy.unique(outliers_IDs)
3436 #outliers_IDs = numpy.unique(outliers_IDs)
3437 # print(indexesJULIA)
3438 if len(indexesJULIA > 1):
3439 iJ = indexesJULIA
3440 locs = [ (iJ[n]-iJ[n-1]) > 5 for n in range(len(iJ))]
3441 locs_2 = numpy.where(locs)[0]
3442 #print(locs_2, indexesJULIA[locs_2-1])
3443 indexesYagi_up = numpy.append(indexesYagi_up, indexesJULIA[locs_2-1])
3444 indexesYagi_down = numpy.append(indexesYagi_down, indexesJULIA[locs_2])
3445
3446
3447 indexesYagi_up = numpy.append(indexesYagi_up,indexesJULIA[-1])
3448 indexesYagi_down = numpy.append(indexesYagi_down,indexesJULIA[0])
3449
3450 indexesYagi_up = numpy.unique(indexesYagi_up)
3451 indexesYagi_down = numpy.unique(indexesYagi_down)
3452
3453
3454 aux_ind = [ numpy.arange( (self.offYagi + k)+1, (self.offYagi + k + self.nTxYagi)+1, 1, dtype=int) for k in indexesYagi_up]
3455 indexesYagi_up = (numpy.asarray(aux_ind)).flatten()
3456
3457 aux_ind2 = [ numpy.arange( (k - self.nTxYagi)+1, k+1 , 1, dtype=int) for k in indexesYagi_down]
3458 indexesYagi_down = (numpy.asarray(aux_ind2)).flatten()
3459
3460 indexesYagi = numpy.append(indexesYagi,indexesYagi_up)
3461 indexesYagi = numpy.append(indexesYagi,indexesYagi_down)
3462
3463
3464 indexesYagi = indexesYagi[ (indexesYagi >= 0) & (indexesYagi<profiles)]
3465 indexesYagi = numpy.unique(indexesYagi)
3363
3466
3364 outs_lines = numpy.unique(indexes)
3467 outs_lines = numpy.unique(indexes)
3365
3468
@@ -3388,30 +3491,49 class RemoveProfileSats2(Operation):
3388 # print("pmax: ",ipmax)
3491 # print("pmax: ",ipmax)
3389
3492
3390
3493
3494
3391
3495
3392 #print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
3496 #print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
3393 outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
3497 outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
3394 #print("outliers final: ", outlier_loc_index)
3498 #print("outliers final: ", outlier_loc_index)
3395
3499
3396 # from matplotlib import pyplot as plt
3500
3397 # x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
3501 if self.debug:
3398 # fig, ax = plt.subplots(nChannels,2,figsize=(8, 6))
3502 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
3399 # for i in range(nChannels):
3503 fig, ax = plt.subplots(nChannels,2,figsize=(8, 6))
3400 # dat = data[i,:,:].real
3504
3401 # dat = 10* numpy.log10((data[i,:,:] * numpy.conjugate(data[i,:,:])).real)
3505 for i in range(nChannels):
3402 # m = numpy.nanmean(dat)
3506 dat = data[i,:,:].real
3403 # o = numpy.nanstd(dat)
3507 dat = 10* numpy.log10((data[i,:,:] * numpy.conjugate(data[i,:,:])).real)
3404 # #print(m, o, x.shape, y.shape)
3508 m = numpy.nanmean(dat)
3405 # #c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
3509 o = numpy.nanstd(dat)
3406 # c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
3510 if nChannels>1:
3407 # ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3511 c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
3408 # #fig.colorbar(c)
3512 ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3409 # ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
3513 #fig.colorbar(c)
3410 # ax[i][1].hist(outs_lines,bins=my_bins)
3514 ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
3411 # plt.show()
3515 ax[i][1].hist(outs_lines,bins=my_bins)
3516 if self.remYagi :
3517 ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m')
3518 else:
3519 c = ax[0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
3520 ax[0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3521 #fig.colorbar(c)
3522 ax[0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
3523
3524 ax[1].hist(outs_lines,bins=my_bins)
3525 if self.remYagi :
3526 ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m')
3527 plt.show()
3528
3529
3530
3412
3531
3532 if self.remYagi and (self.currentTime < self.startTime and self.currentTime < self.endTime):
3533 outlier_loc_index = numpy.append(outlier_loc_index,indexesYagi)
3413
3534
3414 self.outliers_IDs_list = outlier_loc_index
3535 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
3536
3415 #print("outs list: ", self.outliers_IDs_list)
3537 #print("outs list: ", self.outliers_IDs_list)
3416 return self.__buffer_data
3538 return self.__buffer_data
3417
3539
@@ -3467,16 +3589,20 class RemoveProfileSats2(Operation):
3467
3589
3468 return data
3590 return data
3469
3591
3470 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,
3592 def run(self, dataOut, n=None, navg=0.9, nProfilesOut=1, profile_margin=50, th_hist_outlier=15,minHei=None,nBins=10,
3471 th_hist_outlier=15,minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
3593 maxHei=None, minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None,
3594 idate=None,startH=None,endH=None):
3472
3595
3473 if not self.isConfig:
3596 if not self.isConfig:
3474 #print("init p idx: ", dataOut.profileIndex )
3597 #print("init p idx: ", dataOut.profileIndex )
3475 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,
3598 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,minHei=minHei,
3476 minHei=minHei, maxHei=maxHei, minRef=minRef, maxRef=maxRef, thdB=thdB)
3599 nBins=10, maxHei=maxHei, minRef=minRef, maxRef=maxRef, debug=debug, remYagi=remYagi, nProfYagi = nProfYagi,
3600 offYagi=offYagi, minHJULIA=minHJULIA,maxHJULIA=maxHJULIA,idate=idate,startH=startH,endH=endH)
3601
3477 self.isConfig = True
3602 self.isConfig = True
3478
3603
3479 dataBlock = None
3604 dataBlock = None
3605 self.currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
3480
3606
3481 if not dataOut.buffer_empty: #hay datos acumulados
3607 if not dataOut.buffer_empty: #hay datos acumulados
3482
3608
General Comments 0
You need to be logged in to leave comments. Login now