##// END OF EJS Templates
schain
RSS
HDFWriter writes standard Weather parameters/Now we can run a type OTHER Operation even if flagNoData is True
rflores -
r1459:52714d4d854d
parent child
Show More
Side by Side Unified
Context file:
r1459:52714d4d854d Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -1,834 +1,838
1 1 '''
2 2 Created on Jul 3, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 # SUBCHANNELS EN VEZ DE CHANNELS
7 7 # BENCHMARKS -> PROBLEMAS CON ARCHIVOS GRANDES -> INCONSTANTE EN EL TIEMPO
8 8 # ACTUALIZACION DE VERSION
9 9 # HEADERS
10 10 # MODULO DE ESCRITURA
11 11 # METADATA
12 12
13 13 import os
14 14 import time
15 15 import datetime
16 16 import numpy
17 17 import timeit
18 18 from fractions import Fraction
19 19 from time import time
20 20 from time import sleep
21 21
22 22 import schainpy.admin
23 23 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
24 24 from schainpy.model.data.jrodata import Voltage
25 25 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
26 26
27 27 import pickle
28 28 try:
29 29 import digital_rf
30 30 except:
31 31 pass
32 32
33 33
34 34 class DigitalRFReader(ProcessingUnit):
35 35 '''
36 36 classdocs
37 37 '''
38 38
39 39 def __init__(self):
40 40 '''
41 41 Constructor
42 42 '''
43 43
44 44 ProcessingUnit.__init__(self)
45 45
46 46 self.dataOut = Voltage()
47 47 self.__printInfo = True
48 48 self.__flagDiscontinuousBlock = False
49 49 self.__bufferIndex = 9999999
50 50 self.__codeType = 0
51 51 self.__ippKm = None
52 52 self.__nCode = None
53 53 self.__nBaud = None
54 54 self.__code = None
55 55 self.dtype = None
56 56 self.oldAverage = None
57 57 self.path = None
58 58
59 59 def close(self):
60 60 print('Average of writing to digital rf format is ', self.oldAverage * 1000)
61 61 return
62 62
63 63 def __getCurrentSecond(self):
64 64
65 65 return self.__thisUnixSample / self.__sample_rate
66 66
67 67 thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.")
68 68
69 69 def __setFileHeader(self):
70 70 '''
71 71 In this method will be initialized every parameter of dataOut object (header, no data)
72 72 '''
73 73 ippSeconds = 1.0 * self.__nSamples / self.__sample_rate
74 74 if not self.getByBlock:
75 75 nProfiles = 1.0 / ippSeconds # Number of profiles in one second
76 76 else:
77 77 nProfiles = self.nProfileBlocks # Number of profiles in one block
78 78
79 79 try:
80 80 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
81 81 self.__radarControllerHeader)
82 82 except:
83 83 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
84 84 txA=0,
85 85 txB=0,
86 86 nWindows=1,
87 87 nHeights=self.__nSamples,
88 88 firstHeight=self.__firstHeigth,
89 89 deltaHeight=self.__deltaHeigth,
90 90 codeType=self.__codeType,
91 91 nCode=self.__nCode, nBaud=self.__nBaud,
92 92 code=self.__code)
93 93
94 94 try:
95 95 self.dataOut.systemHeaderObj = SystemHeader(self.__systemHeader)
96 96 except:
97 97 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
98 98 nProfiles=nProfiles,
99 99 nChannels=len(
100 100 self.__channelList),
101 101 adcResolution=14)
102 102 self.dataOut.type = "Voltage"
103 103
104 104 self.dataOut.data = None
105 105
106 106 self.dataOut.dtype = self.dtype
107 107
108 108 # self.dataOut.nChannels = 0
109 109
110 110 # self.dataOut.nHeights = 0
111 111
112 112 self.dataOut.nProfiles = int(nProfiles)
113 113
114 114 self.dataOut.heightList = self.__firstHeigth + \
115 115 numpy.arange(self.__nSamples, dtype=numpy.float) * \
116 116 self.__deltaHeigth
117 117
118 118 #self.dataOut.channelList = list(range(self.__num_subchannels))
119 119 self.dataOut.channelList = list(range(len(self.__channelList)))
120 120 if not self.getByBlock:
121 121
122 122 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
123 123 else:
124 124 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights*self.nProfileBlocks
125 125
126 126 # self.dataOut.channelIndexList = None
127 127
128 128 self.dataOut.flagNoData = True
129 129 if not self.getByBlock:
130 130 self.dataOut.flagDataAsBlock = False
131 131 else:
132 132 self.dataOut.flagDataAsBlock = True
133 133 # Set to TRUE if the data is discontinuous
134 134 self.dataOut.flagDiscontinuousBlock = False
135 135
136 136 self.dataOut.utctime = None
137 137
138 138 # timezone like jroheader, difference in minutes between UTC and localtime
139 139 self.dataOut.timeZone = self.__timezone / 60
140 140
141 141 self.dataOut.dstFlag = 0
142 142
143 143 self.dataOut.errorCount = 0
144 144
145 145 try:
146 146 self.dataOut.nCohInt = self.fixed_metadata_dict.get(
147 147 'nCohInt', self.nCohInt)
148 148
149 149 # asumo que la data esta decodificada
150 150 self.dataOut.flagDecodeData = self.fixed_metadata_dict.get(
151 151 'flagDecodeData', self.flagDecodeData)
152 152
153 153 # asumo que la data esta sin flip
154 154 self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData']
155 155
156 156 self.dataOut.flagShiftFFT = self.fixed_metadata_dict['flagShiftFFT']
157 157
158 158 self.dataOut.useLocalTime = self.fixed_metadata_dict['useLocalTime']
159 159 except:
160 160 pass
161 161
162 162 self.dataOut.ippSeconds = ippSeconds
163 163
164 164 # Time interval between profiles
165 165 # self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
166 166
167 167 self.dataOut.frequency = self.__frequency
168 168
169 169 self.dataOut.realtime = self.__online
170 170
171 171 def findDatafiles(self, path, startDate=None, endDate=None):
172 172
173 173 if not os.path.isdir(path):
174 174 return []
175 175
176 176 try:
177 177 digitalReadObj = digital_rf.DigitalRFReader(
178 178 path, load_all_metadata=True)
179 179 except:
180 180 digitalReadObj = digital_rf.DigitalRFReader(path)
181 181
182 182 channelNameList = digitalReadObj.get_channels()
183 183
184 184 if not channelNameList:
185 185 return []
186 186
187 187 metadata_dict = digitalReadObj.get_rf_file_metadata(channelNameList[0])
188 188
189 189 sample_rate = metadata_dict['sample_rate'][0]
190 190
191 191 this_metadata_file = digitalReadObj.get_metadata(channelNameList[0])
192 192
193 193 try:
194 194 timezone = this_metadata_file['timezone'].value
195 195 except:
196 196 timezone = 0
197 197
198 198 startUTCSecond, endUTCSecond = digitalReadObj.get_bounds(
199 199 channelNameList[0]) / sample_rate - timezone
200 200
201 201 startDatetime = datetime.datetime.utcfromtimestamp(startUTCSecond)
202 202 endDatatime = datetime.datetime.utcfromtimestamp(endUTCSecond)
203 203
204 204 if not startDate:
205 205 startDate = startDatetime.date()
206 206
207 207 if not endDate:
208 208 endDate = endDatatime.date()
209 209
210 210 dateList = []
211 211
212 212 thisDatetime = startDatetime
213 213
214 214 while(thisDatetime <= endDatatime):
215 215
216 216 thisDate = thisDatetime.date()
217 217
218 218 if thisDate < startDate:
219 219 continue
220 220
221 221 if thisDate > endDate:
222 222 break
223 223
224 224 dateList.append(thisDate)
225 225 thisDatetime += datetime.timedelta(1)
226 226
227 227 return dateList
228 228
229 229 def setup(self, path=None,
230 230 startDate=None,
231 231 endDate=None,
232 232 startTime=datetime.time(0, 0, 0),
233 233 endTime=datetime.time(23, 59, 59),
234 234 channelList=None,
235 235 nSamples=None,
236 236 online=False,
237 237 delay=60,
238 238 buffer_size=1024,
239 239 ippKm=None,
240 240 nCohInt=1,
241 241 nCode=1,
242 242 nBaud=1,
243 243 flagDecodeData=False,
244 244 code=numpy.ones((1, 1), dtype=numpy.int),
245 245 getByBlock=0,
246 246 nProfileBlocks=1,
247 247 **kwargs):
248 248 '''
249 249 In this method we should set all initial parameters.
250 250
251 251 Inputs:
252 252 path
253 253 startDate
254 254 endDate
255 255 startTime
256 256 endTime
257 257 set
258 258 expLabel
259 259 ext
260 260 online
261 261 delay
262 262 '''
263 263 self.path = path
264 264 self.nCohInt = nCohInt
265 265 self.flagDecodeData = flagDecodeData
266 266 self.i = 0
267 267
268 268 self.getByBlock = getByBlock
269 269 self.nProfileBlocks = nProfileBlocks
270 270 if not os.path.isdir(path):
271 271 raise ValueError("[Reading] Directory %s does not exist" % path)
272 272
273 273 try:
274 274 self.digitalReadObj = digital_rf.DigitalRFReader(
275 275 path, load_all_metadata=True)
276 276 except:
277 277 self.digitalReadObj = digital_rf.DigitalRFReader(path)
278 278
279 279 channelNameList = self.digitalReadObj.get_channels()
280 280
281 281 if not channelNameList:
282 282 raise ValueError("[Reading] Directory %s does not have any files" % path)
283 283
284 284 if not channelList:
285 285 channelList = list(range(len(channelNameList)))
286 286
287 287 ########## Reading metadata ######################
288 288
289 289 top_properties = self.digitalReadObj.get_properties(
290 290 channelNameList[channelList[0]])
291 291
292 292 self.__num_subchannels = top_properties['num_subchannels']
293 293 self.__sample_rate = 1.0 * \
294 294 top_properties['sample_rate_numerator'] / \
295 295 top_properties['sample_rate_denominator']
296 296 # self.__samples_per_file = top_properties['samples_per_file'][0]
297 297 self.__deltaHeigth = 1e6 * 0.15 / self.__sample_rate # why 0.15?
298 298
299 299 this_metadata_file = self.digitalReadObj.get_digital_metadata(
300 300 channelNameList[channelList[0]])
301 301 metadata_bounds = this_metadata_file.get_bounds()
302 302 self.fixed_metadata_dict = this_metadata_file.read(
303 303 metadata_bounds[0])[metadata_bounds[0]] # GET FIRST HEADER
304 304
305 305 try:
306 306 self.__processingHeader = self.fixed_metadata_dict['processingHeader']
307 307 self.__radarControllerHeader = self.fixed_metadata_dict['radarControllerHeader']
308 308 self.__systemHeader = self.fixed_metadata_dict['systemHeader']
309 309 self.dtype = pickle.loads(self.fixed_metadata_dict['dtype'])
310 310 except:
311 311 pass
312 312
313 313 self.__frequency = None
314 314
315 315 self.__frequency = self.fixed_metadata_dict.get('frequency', 1)
316 316
317 317 self.__timezone = self.fixed_metadata_dict.get('timezone', 18000)
318 318
319 319 try:
320 320 nSamples = self.fixed_metadata_dict['nSamples']
321 321 except:
322 322 nSamples = None
323 323
324 324 self.__firstHeigth = 0
325 325
326 326 try:
327 327 codeType = self.__radarControllerHeader['codeType']
328 328 except:
329 329 codeType = 0
330 330
331 331 try:
332 332 if codeType:
333 333 nCode = self.__radarControllerHeader['nCode']
334 334 nBaud = self.__radarControllerHeader['nBaud']
335 335 code = self.__radarControllerHeader['code']
336 336 except:
337 337 pass
338 338
339 339 if not ippKm:
340 340 try:
341 341 # seconds to km
342 342 ippKm = self.__radarControllerHeader['ipp']
343 343 except:
344 344 ippKm = None
345 345 ####################################################
346 346 self.__ippKm = ippKm
347 347 startUTCSecond = None
348 348 endUTCSecond = None
349 349
350 350 if startDate:
351 351 startDatetime = datetime.datetime.combine(startDate, startTime)
352 352 startUTCSecond = (
353 353 startDatetime - datetime.datetime(1970, 1, 1)).total_seconds() + self.__timezone
354 354
355 355 if endDate:
356 356 endDatetime = datetime.datetime.combine(endDate, endTime)
357 357 endUTCSecond = (endDatetime - datetime.datetime(1970,
358 358 1, 1)).total_seconds() + self.__timezone
359 359
360 360
361 361 #print(startUTCSecond,endUTCSecond)
362 362 start_index, end_index = self.digitalReadObj.get_bounds(
363 363 channelNameList[channelList[0]])
364 364
365 365 #print("*****",start_index,end_index)
366 366 if not startUTCSecond:
367 367 startUTCSecond = start_index / self.__sample_rate
368 368
369 369 if start_index > startUTCSecond * self.__sample_rate:
370 370 startUTCSecond = start_index / self.__sample_rate
371 371
372 372 if not endUTCSecond:
373 373 endUTCSecond = end_index / self.__sample_rate
374 374 if end_index < endUTCSecond * self.__sample_rate:
375 375 endUTCSecond = end_index / self.__sample_rate #Check UTC and LT time
376 376 if not nSamples:
377 377 if not ippKm:
378 378 raise ValueError("[Reading] nSamples or ippKm should be defined")
379 379 nSamples = int(ippKm / (1e6 * 0.15 / self.__sample_rate))
380 380
381 381 channelBoundList = []
382 382 channelNameListFiltered = []
383 383
384 384 for thisIndexChannel in channelList:
385 385 thisChannelName = channelNameList[thisIndexChannel]
386 386 start_index, end_index = self.digitalReadObj.get_bounds(
387 387 thisChannelName)
388 388 channelBoundList.append((start_index, end_index))
389 389 channelNameListFiltered.append(thisChannelName)
390 390
391 391 self.profileIndex = 0
392 392 self.i = 0
393 393 self.__delay = delay
394 394
395 395 self.__codeType = codeType
396 396 self.__nCode = nCode
397 397 self.__nBaud = nBaud
398 398 self.__code = code
399 399
400 400 self.__datapath = path
401 401 self.__online = online
402 402 self.__channelList = channelList
403 403 self.__channelNameList = channelNameListFiltered
404 404 self.__channelBoundList = channelBoundList
405 405 self.__nSamples = nSamples
406 406 if self.getByBlock:
407 407 nSamples = nSamples*nProfileBlocks
408 408
409 409
410 410 self.__samples_to_read = int(nSamples) # FIJO: AHORA 40
411 #self.__samples_to_read = int(1000000) # FIJO: AHORA 40
411 412 self.__nChannels = len(self.__channelList)
412 413 #print("------------------------------------------")
413 414 #print("self.__samples_to_read",self.__samples_to_read)
414 415 #print("self.__nSamples",self.__nSamples)
415 416 # son iguales y el buffer_index da 0
416 417 self.__startUTCSecond = startUTCSecond
417 418 self.__endUTCSecond = endUTCSecond
418 419
419 420 self.__timeInterval = 1.0 * self.__samples_to_read / \
420 421 self.__sample_rate # Time interval
421 422
422 423 if online:
423 424 # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read)
424 425 startUTCSecond = numpy.floor(endUTCSecond)
425 426
426 427 # por que en el otro metodo lo primero q se hace es sumar samplestoread
427 428 self.__thisUnixSample = int(startUTCSecond * self.__sample_rate) - self.__samples_to_read
428 429
429 430 #self.__data_buffer = numpy.zeros(
430 431 # (self.__num_subchannels, self.__samples_to_read), dtype=numpy.complex)
431 432 self.__data_buffer = numpy.zeros((int(len(channelList)), self.__samples_to_read), dtype=numpy.complex)
432 433
433 434
434 435 self.__setFileHeader()
435 436 self.isConfig = True
436 437
437 438 print("[Reading] Digital RF Data was found from %s to %s " % (
438 439 datetime.datetime.utcfromtimestamp(
439 440 self.__startUTCSecond - self.__timezone),
440 441 datetime.datetime.utcfromtimestamp(
441 442 self.__endUTCSecond - self.__timezone)
442 443 ))
443 444
444 445 print("[Reading] Starting process from %s to %s" % (datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
445 446 datetime.datetime.utcfromtimestamp(
446 447 endUTCSecond - self.__timezone)
447 448 ))
448 449 self.oldAverage = None
449 450 self.count = 0
450 451 self.executionTime = 0
451 452
452 453 def __reload(self):
453 454 # print
454 455 # print "%s not in range [%s, %s]" %(
455 456 # datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
456 457 # datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
457 458 # datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
458 459 # )
459 460 print("[Reading] reloading metadata ...")
460 461
461 462 try:
462 463 self.digitalReadObj.reload(complete_update=True)
463 464 except:
464 465 self.digitalReadObj = digital_rf.DigitalRFReader(self.path)
465 466
466 467 start_index, end_index = self.digitalReadObj.get_bounds(
467 468 self.__channelNameList[self.__channelList[0]])
468 469
469 470 if start_index > self.__startUTCSecond * self.__sample_rate:
470 471 self.__startUTCSecond = 1.0 * start_index / self.__sample_rate
471 472
472 473 if end_index > self.__endUTCSecond * self.__sample_rate:
473 474 self.__endUTCSecond = 1.0 * end_index / self.__sample_rate
474 475 print()
475 476 print("[Reading] New timerange found [%s, %s] " % (
476 477 datetime.datetime.utcfromtimestamp(
477 478 self.__startUTCSecond - self.__timezone),
478 479 datetime.datetime.utcfromtimestamp(
479 480 self.__endUTCSecond - self.__timezone)
480 481 ))
481 482
482 483 return True
483 484
484 485 return False
485 486
486 487 def timeit(self, toExecute):
487 488 t0 = time.time()
488 489 toExecute()
489 490 self.executionTime = time.time() - t0
490 491 if self.oldAverage is None:
491 492 self.oldAverage = self.executionTime
492 493 self.oldAverage = (self.executionTime + self.count *
493 494 self.oldAverage) / (self.count + 1.0)
494 495 self.count = self.count + 1.0
495 496 return
496 497
497 498 def __readNextBlock(self, seconds=30, volt_scale=1):
498 499 '''
499 500 '''
500 501
501 502 # Set the next data
502 503 self.__flagDiscontinuousBlock = False
503 504 self.__thisUnixSample += self.__samples_to_read
504 505
505 506 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
506 507 print ("[Reading] There are no more data into selected time-range")
507 508 if self.__online:
508 509 sleep(3)
509 510 self.__reload()
510 511 else:
511 512 return False
512 513
513 514 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
514 515 return False
515 516 self.__thisUnixSample -= self.__samples_to_read
516 517
517 518 indexChannel = 0
518 519
519 520 dataOk = False
520 521
521 522 for thisChannelName in self.__channelNameList: # TODO VARIOS CHANNELS?
522 523 for indexSubchannel in range(self.__num_subchannels):
523 524 try:
524 525 t0 = time()
526 #print("Unitindex",self.__thisUnixSample)
527 #print("__samples_to_read",self.__samples_to_read)
525 528 result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,
526 529 self.__samples_to_read,
527 530 thisChannelName, sub_channel=indexSubchannel)
528 531 self.executionTime = time() - t0
529 532 if self.oldAverage is None:
530 533 self.oldAverage = self.executionTime
531 534 self.oldAverage = (
532 535 self.executionTime + self.count * self.oldAverage) / (self.count + 1.0)
533 536 self.count = self.count + 1.0
534 537
535 538 except IOError as e:
536 539 # read next profile
537 540 self.__flagDiscontinuousBlock = True
538 541 print("[Reading] %s" % datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e)
539 542 break
540 543
541 544 if result.shape[0] != self.__samples_to_read:
542 545 self.__flagDiscontinuousBlock = True
543 546 print("[Reading] %s: Too few samples were found, just %d/%d samples" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
544 547 result.shape[0],
545 548 self.__samples_to_read))
546 549 break
547 550
548 551 self.__data_buffer[indexChannel, :] = result * volt_scale
549 552 indexChannel+=1
550 553
551 554 dataOk = True
552 555
553 556 self.__utctime = self.__thisUnixSample / self.__sample_rate
554 557
555 558 if not dataOk:
556 559 return False
557 560
558 561 print("[Reading] %s: %d samples <> %f sec" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
559 562 self.__samples_to_read,
560 563 self.__timeInterval))
561 564
562 565 self.__bufferIndex = 0
563 566
564 567 return True
565 568
566 569 def __isBufferEmpty(self):
567 570
568 571 return self.__bufferIndex > self.__samples_to_read - self.__nSamples # 40960 - 40
569 572
570 573 def getData(self, seconds=30, nTries=5):
571 574 '''
572 575 This method gets the data from files and put the data into the dataOut object
573 576
574 577 In addition, increase el the buffer counter in one.
575 578
576 579 Return:
577 580 data : retorna un perfil de voltages (alturas * canales) copiados desde el
578 581 buffer. Si no hay mas archivos a leer retorna None.
579 582
580 583 Affected:
581 584 self.dataOut
582 585 self.profileIndex
583 586 self.flagDiscontinuousBlock
584 587 self.flagIsNewBlock
585 588 '''
586 589 #print("getdata")
587 590 err_counter = 0
588 591 self.dataOut.flagNoData = True
589 592
590 593
591 594 if self.__isBufferEmpty():
592 595 #print("hi")
593 596 self.__flagDiscontinuousBlock = False
594 597
595 598 while True:
596 599 if self.__readNextBlock():
597 600 break
598 601 if self.__thisUnixSample > self.__endUTCSecond * self.__sample_rate:
599 602 raise schainpy.admin.SchainError('Error')
600 603 return
601 604
602 605 if self.__flagDiscontinuousBlock:
603 606 raise schainpy.admin.SchainError('discontinuous block found')
604 607 return
605 608
606 609 if not self.__online:
607 610 raise schainpy.admin.SchainError('Online?')
608 611 return
609 612
610 613 err_counter += 1
611 614 if err_counter > nTries:
612 615 raise schainpy.admin.SchainError('Max retrys reach')
613 616 return
614 617
615 618 print('[Reading] waiting %d seconds to read a new block' % seconds)
616 619 sleep(seconds)
617 620
618 621
619 622 if not self.getByBlock:
620 623
621 624 #print("self.__bufferIndex",self.__bufferIndex)# este valor siempre es cero aparentemente
622 625 self.dataOut.data = self.__data_buffer[:, self.__bufferIndex:self.__bufferIndex + self.__nSamples]
623 626 self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
624 627 self.dataOut.flagNoData = False
625 628 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
626 629 self.dataOut.profileIndex = self.profileIndex
627 630
628 631 self.__bufferIndex += self.__nSamples
629 632 self.profileIndex += 1
630 633
631 634 if self.profileIndex == self.dataOut.nProfiles:
632 635 self.profileIndex = 0
633 636 else:
634 637 # ojo debo anadir el readNextBLock y el __isBufferEmpty(
635 638 self.dataOut.flagNoData = False
636 639 buffer = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex + self.__samples_to_read]
640 #print("test",self.__bufferIndex)
637 641 buffer = buffer.reshape((self.__nChannels, self.nProfileBlocks, int(self.__samples_to_read/self.nProfileBlocks)))
638 642 self.dataOut.nProfileBlocks = self.nProfileBlocks
639 643 self.dataOut.data = buffer
640 644 self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
641 645 self.profileIndex += self.__samples_to_read
642 646 self.__bufferIndex += self.__samples_to_read
643 647 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
644 648 return True
645 649
646 650
647 651 def printInfo(self):
648 652 '''
649 653 '''
650 654 if self.__printInfo == False:
651 655 return
652 656
653 657 # self.systemHeaderObj.printInfo()
654 658 # self.radarControllerHeaderObj.printInfo()
655 659
656 660 self.__printInfo = False
657 661
658 662 def printNumberOfBlock(self):
659 663 '''
660 664 '''
661 665 return
662 666 # print self.profileIndex
663 667
664 668 def run(self, **kwargs):
665 669 '''
666 670 This method will be called many times so here you should put all your code
667 671 '''
668 672
669 673 if not self.isConfig:
670 674 self.setup(**kwargs)
671 675
672 676 self.getData(seconds=self.__delay)
673 677
674 678 return
675 679
676 680 @MPDecorator
677 681 class DigitalRFWriter(Operation):
678 682 '''
679 683 classdocs
680 684 '''
681 685
682 686 def __init__(self, **kwargs):
683 687 '''
684 688 Constructor
685 689 '''
686 690 Operation.__init__(self, **kwargs)
687 691 self.metadata_dict = {}
688 692 self.dataOut = None
689 693 self.dtype = None
690 694 self.oldAverage = 0
691 695
692 696 def setHeader(self):
693 697
694 698 self.metadata_dict['frequency'] = self.dataOut.frequency
695 699 self.metadata_dict['timezone'] = self.dataOut.timeZone
696 700 self.metadata_dict['dtype'] = pickle.dumps(self.dataOut.dtype)
697 701 self.metadata_dict['nProfiles'] = self.dataOut.nProfiles
698 702 self.metadata_dict['heightList'] = self.dataOut.heightList
699 703 self.metadata_dict['channelList'] = self.dataOut.channelList
700 704 self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData
701 705 self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData
702 706 self.metadata_dict['flagShiftFFT'] = self.dataOut.flagShiftFFT
703 707 self.metadata_dict['useLocalTime'] = self.dataOut.useLocalTime
704 708 self.metadata_dict['nCohInt'] = self.dataOut.nCohInt
705 709 self.metadata_dict['type'] = self.dataOut.type
706 710 self.metadata_dict['flagDataAsBlock']= getattr(
707 711 self.dataOut, 'flagDataAsBlock', None) # chequear
708 712
709 713 def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'):
710 714 '''
711 715 In this method we should set all initial parameters.
712 716 Input:
713 717 dataOut: Input data will also be outputa data
714 718 '''
715 719 self.setHeader()
716 720 self.__ippSeconds = dataOut.ippSeconds
717 721 self.__deltaH = dataOut.getDeltaH()
718 722 self.__sample_rate = 1e6 * 0.15 / self.__deltaH
719 723 self.__dtype = dataOut.dtype
720 724 if len(dataOut.dtype) == 2:
721 725 self.__dtype = dataOut.dtype[0]
722 726 self.__nSamples = dataOut.systemHeaderObj.nSamples
723 727 self.__nProfiles = dataOut.nProfiles
724 728
725 729 if self.dataOut.type != 'Voltage':
726 730 raise 'Digital RF cannot be used with this data type'
727 731 self.arr_data = numpy.ones((1, dataOut.nFFTPoints * len(
728 732 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
729 733 else:
730 734 self.arr_data = numpy.ones((self.__nSamples, len(
731 735 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
732 736
733 737 file_cadence_millisecs = 1000
734 738
735 739 sample_rate_fraction = Fraction(self.__sample_rate).limit_denominator()
736 740 sample_rate_numerator = int(sample_rate_fraction.numerator)
737 741 sample_rate_denominator = int(sample_rate_fraction.denominator)
738 742 start_global_index = dataOut.utctime * self.__sample_rate
739 743
740 744 uuid = 'prueba'
741 745 compression_level = 0
742 746 checksum = False
743 747 is_complex = True
744 748 num_subchannels = len(dataOut.channelList)
745 749 is_continuous = True
746 750 marching_periods = False
747 751
748 752 self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence,
749 753 fileCadence, start_global_index,
750 754 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
751 755 is_complex, num_subchannels, is_continuous, marching_periods)
752 756 metadata_dir = os.path.join(path, 'metadata')
753 757 os.system('mkdir %s' % (metadata_dir))
754 758 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000
755 759 sample_rate_numerator, sample_rate_denominator,
756 760 metadataFile)
757 761 self.isConfig = True
758 762 self.currentSample = 0
759 763 self.oldAverage = 0
760 764 self.count = 0
761 765 return
762 766
763 767 def writeMetadata(self):
764 768 start_idx = self.__sample_rate * self.dataOut.utctime
765 769
766 770 self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict(
767 771 )
768 772 self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict(
769 773 )
770 774 self.metadata_dict['systemHeader'] = self.dataOut.systemHeaderObj.getAsDict(
771 775 )
772 776 self.digitalMetadataWriteObj.write(start_idx, self.metadata_dict)
773 777 return
774 778
775 779 def timeit(self, toExecute):
776 780 t0 = time()
777 781 toExecute()
778 782 self.executionTime = time() - t0
779 783 if self.oldAverage is None:
780 784 self.oldAverage = self.executionTime
781 785 self.oldAverage = (self.executionTime + self.count *
782 786 self.oldAverage) / (self.count + 1.0)
783 787 self.count = self.count + 1.0
784 788 return
785 789
786 790 def writeData(self):
787 791 if self.dataOut.type != 'Voltage':
788 792 raise 'Digital RF cannot be used with this data type'
789 793 for channel in self.dataOut.channelList:
790 794 for i in range(self.dataOut.nFFTPoints):
791 795 self.arr_data[1][channel * self.dataOut.nFFTPoints +
792 796 i]['r'] = self.dataOut.data[channel][i].real
793 797 self.arr_data[1][channel * self.dataOut.nFFTPoints +
794 798 i]['i'] = self.dataOut.data[channel][i].imag
795 799 else:
796 800 for i in range(self.dataOut.systemHeaderObj.nSamples):
797 801 for channel in self.dataOut.channelList:
798 802 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
799 803 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
800 804
801 805 def f(): return self.digitalWriteObj.rf_write(self.arr_data)
802 806 self.timeit(f)
803 807
804 808 return
805 809
806 810 def run(self, dataOut, frequency=49.92e6, path=None, fileCadence=1000, dirCadence=36000, metadataCadence=1, **kwargs):
807 811 '''
808 812 This method will be called many times so here you should put all your code
809 813 Inputs:
810 814 dataOut: object with the data
811 815 '''
812 816 # print dataOut.__dict__
813 817 self.dataOut = dataOut
814 818 if not self.isConfig:
815 819 self.setup(dataOut, path, frequency, fileCadence,
816 820 dirCadence, metadataCadence, **kwargs)
817 821 self.writeMetadata()
818 822
819 823 self.writeData()
820 824
821 825 ## self.currentSample += 1
822 826 # if self.dataOut.flagDataAsBlock or self.currentSample == 1:
823 827 # self.writeMetadata()
824 828 ## if self.currentSample == self.__nProfiles: self.currentSample = 0
825 829
826 830 return dataOut# en la version 2.7 no aparece este return
827 831
828 832 def close(self):
829 833 print('[Writing] - Closing files ')
830 834 print('Average of writing to digital rf format is ', self.oldAverage * 1000)
831 835 try:
832 836 self.digitalWriteObj.close()
833 837 except:
834 838 pass
Show file before | Show file after | Hide comments
@@ -1,714 +1,731
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85 self.dataOut = Parameters()
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 if not self.ext.startswith('.'):
102 102 self.ext = '.{}'.format(self.ext)
103 103
104 104 if self.online:
105 105 log.log("Searching files in online mode...", self.name)
106 106
107 107 for nTries in range(self.nTries):
108 108 fullpath = self.searchFilesOnLine(self.path, self.startDate,
109 109 self.endDate, self.expLabel, self.ext, self.walk,
110 110 self.filefmt, self.folderfmt)
111 111 try:
112 112 fullpath = next(fullpath)
113 113 except:
114 114 fullpath = None
115 115
116 116 if fullpath:
117 117 break
118 118
119 119 log.warning(
120 120 'Waiting {} sec for a valid file in {}: try {} ...'.format(
121 121 self.delay, self.path, nTries + 1),
122 122 self.name)
123 123 time.sleep(self.delay)
124 124
125 125 if not(fullpath):
126 126 raise schainpy.admin.SchainError(
127 127 'There isn\'t any valid file in {}'.format(self.path))
128 128
129 129 pathname, filename = os.path.split(fullpath)
130 130 self.year = int(filename[1:5])
131 131 self.doy = int(filename[5:8])
132 132 self.set = int(filename[8:11]) - 1
133 133 else:
134 134 log.log("Searching files in {}".format(self.path), self.name)
135 135 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
136 136 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
137 137
138 138 self.setNextFile()
139 139
140 140 return
141 141
142 142 def readFirstHeader(self):
143 143 '''Read metadata and data'''
144 144
145 145 self.__readMetadata()
146 146 self.__readData()
147 147 self.__setBlockList()
148 148
149 149 if 'type' in self.meta:
150 150 self.dataOut = eval(self.meta['type'])()
151 151
152 152 for attr in self.meta:
153 153 setattr(self.dataOut, attr, self.meta[attr])
154 154
155 155 self.blockIndex = 0
156 156
157 157 return
158 158
159 159 def __setBlockList(self):
160 160 '''
161 161 Selects the data within the times defined
162 162
163 163 self.fp
164 164 self.startTime
165 165 self.endTime
166 166 self.blockList
167 167 self.blocksPerFile
168 168
169 169 '''
170 170
171 171 startTime = self.startTime
172 172 endTime = self.endTime
173 173 thisUtcTime = self.data['utctime'] + self.utcoffset
174 174 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
175 175 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
176 176
177 177 thisDate = thisDatetime.date()
178 178 thisTime = thisDatetime.time()
179 179
180 180 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
181 181 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
182 182
183 183 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
184 184
185 185 self.blockList = ind
186 186 self.blocksPerFile = len(ind)
187 187 return
188 188
189 189 def __readMetadata(self):
190 190 '''
191 191 Reads Metadata
192 192 '''
193 193
194 194 meta = {}
195 195
196 196 if self.description:
197 197 for key, value in self.description['Metadata'].items():
198 198 meta[key] = self.fp[value][()]
199 199 else:
200 200 grp = self.fp['Metadata']
201 201 for name in grp:
202 202 meta[name] = grp[name][()]
203 203
204 204 if self.extras:
205 205 for key, value in self.extras.items():
206 206 meta[key] = value
207 207 self.meta = meta
208 208
209 209 return
210 210
211 211 def __readData(self):
212 212
213 213 data = {}
214 214
215 215 if self.description:
216 216 for key, value in self.description['Data'].items():
217 217 if isinstance(value, str):
218 218 if isinstance(self.fp[value], h5py.Dataset):
219 219 data[key] = self.fp[value][()]
220 220 elif isinstance(self.fp[value], h5py.Group):
221 221 array = []
222 222 for ch in self.fp[value]:
223 223 array.append(self.fp[value][ch][()])
224 224 data[key] = numpy.array(array)
225 225 elif isinstance(value, list):
226 226 array = []
227 227 for ch in value:
228 228 array.append(self.fp[ch][()])
229 229 data[key] = numpy.array(array)
230 230 else:
231 231 grp = self.fp['Data']
232 232 for name in grp:
233 233 if isinstance(grp[name], h5py.Dataset):
234 234 array = grp[name][()]
235 235 elif isinstance(grp[name], h5py.Group):
236 236 array = []
237 237 for ch in grp[name]:
238 238 array.append(grp[name][ch][()])
239 239 array = numpy.array(array)
240 240 else:
241 241 log.warning('Unknown type: {}'.format(name))
242 242
243 243 if name in self.description:
244 244 key = self.description[name]
245 245 else:
246 246 key = name
247 247 data[key] = array
248 248
249 249 self.data = data
250 250 return
251 251
252 252 def getData(self):
253 253
254 254 for attr in self.data:
255 255 if self.data[attr].ndim == 1:
256 256 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
257 257 else:
258 258 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
259 259
260 260 self.dataOut.flagNoData = False
261 261 self.blockIndex += 1
262 262
263 263 log.log("Block No. {}/{} -> {}".format(
264 264 self.blockIndex,
265 265 self.blocksPerFile,
266 266 self.dataOut.datatime.ctime()), self.name)
267 267
268 268 return
269 269
270 270 def run(self, **kwargs):
271 271
272 272 if not(self.isConfig):
273 273 self.setup(**kwargs)
274 274 self.isConfig = True
275 275
276 276 if self.blockIndex == self.blocksPerFile:
277 277 self.setNextFile()
278 278
279 279 self.getData()
280 280
281 281 return
282 282
283 283 @MPDecorator
284 284 class HDFWriter(Operation):
285 285 """Operation to write HDF5 files.
286 286
287 287 The HDF5 file contains by default two groups Data and Metadata where
288 288 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
289 289 parameters, data attributes are normaly time dependent where the metadata
290 290 are not.
291 291 It is possible to customize the structure of the HDF5 file with the
292 292 optional description parameter see the examples.
293 293
294 294 Parameters:
295 295 -----------
296 296 path : str
297 297 Path where files will be saved.
298 298 blocksPerFile : int
299 299 Number of blocks per file
300 300 metadataList : list
301 301 List of the dataOut attributes that will be saved as metadata
302 302 dataList : int
303 303 List of the dataOut attributes that will be saved as data
304 304 setType : bool
305 305 If True the name of the files corresponds to the timestamp of the data
306 306 description : dict, optional
307 307 Dictionary with the desired description of the HDF5 file
308 308
309 309 Examples
310 310 --------
311 311
312 312 desc = {
313 313 'data_output': {'winds': ['z', 'w', 'v']},
314 314 'utctime': 'timestamps',
315 315 'heightList': 'heights'
316 316 }
317 317 desc = {
318 318 'data_output': ['z', 'w', 'v'],
319 319 'utctime': 'timestamps',
320 320 'heightList': 'heights'
321 321 }
322 322 desc = {
323 323 'Data': {
324 324 'data_output': 'winds',
325 325 'utctime': 'timestamps'
326 326 },
327 327 'Metadata': {
328 328 'heightList': 'heights'
329 329 }
330 330 }
331 331
332 332 writer = proc_unit.addOperation(name='HDFWriter')
333 333 writer.addParameter(name='path', value='/path/to/file')
334 334 writer.addParameter(name='blocksPerFile', value='32')
335 335 writer.addParameter(name='metadataList', value='heightList,timeZone')
336 336 writer.addParameter(name='dataList',value='data_output,utctime')
337 337 # writer.addParameter(name='description',value=json.dumps(desc))
338 338
339 339 """
340 340
341 341 ext = ".hdf5"
342 342 optchar = "D"
343 343 filename = None
344 344 path = None
345 345 setFile = None
346 346 fp = None
347 347 firsttime = True
348 348 #Configurations
349 349 blocksPerFile = None
350 350 blockIndex = None
351 351 dataOut = None
352 352 #Data Arrays
353 353 dataList = None
354 354 metadataList = None
355 355 currentDay = None
356 356 lastTime = None
357 357 last_Azipos = None
358 358 last_Elepos = None
359 359 mode = None
360 360 #-----------------------
361 361 Typename = None
362 362
363 363
364 364
365 365 def __init__(self):
366 366
367 367 Operation.__init__(self)
368 368 return
369 369
370 370
371 371 def set_kwargs(self, **kwargs):
372 372
373 373 for key, value in kwargs.items():
374 374 setattr(self, key, value)
375 375
376 376 def set_kwargs_obj(self,obj, **kwargs):
377 377
378 378 for key, value in kwargs.items():
379 379 setattr(obj, key, value)
380 380
381 381 def generalFlag(self):
382 382 ####rint("GENERALFLAG")
383 383 if self.mode== "weather":
384 384 if self.last_Azipos == None:
385 385 tmp = self.dataOut.azimuth
386 386 ####print("ang azimuth writer",tmp)
387 387 self.last_Azipos = tmp
388 388 flag = False
389 389 return flag
390 390 ####print("ang_azimuth writer",self.dataOut.azimuth)
391 391 result = self.dataOut.azimuth - self.last_Azipos
392 392 self.last_Azipos = self.dataOut.azimuth
393 393 if result<0:
394 394 flag = True
395 395 return flag
396 396
397 397 def generalFlag_vRF(self):
398 398 ####rint("GENERALFLAG")
399 399
400 400 try:
401 401 self.dataOut.flagBlock360Done
402 402 return self.dataOut.flagBlock360Done
403 403 except:
404 404 return 0
405 405
406 406
407 407 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None,type_data=None,**kwargs):
408 408 self.path = path
409 409 self.blocksPerFile = blocksPerFile
410 410 self.metadataList = metadataList
411 411 self.dataList = [s.strip() for s in dataList]
412 412 self.setType = setType
413 413 if self.mode == "weather":
414 414 self.setType = "weather"
415 415 self.set_kwargs(**kwargs)
416 416 self.set_kwargs_obj(self.dataOut,**kwargs)
417 417
418 418
419 419 self.description = description
420 420 self.type_data=type_data
421 421
422 422 if self.metadataList is None:
423 423 self.metadataList = self.dataOut.metadata_list
424 424
425 425 tableList = []
426 426 dsList = []
427 427
428 428 for i in range(len(self.dataList)):
429 429 dsDict = {}
430 430 if hasattr(self.dataOut, self.dataList[i]):
431 431 dataAux = getattr(self.dataOut, self.dataList[i])
432 432 dsDict['variable'] = self.dataList[i]
433 433 else:
434 434 log.warning('Attribute {} not found in dataOut', self.name)
435 435 continue
436 436
437 437 if dataAux is None:
438 438 continue
439 439 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
440 440 dsDict['nDim'] = 0
441 441 else:
442 442 dsDict['nDim'] = len(dataAux.shape)
443 443 dsDict['shape'] = dataAux.shape
444 444 dsDict['dsNumber'] = dataAux.shape[0]
445 445 dsDict['dtype'] = dataAux.dtype
446 446 dsList.append(dsDict)
447 447
448 448 self.dsList = dsList
449 449 self.currentDay = self.dataOut.datatime.date()
450 450
451 451 def timeFlag(self):
452 452 currentTime = self.dataOut.utctime
453 453 timeTuple = time.localtime(currentTime)
454 454 dataDay = timeTuple.tm_yday
455 455
456 456 if self.lastTime is None:
457 457 self.lastTime = currentTime
458 458 self.currentDay = dataDay
459 459 return False
460 460
461 461 timeDiff = currentTime - self.lastTime
462 462
463 463 #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
464 464 if dataDay != self.currentDay:
465 465 self.currentDay = dataDay
466 466 return True
467 467 elif timeDiff > 3*60*60:
468 468 self.lastTime = currentTime
469 469 return True
470 470 else:
471 471 self.lastTime = currentTime
472 472 return False
473 473
474 474 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
475 475 dataList=[], setType=None, description={},mode= None,type_data=None,Reset = False,**kwargs):
476 476
477 477 if Reset:
478 478 self.isConfig = False
479 479 self.closeFile()
480 480 self.lastTime = None
481 481 self.blockIndex = 0
482 482
483 483 self.dataOut = dataOut
484 484 self.mode = mode
485 self.var = dataList[0]
486
485 487 if not(self.isConfig):
486 488 self.setup(path=path, blocksPerFile=blocksPerFile,
487 489 metadataList=metadataList, dataList=dataList,
488 490 setType=setType, description=description,type_data=type_data,**kwargs)
489 491
490 492 self.isConfig = True
491 493 self.setNextFile()
492 494
493 495 self.putData()
494 496 return
495 497
496 498 def setNextFile(self):
497 499 ###print("HELLO WORLD--------------------------------")
498 500 ext = self.ext
499 501 path = self.path
500 502 setFile = self.setFile
501 503 type_data = self.type_data
502 504
503 505 timeTuple = time.localtime(self.dataOut.utctime)
504 506 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
505 507 fullpath = os.path.join(path, subfolder)
506 508
507 509 if os.path.exists(fullpath):
508 510 filesList = os.listdir(fullpath)
509 511 filesList = [k for k in filesList if k.startswith(self.optchar)]
510 512 if len( filesList ) > 0:
511 513 filesList = sorted(filesList, key=str.lower)
512 514 filen = filesList[-1]
513 515 # el filename debera tener el siguiente formato
514 516 # 0 1234 567 89A BCDE (hex)
515 517 # x YYYY DDD SSS .ext
516 518 if isNumber(filen[8:11]):
517 519 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
518 520 else:
519 521 setFile = -1
520 522 else:
521 523 setFile = -1 #inicializo mi contador de seteo
522 524 else:
523 525 os.makedirs(fullpath)
524 526 setFile = -1 #inicializo mi contador de seteo
525 527
526 528 ###print("**************************",self.setType)
527 529 if self.setType is None:
528 530 setFile += 1
529 531 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
530 532 timeTuple.tm_year,
531 533 timeTuple.tm_yday,
532 534 setFile,
533 535 ext )
534 536 elif self.setType == "weather":
535 print("HOLA AMIGOS")
536 wr_exp = self.dataOut.wr_exp
537 if wr_exp== "PPI":
538 wr_type = 'E'
539 ang_ = numpy.mean(self.dataOut.elevation)
540 else:
541 wr_type = 'A'
542 ang_ = numpy.mean(self.dataOut.azimuth)
543
544 wr_writer = '%s%s%2.1f%s'%('-',
545 wr_type,
546 ang_,
547 '-')
548 ###print("wr_writer********************",wr_writer)
549 file = '%s%4.4d%2.2d%2.2d%s%2.2d%2.2d%2.2d%s%s%s' % (self.optchar,
550 timeTuple.tm_year,
551 timeTuple.tm_mon,
552 timeTuple.tm_mday,
553 '-',
554 timeTuple.tm_hour,
555 timeTuple.tm_min,
556 timeTuple.tm_sec,
557 wr_writer,
558 type_data,
559 ext )
560 ###print("FILENAME", file)
561 537
538 if self.var.lower() == 'Zdb'.lower():
539 wr_type = 'Z'
540 elif self.var.lower() == 'Zdb_D'.lower():
541 wr_type = 'D'
542 elif self.var.lower() == 'PhiD_P'.lower():
543 wr_type = 'P'
544 elif self.var.lower() == 'RhoHV_R'.lower():
545 wr_type = 'R'
546 elif self.var.lower() == 'velRadial_V'.lower():
547 wr_type = 'V'
548 elif self.var.lower() == 'Sigmav_W'.lower():
549 wr_type = 'S'
550 elif self.var.lower() == 'dataPP_POWER'.lower():
551 wr_type = 'Pow'
552 elif self.var.lower() == 'dataPP_DOP'.lower():
553 wr_type = 'Dop'
554
555
556 #Z_SOPHy_El10.0_20200505_14:02:15.h5
557 #Z_SOPHy_Az40.0_20200505_14:02:15.h5
558 if self.dataOut.flagMode == 1: #'AZI' #PPI
559 ang_type = 'El'
560 ang_ = round(numpy.mean(self.dataOut.data_ele),1)
561 elif self.dataOut.flagMode == 0: #'ELE' #RHI
562 ang_type = 'Az'
563 ang_ = round(numpy.mean(self.dataOut.data_azi),1)
564
565 file = '%s%s%s%2.1f%s%2.2d%2.2d%2.2d%s%2.2d%2.2d%2.2d%s' % (wr_type,
566 '_SOPHy_',
567 ang_type,
568 ang_,
569 '_',
570 timeTuple.tm_year,
571 timeTuple.tm_mon,
572 timeTuple.tm_mday,
573 '_',
574 timeTuple.tm_hour,
575 timeTuple.tm_min,
576 timeTuple.tm_sec,
577 ext )
562 578
563 579 else:
564 580 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
565 581 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
566 582 timeTuple.tm_year,
567 583 timeTuple.tm_yday,
568 584 setFile,
569 585 ext )
570 586
571 587 self.filename = os.path.join( path, subfolder, file )
572 588
573 589 #Setting HDF5 File
574
590 print("filename",self.filename)
575 591 self.fp = h5py.File(self.filename, 'w')
576 592 #write metadata
577 593 self.writeMetadata(self.fp)
578 594 #Write data
579 595 self.writeData(self.fp)
580 596
581 597 def getLabel(self, name, x=None):
582 598
583 599 if x is None:
584 600 if 'Data' in self.description:
585 601 data = self.description['Data']
586 602 if 'Metadata' in self.description:
587 603 data.update(self.description['Metadata'])
588 604 else:
589 605 data = self.description
590 606 if name in data:
591 607 if isinstance(data[name], str):
592 608 return data[name]
593 609 elif isinstance(data[name], list):
594 610 return None
595 611 elif isinstance(data[name], dict):
596 612 for key, value in data[name].items():
597 613 return key
598 614 return name
599 615 else:
600 616 if 'Metadata' in self.description:
601 617 meta = self.description['Metadata']
602 618 else:
603 619 meta = self.description
604 620 if name in meta:
605 621 if isinstance(meta[name], list):
606 622 return meta[name][x]
607 623 elif isinstance(meta[name], dict):
608 624 for key, value in meta[name].items():
609 625 return value[x]
610 626 if 'cspc' in name:
611 627 return 'pair{:02d}'.format(x)
612 628 else:
613 629 return 'channel{:02d}'.format(x)
614 630
615 631 def writeMetadata(self, fp):
616 632
617 633 if self.description:
618 634 if 'Metadata' in self.description:
619 635 grp = fp.create_group('Metadata')
620 636 else:
621 637 grp = fp
622 638 else:
623 639 grp = fp.create_group('Metadata')
624 640
625 641 for i in range(len(self.metadataList)):
626 642 if not hasattr(self.dataOut, self.metadataList[i]):
627 643 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
628 644 continue
629 645 value = getattr(self.dataOut, self.metadataList[i])
630 646 if isinstance(value, bool):
631 647 if value is True:
632 648 value = 1
633 649 else:
634 650 value = 0
635 651 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
636 652 return
637 653
638 654 def writeData(self, fp):
639 655
640 656 if self.description:
641 657 if 'Data' in self.description:
642 658 grp = fp.create_group('Data')
643 659 else:
644 660 grp = fp
645 661 else:
646 662 grp = fp.create_group('Data')
647 663
648 664 dtsets = []
649 665 data = []
650 666
651 667 for dsInfo in self.dsList:
652 668 if dsInfo['nDim'] == 0:
653 669 ds = grp.create_dataset(
654 670 self.getLabel(dsInfo['variable']),
655 671 (self.blocksPerFile, ),
656 672 chunks=True,
657 673 dtype=numpy.float64)
658 674 dtsets.append(ds)
659 675 data.append((dsInfo['variable'], -1))
660 676 else:
661 677 label = self.getLabel(dsInfo['variable'])
662 678 if label is not None:
663 679 sgrp = grp.create_group(label)
664 680 else:
665 681 sgrp = grp
666 682 for i in range(dsInfo['dsNumber']):
667 683 ds = sgrp.create_dataset(
668 684 self.getLabel(dsInfo['variable'], i),
669 685 (self.blocksPerFile, ) + dsInfo['shape'][1:],
670 686 chunks=True,
671 687 dtype=dsInfo['dtype'])
672 688 dtsets.append(ds)
673 689 data.append((dsInfo['variable'], i))
674 690 fp.flush()
675 691
676 692 log.log('Creating file: {}'.format(fp.filename), self.name)
677 693
678 694 self.ds = dtsets
679 695 self.data = data
680 696 self.firsttime = True
681 697 self.blockIndex = 0
682 698 return
683 699
684 700 def putData(self):
701
685 702 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():# or self.generalFlag_vRF():
686 703 self.closeFile()
687 704 self.setNextFile()
688 705
689 706 for i, ds in enumerate(self.ds):
690 707 attr, ch = self.data[i]
691 708 if ch == -1:
692 709 ds[self.blockIndex] = getattr(self.dataOut, attr)
693 710 else:
694 711 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
695 712
696 713 self.fp.flush()
697 714 self.blockIndex += 1
698 715 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
699 716
700 717 return
701 718
702 719 def closeFile(self):
703 720
704 721 if self.blockIndex != self.blocksPerFile:
705 722 for ds in self.ds:
706 723 ds.resize(self.blockIndex, axis=0)
707 724
708 725 if self.fp:
709 726 self.fp.flush()
710 727 self.fp.close()
711 728
712 729 def close(self):
713 730
714 731 self.closeFile()
Show file before | Show file after | Hide comments
@@ -1,226 +1,236
1 1 '''
2 2 Base clases to create Processing units and operations, the MPDecorator
3 3 must be used in plotting and writing operations to allow to run as an
4 4 external process.
5 5 '''
6 6
7 7 import os
8 8 import inspect
9 9 import zmq
10 10 import time
11 11 import pickle
12 12 import traceback
13 13 from threading import Thread
14 14 from multiprocessing import Process, Queue
15 15 from schainpy.utils import log
16 16
17 17 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
18 18
19 19 class ProcessingUnit(object):
20 20 '''
21 21 Base class to create Signal Chain Units
22 22 '''
23 23
24 24 proc_type = 'processing'
25 25
26 26 def __init__(self):
27 27
28 28 self.dataIn = None
29 29 self.dataOut = None
30 30 self.isConfig = False
31 31 self.operations = []
32 32 self.name = 'Test'
33 33 self.inputs = []
34 34
35 35 def setInput(self, unit):
36 36
37 37 attr = 'dataIn'
38 38 for i, u in enumerate(unit):
39 39 if i==0:
40 40 self.dataIn = u.dataOut
41 41 self.inputs.append('dataIn')
42 42 else:
43 43 setattr(self, 'dataIn{}'.format(i), u.dataOut)
44 44 self.inputs.append('dataIn{}'.format(i))
45 45
46 46 def getAllowedArgs(self):
47 47 if hasattr(self, '__attrs__'):
48 48 return self.__attrs__
49 49 else:
50 50 return inspect.getargspec(self.run).args
51 51
52 52 def addOperation(self, conf, operation):
53 53 '''
54 54 '''
55 55
56 56 self.operations.append((operation, conf.type, conf.getKwargs()))
57 57
58 58 def getOperationObj(self, objId):
59 59
60 60 if objId not in list(self.operations.keys()):
61 61 return None
62 62
63 63 return self.operations[objId]
64 64
65 65 def call(self, **kwargs):
66 66 '''
67 67 '''
68 68
69 69 try:
70 70 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
71 71 return self.dataIn.isReady()
72 72 elif self.dataIn is None or not self.dataIn.error:
73 73 self.run(**kwargs)
74 74 elif self.dataIn.error:
75 75 self.dataOut.error = self.dataIn.error
76 76 self.dataOut.flagNoData = True
77 77 except:
78 78 err = traceback.format_exc()
79 79 if 'SchainWarning' in err:
80 80 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
81 81 elif 'SchainError' in err:
82 82 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
83 83 else:
84 84 log.error(err, self.name)
85 85 self.dataOut.error = True
86 86 ##### correcion de la declaracion Out
87 87 for op, optype, opkwargs in self.operations:
88 88 aux = self.dataOut.copy()
89 if optype == 'other' and not self.dataOut.flagNoData:
89 '''
90 print("op",op)
91 try:
92 print("runNextOp",self.dataOut.runNextOp)
93 except:
94 pass
95 '''
96 if not hasattr(self.dataOut, 'runNextOp'):
97 self.dataOut.runNextOp = False
98 if optype == 'other' and (not self.dataOut.flagNoData or self.dataOut.runNextOp):
99 #if optype == 'other' and not self.dataOut.flagNoData:
90 100 self.dataOut = op.run(self.dataOut, **opkwargs)
91 101 elif optype == 'external' and not self.dataOut.flagNoData:
92 102 #op.queue.put(self.dataOut)
93 103 op.queue.put(aux)
94 104 elif optype == 'external' and self.dataOut.error:
95 105 #op.queue.put(self.dataOut)
96 106 op.queue.put(aux)
97 107
98 108 try:
99 109 if self.dataOut.runNextUnit:
100 110 runNextUnit = self.dataOut.runNextUnit
101 111
102 112 else:
103 113 runNextUnit = self.dataOut.isReady()
104 114 except:
105 115 runNextUnit = self.dataOut.isReady()
106 116
107 117 return 'Error' if self.dataOut.error else runNextUnit
108 118
109 119 def setup(self):
110 120
111 121 raise NotImplementedError
112 122
113 123 def run(self):
114 124
115 125 raise NotImplementedError
116 126
117 127 def close(self):
118 128
119 129 return
120 130
121 131
122 132 class Operation(object):
123 133
124 134 '''
125 135 '''
126 136
127 137 proc_type = 'operation'
128 138
129 139 def __init__(self):
130 140
131 141 self.id = None
132 142 self.isConfig = False
133 143
134 144 if not hasattr(self, 'name'):
135 145 self.name = self.__class__.__name__
136 146
137 147 def getAllowedArgs(self):
138 148 if hasattr(self, '__attrs__'):
139 149 return self.__attrs__
140 150 else:
141 151 return inspect.getargspec(self.run).args
142 152
143 153 def setup(self):
144 154
145 155 self.isConfig = True
146 156
147 157 raise NotImplementedError
148 158
149 159 def run(self, dataIn, **kwargs):
150 160 """
151 161 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
152 162 atributos del objeto dataIn.
153 163
154 164 Input:
155 165
156 166 dataIn : objeto del tipo JROData
157 167
158 168 Return:
159 169
160 170 None
161 171
162 172 Affected:
163 173 __buffer : buffer de recepcion de datos.
164 174
165 175 """
166 176 if not self.isConfig:
167 177 self.setup(**kwargs)
168 178
169 179 raise NotImplementedError
170 180
171 181 def close(self):
172 182
173 183 return
174 184
175 185
176 186 def MPDecorator(BaseClass):
177 187 """
178 188 Multiprocessing class decorator
179 189
180 190 This function add multiprocessing features to a BaseClass.
181 191 """
182 192
183 193 class MPClass(BaseClass, Process):
184 194
185 195 def __init__(self, *args, **kwargs):
186 196 super(MPClass, self).__init__()
187 197 Process.__init__(self)
188 198
189 199 self.args = args
190 200 self.kwargs = kwargs
191 201 self.t = time.time()
192 202 self.op_type = 'external'
193 203 self.name = BaseClass.__name__
194 204 self.__doc__ = BaseClass.__doc__
195 205
196 206 if 'plot' in self.name.lower() and not self.name.endswith('_'):
197 207 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
198 208
199 209 self.start_time = time.time()
200 210 self.err_queue = args[3]
201 211 self.queue = Queue(maxsize=QUEUE_SIZE)
202 212 self.myrun = BaseClass.run
203 213
204 214 def run(self):
205 215
206 216 while True:
207 217
208 218 dataOut = self.queue.get()
209 219
210 220 if not dataOut.error:
211 221 try:
212 222 BaseClass.run(self, dataOut, **self.kwargs)
213 223 except:
214 224 err = traceback.format_exc()
215 225 log.error(err, self.name)
216 226 else:
217 227 break
218 228
219 229 self.close()
220 230
221 231 def close(self):
222 232
223 233 BaseClass.close(self)
224 234 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
225 235
226 236 return MPClass
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,1860 +1,1861
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei and maxHei:
168 168
169 169 if (minHei < self.dataOut.heightList[0]):
170 170 minHei = self.dataOut.heightList[0]
171 171
172 172 if (maxHei > self.dataOut.heightList[-1]):
173 173 maxHei = self.dataOut.heightList[-1]
174 174
175 175 minIndex = 0
176 176 maxIndex = 0
177 177 heights = self.dataOut.heightList
178 178
179 179 inda = numpy.where(heights >= minHei)
180 180 indb = numpy.where(heights <= maxHei)
181 181
182 182 try:
183 183 minIndex = inda[0][0]
184 184 except:
185 185 minIndex = 0
186 186
187 187 try:
188 188 maxIndex = indb[0][-1]
189 189 except:
190 190 maxIndex = len(heights)
191 191
192 192 self.selectHeightsByIndex(minIndex, maxIndex)
193 193
194 194 return self.dataOut
195 195
196 196 def selectHeightsByIndex(self, minIndex, maxIndex):
197 197 """
198 198 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
199 199 minIndex <= index <= maxIndex
200 200
201 201 Input:
202 202 minIndex : valor de indice minimo de altura a considerar
203 203 maxIndex : valor de indice maximo de altura a considerar
204 204
205 205 Affected:
206 206 self.dataOut.data
207 207 self.dataOut.heightList
208 208
209 209 Return:
210 210 1 si el metodo se ejecuto con exito caso contrario devuelve 0
211 211 """
212 212
213 213 if self.dataOut.type == 'Voltage':
214 214 if (minIndex < 0) or (minIndex > maxIndex):
215 215 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
216 216
217 217 if (maxIndex >= self.dataOut.nHeights):
218 218 maxIndex = self.dataOut.nHeights
219 219 #print("shapeeee",self.dataOut.data.shape)
220 220 #voltage
221 221 if self.dataOut.flagDataAsBlock:
222 222 """
223 223 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
224 224 """
225 225 data = self.dataOut.data[:,:, minIndex:maxIndex]
226 226 else:
227 227 data = self.dataOut.data[:, minIndex:maxIndex]
228 228
229 229 # firstHeight = self.dataOut.heightList[minIndex]
230 230
231 231 self.dataOut.data = data
232 232 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
233 233
234 234 if self.dataOut.nHeights <= 1:
235 235 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
236 236 elif self.dataOut.type == 'Spectra':
237 237 if (minIndex < 0) or (minIndex > maxIndex):
238 238 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
239 239 minIndex, maxIndex))
240 240
241 241 if (maxIndex >= self.dataOut.nHeights):
242 242 maxIndex = self.dataOut.nHeights - 1
243 243
244 244 # Spectra
245 245 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
246 246
247 247 data_cspc = None
248 248 if self.dataOut.data_cspc is not None:
249 249 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_dc = None
252 252 if self.dataOut.data_dc is not None:
253 253 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
254 254
255 255 self.dataOut.data_spc = data_spc
256 256 self.dataOut.data_cspc = data_cspc
257 257 self.dataOut.data_dc = data_dc
258 258
259 259 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
260 260
261 261 return 1
262 262
263 263
264 264 class filterByHeights(Operation):
265 265
266 266 def run(self, dataOut, window):
267 267
268 268 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
269 269
270 270 if window == None:
271 271 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
272 272
273 273 newdelta = deltaHeight * window
274 274 r = dataOut.nHeights % window
275 275 newheights = (dataOut.nHeights-r)/window
276 276
277 277 if newheights <= 1:
278 278 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
279 279
280 280 if dataOut.flagDataAsBlock:
281 281 """
282 282 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
283 283 """
284 284 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
285 285 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
286 286 buffer = numpy.sum(buffer,3)
287 287
288 288 else:
289 289 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
290 290 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
291 291 buffer = numpy.sum(buffer,2)
292 292
293 293 dataOut.data = buffer
294 294 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
295 295 dataOut.windowOfFilter = window
296 296
297 297 return dataOut
298 298
299 299
300 300 class setH0(Operation):
301 301
302 302 def run(self, dataOut, h0, deltaHeight = None):
303 303
304 304 if not deltaHeight:
305 305 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
306 306
307 307 nHeights = dataOut.nHeights
308 308
309 309 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
310 310
311 311 dataOut.heightList = newHeiRange
312 312
313 313 return dataOut
314 314
315 315
316 316 class deFlip(Operation):
317 317
318 318 def run(self, dataOut, channelList = []):
319 319
320 320 data = dataOut.data.copy()
321 321
322 322 if dataOut.flagDataAsBlock:
323 323 flip = self.flip
324 324 profileList = list(range(dataOut.nProfiles))
325 325
326 326 if not channelList:
327 327 for thisProfile in profileList:
328 328 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
329 329 flip *= -1.0
330 330 else:
331 331 for thisChannel in channelList:
332 332 if thisChannel not in dataOut.channelList:
333 333 continue
334 334
335 335 for thisProfile in profileList:
336 336 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
337 337 flip *= -1.0
338 338
339 339 self.flip = flip
340 340
341 341 else:
342 342 if not channelList:
343 343 data[:,:] = data[:,:]*self.flip
344 344 else:
345 345 for thisChannel in channelList:
346 346 if thisChannel not in dataOut.channelList:
347 347 continue
348 348
349 349 data[thisChannel,:] = data[thisChannel,:]*self.flip
350 350
351 351 self.flip *= -1.
352 352
353 353 dataOut.data = data
354 354
355 355 return dataOut
356 356
357 357
358 358 class setAttribute(Operation):
359 359 '''
360 360 Set an arbitrary attribute(s) to dataOut
361 361 '''
362 362
363 363 def __init__(self):
364 364
365 365 Operation.__init__(self)
366 366 self._ready = False
367 367
368 368 def run(self, dataOut, **kwargs):
369 369
370 370 for key, value in kwargs.items():
371 371 setattr(dataOut, key, value)
372 372
373 373 return dataOut
374 374
375 375
376 376 @MPDecorator
377 377 class printAttribute(Operation):
378 378 '''
379 379 Print an arbitrary attribute of dataOut
380 380 '''
381 381
382 382 def __init__(self):
383 383
384 384 Operation.__init__(self)
385 385
386 386 def run(self, dataOut, attributes):
387 387
388 388 if isinstance(attributes, str):
389 389 attributes = [attributes]
390 390 for attr in attributes:
391 391 if hasattr(dataOut, attr):
392 392 log.log(getattr(dataOut, attr), attr)
393 393
394 394
395 395 class interpolateHeights(Operation):
396 396
397 397 def run(self, dataOut, topLim, botLim):
398 398 #69 al 72 para julia
399 399 #82-84 para meteoros
400 400 if len(numpy.shape(dataOut.data))==2:
401 401 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
402 402 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
403 403 #dataOut.data[:,botLim:limSup+1] = sampInterp
404 404 dataOut.data[:,botLim:topLim+1] = sampInterp
405 405 else:
406 406 nHeights = dataOut.data.shape[2]
407 407 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
408 408 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
409 409 f = interpolate.interp1d(x, y, axis = 2)
410 410 xnew = numpy.arange(botLim,topLim+1)
411 411 ynew = f(xnew)
412 412 dataOut.data[:,:,botLim:topLim+1] = ynew
413 413
414 414 return dataOut
415 415
416 416
417 417 class CohInt(Operation):
418 418
419 419 isConfig = False
420 420 __profIndex = 0
421 421 __byTime = False
422 422 __initime = None
423 423 __lastdatatime = None
424 424 __integrationtime = None
425 425 __buffer = None
426 426 __bufferStride = []
427 427 __dataReady = False
428 428 __profIndexStride = 0
429 429 __dataToPutStride = False
430 430 n = None
431 431
432 432 def __init__(self, **kwargs):
433 433
434 434 Operation.__init__(self, **kwargs)
435 435
436 436 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
437 437 """
438 438 Set the parameters of the integration class.
439 439
440 440 Inputs:
441 441
442 442 n : Number of coherent integrations
443 443 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
444 444 overlapping :
445 445 """
446 446
447 447 self.__initime = None
448 448 self.__lastdatatime = 0
449 449 self.__buffer = None
450 450 self.__dataReady = False
451 451 self.byblock = byblock
452 452 self.stride = stride
453 453
454 454 if n == None and timeInterval == None:
455 455 raise ValueError("n or timeInterval should be specified ...")
456 456
457 457 if n != None:
458 458 self.n = n
459 459 self.__byTime = False
460 460 else:
461 461 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
462 462 self.n = 9999
463 463 self.__byTime = True
464 464
465 465 if overlapping:
466 466 self.__withOverlapping = True
467 467 self.__buffer = None
468 468 else:
469 469 self.__withOverlapping = False
470 470 self.__buffer = 0
471 471
472 472 self.__profIndex = 0
473 473
474 474 def putData(self, data):
475 475
476 476 """
477 477 Add a profile to the __buffer and increase in one the __profileIndex
478 478
479 479 """
480 480
481 481 if not self.__withOverlapping:
482 482 self.__buffer += data.copy()
483 483 self.__profIndex += 1
484 484 return
485 485
486 486 #Overlapping data
487 487 nChannels, nHeis = data.shape
488 488 data = numpy.reshape(data, (1, nChannels, nHeis))
489 489
490 490 #If the buffer is empty then it takes the data value
491 491 if self.__buffer is None:
492 492 self.__buffer = data
493 493 self.__profIndex += 1
494 494 return
495 495
496 496 #If the buffer length is lower than n then stakcing the data value
497 497 if self.__profIndex < self.n:
498 498 self.__buffer = numpy.vstack((self.__buffer, data))
499 499 self.__profIndex += 1
500 500 return
501 501
502 502 #If the buffer length is equal to n then replacing the last buffer value with the data value
503 503 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
504 504 self.__buffer[self.n-1] = data
505 505 self.__profIndex = self.n
506 506 return
507 507
508 508
509 509 def pushData(self):
510 510 """
511 511 Return the sum of the last profiles and the profiles used in the sum.
512 512
513 513 Affected:
514 514
515 515 self.__profileIndex
516 516
517 517 """
518 518
519 519 if not self.__withOverlapping:
520 520 data = self.__buffer
521 521 n = self.__profIndex
522 522
523 523 self.__buffer = 0
524 524 self.__profIndex = 0
525 525
526 526 return data, n
527 527
528 528 #Integration with Overlapping
529 529 data = numpy.sum(self.__buffer, axis=0)
530 530 # print data
531 531 # raise
532 532 n = self.__profIndex
533 533
534 534 return data, n
535 535
536 536 def byProfiles(self, data):
537 537
538 538 self.__dataReady = False
539 539 avgdata = None
540 540 # n = None
541 541 # print data
542 542 # raise
543 543 self.putData(data)
544 544
545 545 if self.__profIndex == self.n:
546 546 avgdata, n = self.pushData()
547 547 self.__dataReady = True
548 548
549 549 return avgdata
550 550
551 551 def byTime(self, data, datatime):
552 552
553 553 self.__dataReady = False
554 554 avgdata = None
555 555 n = None
556 556
557 557 self.putData(data)
558 558
559 559 if (datatime - self.__initime) >= self.__integrationtime:
560 560 avgdata, n = self.pushData()
561 561 self.n = n
562 562 self.__dataReady = True
563 563
564 564 return avgdata
565 565
566 566 def integrateByStride(self, data, datatime):
567 567 # print data
568 568 if self.__profIndex == 0:
569 569 self.__buffer = [[data.copy(), datatime]]
570 570 else:
571 571 self.__buffer.append([data.copy(),datatime])
572 572 self.__profIndex += 1
573 573 self.__dataReady = False
574 574
575 575 if self.__profIndex == self.n * self.stride :
576 576 self.__dataToPutStride = True
577 577 self.__profIndexStride = 0
578 578 self.__profIndex = 0
579 579 self.__bufferStride = []
580 580 for i in range(self.stride):
581 581 current = self.__buffer[i::self.stride]
582 582 data = numpy.sum([t[0] for t in current], axis=0)
583 583 avgdatatime = numpy.average([t[1] for t in current])
584 584 # print data
585 585 self.__bufferStride.append((data, avgdatatime))
586 586
587 587 if self.__dataToPutStride:
588 588 self.__dataReady = True
589 589 self.__profIndexStride += 1
590 590 if self.__profIndexStride == self.stride:
591 591 self.__dataToPutStride = False
592 592 # print self.__bufferStride[self.__profIndexStride - 1]
593 593 # raise
594 594 return self.__bufferStride[self.__profIndexStride - 1]
595 595
596 596
597 597 return None, None
598 598
599 599 def integrate(self, data, datatime=None):
600 600
601 601 if self.__initime == None:
602 602 self.__initime = datatime
603 603
604 604 if self.__byTime:
605 605 avgdata = self.byTime(data, datatime)
606 606 else:
607 607 avgdata = self.byProfiles(data)
608 608
609 609
610 610 self.__lastdatatime = datatime
611 611
612 612 if avgdata is None:
613 613 return None, None
614 614
615 615 avgdatatime = self.__initime
616 616
617 617 deltatime = datatime - self.__lastdatatime
618 618
619 619 if not self.__withOverlapping:
620 620 self.__initime = datatime
621 621 else:
622 622 self.__initime += deltatime
623 623
624 624 return avgdata, avgdatatime
625 625
626 626 def integrateByBlock(self, dataOut):
627 627
628 628 times = int(dataOut.data.shape[1]/self.n)
629 629 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
630 630
631 631 id_min = 0
632 632 id_max = self.n
633 633
634 634 for i in range(times):
635 635 junk = dataOut.data[:,id_min:id_max,:]
636 636 avgdata[:,i,:] = junk.sum(axis=1)
637 637 id_min += self.n
638 638 id_max += self.n
639 639
640 640 timeInterval = dataOut.ippSeconds*self.n
641 641 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
642 642 self.__dataReady = True
643 643 return avgdata, avgdatatime
644 644
645 645 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
646 646
647 647 if not self.isConfig:
648 648 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
649 649 self.isConfig = True
650 650
651 651 if dataOut.flagDataAsBlock:
652 652 """
653 653 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
654 654 """
655 655 avgdata, avgdatatime = self.integrateByBlock(dataOut)
656 656 dataOut.nProfiles /= self.n
657 657 else:
658 658 if stride is None:
659 659 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
660 660 else:
661 661 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
662 662
663 663
664 664 # dataOut.timeInterval *= n
665 665 dataOut.flagNoData = True
666 666
667 667 if self.__dataReady:
668 668 dataOut.data = avgdata
669 669 if not dataOut.flagCohInt:
670 670 dataOut.nCohInt *= self.n
671 671 dataOut.flagCohInt = True
672 dataOut.utctime = avgdatatime
672 ####################################dataOut.utctime = avgdatatime
673 673 # print avgdata, avgdatatime
674 674 # raise
675 675 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
676 676 dataOut.flagNoData = False
677 677 return dataOut
678 678
679 679 class Decoder(Operation):
680 680
681 681 isConfig = False
682 682 __profIndex = 0
683 683
684 684 code = None
685 685
686 686 nCode = None
687 687 nBaud = None
688 688
689 689 def __init__(self, **kwargs):
690 690
691 691 Operation.__init__(self, **kwargs)
692 692
693 693 self.times = None
694 694 self.osamp = None
695 695 # self.__setValues = False
696 696 self.isConfig = False
697 697 self.setupReq = False
698 698 def setup(self, code, osamp, dataOut):
699 699
700 700 self.__profIndex = 0
701 701
702 702 self.code = code
703 703
704 704 self.nCode = len(code)
705 705 self.nBaud = len(code[0])
706 706
707 707 if (osamp != None) and (osamp >1):
708 708 self.osamp = osamp
709 709 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
710 710 self.nBaud = self.nBaud*self.osamp
711 711
712 712 self.__nChannels = dataOut.nChannels
713 713 self.__nProfiles = dataOut.nProfiles
714 714 self.__nHeis = dataOut.nHeights
715 715
716 716 if self.__nHeis < self.nBaud:
717 717 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
718 718
719 719 #Frequency
720 720 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
721 721
722 722 __codeBuffer[:,0:self.nBaud] = self.code
723 723
724 724 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
725 725
726 726 if dataOut.flagDataAsBlock:
727 727
728 728 self.ndatadec = self.__nHeis #- self.nBaud + 1
729 729
730 730 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
731 731
732 732 else:
733 733
734 734 #Time
735 735 self.ndatadec = self.__nHeis #- self.nBaud + 1
736 736
737 737 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
738 738
739 739 def __convolutionInFreq(self, data):
740 740
741 741 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
742 742
743 743 fft_data = numpy.fft.fft(data, axis=1)
744 744
745 745 conv = fft_data*fft_code
746 746
747 747 data = numpy.fft.ifft(conv,axis=1)
748 748
749 749 return data
750 750
751 751 def __convolutionInFreqOpt(self, data):
752 752
753 753 raise NotImplementedError
754 754
755 755 def __convolutionInTime(self, data):
756 756
757 757 code = self.code[self.__profIndex]
758 758 for i in range(self.__nChannels):
759 759 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
760 760
761 761 return self.datadecTime
762 762
763 763 def __convolutionByBlockInTime(self, data):
764 764
765 765 repetitions = int(self.__nProfiles / self.nCode)
766 766 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
767 767 junk = junk.flatten()
768 768 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
769 769 profilesList = range(self.__nProfiles)
770 770
771 771 for i in range(self.__nChannels):
772 772 for j in profilesList:
773 773 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
774 774 return self.datadecTime
775 775
776 776 def __convolutionByBlockInFreq(self, data):
777 777
778 778 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
779 779
780 780
781 781 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
782 782
783 783 fft_data = numpy.fft.fft(data, axis=2)
784 784
785 785 conv = fft_data*fft_code
786 786
787 787 data = numpy.fft.ifft(conv,axis=2)
788 788
789 789 return data
790 790
791 791
792 792 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
793 793
794 794 if dataOut.flagDecodeData:
795 795 print("This data is already decoded, recoding again ...")
796 796
797 797 if not self.isConfig:
798 798
799 799 if code is None:
800 800 if dataOut.code is None:
801 801 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
802 802
803 803 code = dataOut.code
804 804 else:
805 805 code = numpy.array(code).reshape(nCode,nBaud)
806 806 self.setup(code, osamp, dataOut)
807 807
808 808 self.isConfig = True
809 809
810 810 if mode == 3:
811 811 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
812 812
813 813 if times != None:
814 814 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
815 815
816 816 if self.code is None:
817 817 print("Fail decoding: Code is not defined.")
818 818 return
819 819
820 820 self.__nProfiles = dataOut.nProfiles
821 821 datadec = None
822 822
823 823 if mode == 3:
824 824 mode = 0
825 825
826 826 if dataOut.flagDataAsBlock:
827 827 """
828 828 Decoding when data have been read as block,
829 829 """
830 830
831 831 if mode == 0:
832 832 datadec = self.__convolutionByBlockInTime(dataOut.data)
833 833 if mode == 1:
834 834 datadec = self.__convolutionByBlockInFreq(dataOut.data)
835 835 else:
836 836 """
837 837 Decoding when data have been read profile by profile
838 838 """
839 839 if mode == 0:
840 840 datadec = self.__convolutionInTime(dataOut.data)
841 841
842 842 if mode == 1:
843 843 datadec = self.__convolutionInFreq(dataOut.data)
844 844
845 845 if mode == 2:
846 846 datadec = self.__convolutionInFreqOpt(dataOut.data)
847 847
848 848 if datadec is None:
849 849 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
850 850
851 851 dataOut.code = self.code
852 852 dataOut.nCode = self.nCode
853 853 dataOut.nBaud = self.nBaud
854 854
855 855 dataOut.data = datadec
856 856
857 857 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
858 858
859 859 dataOut.flagDecodeData = True #asumo q la data esta decodificada
860 860
861 861 if self.__profIndex == self.nCode-1:
862 862 self.__profIndex = 0
863 863 return dataOut
864 864
865 865 self.__profIndex += 1
866 866
867 867 return dataOut
868 868 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
869 869
870 870
871 871 class ProfileConcat(Operation):
872 872
873 873 isConfig = False
874 874 buffer = None
875 875
876 876 def __init__(self, **kwargs):
877 877
878 878 Operation.__init__(self, **kwargs)
879 879 self.profileIndex = 0
880 880
881 881 def reset(self):
882 882 self.buffer = numpy.zeros_like(self.buffer)
883 883 self.start_index = 0
884 884 self.times = 1
885 885
886 886 def setup(self, data, m, n=1):
887 887 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
888 888 self.nHeights = data.shape[1]#.nHeights
889 889 self.start_index = 0
890 890 self.times = 1
891 891
892 892 def concat(self, data):
893 893
894 894 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
895 895 self.start_index = self.start_index + self.nHeights
896 896
897 897 def run(self, dataOut, m):
898 898 dataOut.flagNoData = True
899 899
900 900 if not self.isConfig:
901 901 self.setup(dataOut.data, m, 1)
902 902 self.isConfig = True
903 903
904 904 if dataOut.flagDataAsBlock:
905 905 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
906 906
907 907 else:
908 908 self.concat(dataOut.data)
909 909 self.times += 1
910 910 if self.times > m:
911 911 dataOut.data = self.buffer
912 912 self.reset()
913 913 dataOut.flagNoData = False
914 914 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
915 915 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
916 916 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
917 917 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
918 918 dataOut.ippSeconds *= m
919 919 return dataOut
920 920
921 921 class ProfileSelector(Operation):
922 922
923 923 profileIndex = None
924 924 # Tamanho total de los perfiles
925 925 nProfiles = None
926 926
927 927 def __init__(self, **kwargs):
928 928
929 929 Operation.__init__(self, **kwargs)
930 930 self.profileIndex = 0
931 931
932 932 def incProfileIndex(self):
933 933
934 934 self.profileIndex += 1
935 935
936 936 if self.profileIndex >= self.nProfiles:
937 937 self.profileIndex = 0
938 938
939 939 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
940 940
941 941 if profileIndex < minIndex:
942 942 return False
943 943
944 944 if profileIndex > maxIndex:
945 945 return False
946 946
947 947 return True
948 948
949 949 def isThisProfileInList(self, profileIndex, profileList):
950 950
951 951 if profileIndex not in profileList:
952 952 return False
953 953
954 954 return True
955 955
956 956 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
957 957 #print("before",dataOut.data.shape)
958 958 """
959 959 ProfileSelector:
960 960
961 961 Inputs:
962 962 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
963 963
964 964 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
965 965
966 966 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
967 967
968 968 """
969 969
970 970 if rangeList is not None:
971 971 if type(rangeList[0]) not in (tuple, list):
972 972 rangeList = [rangeList]
973 973
974 974 dataOut.flagNoData = True
975 975
976 976 if dataOut.flagDataAsBlock:
977 977 """
978 978 data dimension = [nChannels, nProfiles, nHeis]
979 979 """
980 980 if profileList != None:
981 981 dataOut.data = dataOut.data[:,profileList,:]
982 982
983 983 if profileRangeList != None:
984 984 minIndex = profileRangeList[0]
985 985 maxIndex = profileRangeList[1]
986 986 profileList = list(range(minIndex, maxIndex+1))
987 987
988 988 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
989 989
990 990 if rangeList != None:
991 991
992 992 profileList = []
993 993
994 994 for thisRange in rangeList:
995 995 minIndex = thisRange[0]
996 996 maxIndex = thisRange[1]
997 997
998 998 profileList.extend(list(range(minIndex, maxIndex+1)))
999 999
1000 1000 dataOut.data = dataOut.data[:,profileList,:]
1001 1001
1002 1002 dataOut.nProfiles = len(profileList)
1003 1003 dataOut.profileIndex = dataOut.nProfiles - 1
1004 1004 dataOut.flagNoData = False
1005 1005 #print(dataOut.data.shape)
1006 1006 return dataOut
1007 1007
1008 1008 """
1009 1009 data dimension = [nChannels, nHeis]
1010 1010 """
1011 1011
1012 1012 if profileList != None:
1013 1013
1014 1014 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1015 1015
1016 1016 self.nProfiles = len(profileList)
1017 1017 dataOut.nProfiles = self.nProfiles
1018 1018 dataOut.profileIndex = self.profileIndex
1019 1019 dataOut.flagNoData = False
1020 1020
1021 1021 self.incProfileIndex()
1022 1022 return dataOut
1023 1023
1024 1024 if profileRangeList != None:
1025 1025
1026 1026 minIndex = profileRangeList[0]
1027 1027 maxIndex = profileRangeList[1]
1028 1028
1029 1029 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1030 1030
1031 1031 self.nProfiles = maxIndex - minIndex + 1
1032 1032 dataOut.nProfiles = self.nProfiles
1033 1033 dataOut.profileIndex = self.profileIndex
1034 1034 dataOut.flagNoData = False
1035 1035
1036 1036 self.incProfileIndex()
1037 1037 return dataOut
1038 1038
1039 1039 if rangeList != None:
1040 1040
1041 1041 nProfiles = 0
1042 1042
1043 1043 for thisRange in rangeList:
1044 1044 minIndex = thisRange[0]
1045 1045 maxIndex = thisRange[1]
1046 1046
1047 1047 nProfiles += maxIndex - minIndex + 1
1048 1048
1049 1049 for thisRange in rangeList:
1050 1050
1051 1051 minIndex = thisRange[0]
1052 1052 maxIndex = thisRange[1]
1053 1053
1054 1054 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1055 1055
1056 1056 self.nProfiles = nProfiles
1057 1057 dataOut.nProfiles = self.nProfiles
1058 1058 dataOut.profileIndex = self.profileIndex
1059 1059 dataOut.flagNoData = False
1060 1060
1061 1061 self.incProfileIndex()
1062 1062
1063 1063 break
1064 1064
1065 1065 return dataOut
1066 1066
1067 1067
1068 1068 if beam != None: #beam is only for AMISR data
1069 1069 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1070 1070 dataOut.flagNoData = False
1071 1071 dataOut.profileIndex = self.profileIndex
1072 1072
1073 1073 self.incProfileIndex()
1074 1074
1075 1075 return dataOut
1076 1076
1077 1077 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1078 1078
1079 1079
1080 1080 class Reshaper(Operation):
1081 1081
1082 1082 def __init__(self, **kwargs):
1083 1083
1084 1084 Operation.__init__(self, **kwargs)
1085 1085
1086 1086 self.__buffer = None
1087 1087 self.__nitems = 0
1088 1088
1089 1089 def __appendProfile(self, dataOut, nTxs):
1090 1090
1091 1091 if self.__buffer is None:
1092 1092 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1093 1093 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1094 1094
1095 1095 ini = dataOut.nHeights * self.__nitems
1096 1096 end = ini + dataOut.nHeights
1097 1097
1098 1098 self.__buffer[:, ini:end] = dataOut.data
1099 1099
1100 1100 self.__nitems += 1
1101 1101
1102 1102 return int(self.__nitems*nTxs)
1103 1103
1104 1104 def __getBuffer(self):
1105 1105
1106 1106 if self.__nitems == int(1./self.__nTxs):
1107 1107
1108 1108 self.__nitems = 0
1109 1109
1110 1110 return self.__buffer.copy()
1111 1111
1112 1112 return None
1113 1113
1114 1114 def __checkInputs(self, dataOut, shape, nTxs):
1115 1115
1116 1116 if shape is None and nTxs is None:
1117 1117 raise ValueError("Reshaper: shape of factor should be defined")
1118 1118
1119 1119 if nTxs:
1120 1120 if nTxs < 0:
1121 1121 raise ValueError("nTxs should be greater than 0")
1122 1122
1123 1123 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1124 1124 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1125 1125
1126 1126 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1127 1127
1128 1128 return shape, nTxs
1129 1129
1130 1130 if len(shape) != 2 and len(shape) != 3:
1131 1131 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1132 1132
1133 1133 if len(shape) == 2:
1134 1134 shape_tuple = [dataOut.nChannels]
1135 1135 shape_tuple.extend(shape)
1136 1136 else:
1137 1137 shape_tuple = list(shape)
1138 1138
1139 1139 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1140 1140
1141 1141 return shape_tuple, nTxs
1142 1142
1143 1143 def run(self, dataOut, shape=None, nTxs=None):
1144 1144
1145 1145 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1146 1146
1147 1147 dataOut.flagNoData = True
1148 1148 profileIndex = None
1149 1149
1150 1150 if dataOut.flagDataAsBlock:
1151 1151
1152 1152 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1153 1153 dataOut.flagNoData = False
1154 1154
1155 1155 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1156 1156
1157 1157 else:
1158 1158
1159 1159 if self.__nTxs < 1:
1160 1160
1161 1161 self.__appendProfile(dataOut, self.__nTxs)
1162 1162 new_data = self.__getBuffer()
1163 1163
1164 1164 if new_data is not None:
1165 1165 dataOut.data = new_data
1166 1166 dataOut.flagNoData = False
1167 1167
1168 1168 profileIndex = dataOut.profileIndex*nTxs
1169 1169
1170 1170 else:
1171 1171 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1172 1172
1173 1173 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1174 1174
1175 1175 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1176 1176
1177 1177 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1178 1178
1179 1179 dataOut.profileIndex = profileIndex
1180 1180
1181 1181 dataOut.ippSeconds /= self.__nTxs
1182 1182
1183 1183 return dataOut
1184 1184
1185 1185 class SplitProfiles(Operation):
1186 1186
1187 1187 def __init__(self, **kwargs):
1188 1188
1189 1189 Operation.__init__(self, **kwargs)
1190 1190
1191 1191 def run(self, dataOut, n):
1192 1192
1193 1193 dataOut.flagNoData = True
1194 1194 profileIndex = None
1195 1195
1196 1196 if dataOut.flagDataAsBlock:
1197 1197
1198 1198 #nchannels, nprofiles, nsamples
1199 1199 shape = dataOut.data.shape
1200 1200
1201 1201 if shape[2] % n != 0:
1202 1202 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1203 1203
1204 1204 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1205 1205
1206 1206 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1207 1207 dataOut.flagNoData = False
1208 1208
1209 1209 profileIndex = int(dataOut.nProfiles/n) - 1
1210 1210
1211 1211 else:
1212 1212
1213 1213 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1214 1214
1215 1215 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1216 1216
1217 1217 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1218 1218
1219 1219 dataOut.nProfiles = int(dataOut.nProfiles*n)
1220 1220
1221 1221 dataOut.profileIndex = profileIndex
1222 1222
1223 1223 dataOut.ippSeconds /= n
1224 1224
1225 1225 return dataOut
1226 1226
1227 1227 class CombineProfiles(Operation):
1228 1228 def __init__(self, **kwargs):
1229 1229
1230 1230 Operation.__init__(self, **kwargs)
1231 1231
1232 1232 self.__remData = None
1233 1233 self.__profileIndex = 0
1234 1234
1235 1235 def run(self, dataOut, n):
1236 1236
1237 1237 dataOut.flagNoData = True
1238 1238 profileIndex = None
1239 1239
1240 1240 if dataOut.flagDataAsBlock:
1241 1241
1242 1242 #nchannels, nprofiles, nsamples
1243 1243 shape = dataOut.data.shape
1244 1244 new_shape = shape[0], shape[1]/n, shape[2]*n
1245 1245
1246 1246 if shape[1] % n != 0:
1247 1247 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1248 1248
1249 1249 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1250 1250 dataOut.flagNoData = False
1251 1251
1252 1252 profileIndex = int(dataOut.nProfiles*n) - 1
1253 1253
1254 1254 else:
1255 1255
1256 1256 #nchannels, nsamples
1257 1257 if self.__remData is None:
1258 1258 newData = dataOut.data
1259 1259 else:
1260 1260 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1261 1261
1262 1262 self.__profileIndex += 1
1263 1263
1264 1264 if self.__profileIndex < n:
1265 1265 self.__remData = newData
1266 1266 #continue
1267 1267 return
1268 1268
1269 1269 self.__profileIndex = 0
1270 1270 self.__remData = None
1271 1271
1272 1272 dataOut.data = newData
1273 1273 dataOut.flagNoData = False
1274 1274
1275 1275 profileIndex = dataOut.profileIndex/n
1276 1276
1277 1277
1278 1278 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1279 1279
1280 1280 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1281 1281
1282 1282 dataOut.nProfiles = int(dataOut.nProfiles/n)
1283 1283
1284 1284 dataOut.profileIndex = profileIndex
1285 1285
1286 1286 dataOut.ippSeconds *= n
1287 1287
1288 1288 return dataOut
1289 1289
1290 1290 class PulsePair(Operation):
1291 1291 '''
1292 1292 Function PulsePair(Signal Power, Velocity)
1293 1293 The real component of Lag[0] provides Intensity Information
1294 1294 The imag component of Lag[1] Phase provides Velocity Information
1295 1295
1296 1296 Configuration Parameters:
1297 1297 nPRF = Number of Several PRF
1298 1298 theta = Degree Azimuth angel Boundaries
1299 1299
1300 1300 Input:
1301 1301 self.dataOut
1302 1302 lag[N]
1303 1303 Affected:
1304 1304 self.dataOut.spc
1305 1305 '''
1306 1306 isConfig = False
1307 1307 __profIndex = 0
1308 1308 __initime = None
1309 1309 __lastdatatime = None
1310 1310 __buffer = None
1311 1311 noise = None
1312 1312 __dataReady = False
1313 1313 n = None
1314 1314 __nch = 0
1315 1315 __nHeis = 0
1316 1316 removeDC = False
1317 1317 ipp = None
1318 1318 lambda_ = 0
1319 1319
1320 1320 def __init__(self,**kwargs):
1321 1321 Operation.__init__(self,**kwargs)
1322 1322
1323 1323 def setup(self, dataOut, n = None, removeDC=False):
1324 1324 '''
1325 1325 n= Numero de PRF's de entrada
1326 1326 '''
1327 1327 self.__initime = None
1328 1328 ####print("[INICIO]-setup del METODO PULSE PAIR")
1329 1329 self.__lastdatatime = 0
1330 1330 self.__dataReady = False
1331 1331 self.__buffer = 0
1332 1332 self.__profIndex = 0
1333 1333 self.noise = None
1334 1334 self.__nch = dataOut.nChannels
1335 1335 self.__nHeis = dataOut.nHeights
1336 1336 self.removeDC = removeDC
1337 1337 self.lambda_ = 3.0e8/(9345.0e6)
1338 1338 self.ippSec = dataOut.ippSeconds
1339 1339 self.nCohInt = dataOut.nCohInt
1340 1340 ####print("IPPseconds",dataOut.ippSeconds)
1341 1341 ####print("ELVALOR DE n es:", n)
1342 1342 if n == None:
1343 1343 raise ValueError("n should be specified.")
1344 1344
1345 1345 if n != None:
1346 1346 if n<2:
1347 1347 raise ValueError("n should be greater than 2")
1348 1348
1349 1349 self.n = n
1350 1350 self.__nProf = n
1351 1351
1352 1352 self.__buffer = numpy.zeros((dataOut.nChannels,
1353 1353 n,
1354 1354 dataOut.nHeights),
1355 1355 dtype='complex')
1356 1356
1357 1357 def putData(self,data):
1358 1358 '''
1359 1359 Add a profile to he __buffer and increase in one the __profiel Index
1360 1360 '''
1361 1361 self.__buffer[:,self.__profIndex,:]= data
1362 1362 self.__profIndex += 1
1363 1363 return
1364 1364
1365 1365 def pushData(self,dataOut):
1366 1366 '''
1367 1367 Return the PULSEPAIR and the profiles used in the operation
1368 1368 Affected : self.__profileIndex
1369 1369 '''
1370 1370 #----------------- Remove DC-----------------------------------
1371 1371 if self.removeDC==True:
1372 1372 mean = numpy.mean(self.__buffer,1)
1373 1373 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1374 1374 dc= numpy.tile(tmp,[1,self.__nProf,1])
1375 1375 self.__buffer = self.__buffer - dc
1376 1376 #------------------Calculo de Potencia ------------------------
1377 1377 pair0 = self.__buffer*numpy.conj(self.__buffer)
1378 1378 pair0 = pair0.real
1379 1379 lag_0 = numpy.sum(pair0,1)
1380 1380 #-----------------Calculo de Cscp------------------------------ New
1381 1381 cspc_pair01 = self.__buffer[0]*self.__buffer[1]
1382 1382 #------------------Calculo de Ruido x canal--------------------
1383 1383 self.noise = numpy.zeros(self.__nch)
1384 1384 for i in range(self.__nch):
1385 1385 daux = numpy.sort(pair0[i,:,:],axis= None)
1386 1386 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1387 1387
1388 1388 self.noise = self.noise.reshape(self.__nch,1)
1389 1389 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1390 1390 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1391 1391 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1392 1392 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1393 1393 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1394 1394 #-------------------- Power --------------------------------------------------
1395 1395 data_power = lag_0/(self.n*self.nCohInt)
1396 1396 #--------------------CCF------------------------------------------------------
1397 1397 data_ccf =numpy.sum(cspc_pair01,axis=0)/(self.n*self.nCohInt)
1398 1398 #------------------ Senal --------------------------------------------------
1399 1399 data_intensity = pair0 - noise_buffer
1400 1400 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1401 1401 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1402 1402 for i in range(self.__nch):
1403 1403 for j in range(self.__nHeis):
1404 1404 if data_intensity[i][j] < 0:
1405 1405 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1406 1406
1407 1407 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1408 1408 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1409 1409 lag_1 = numpy.sum(pair1,1)
1410 1410 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1411 1411 data_velocity = (self.lambda_/2.0)*data_freq
1412 1412
1413 1413 #---------------- Potencia promedio estimada de la Senal-----------
1414 1414 lag_0 = lag_0/self.n
1415 1415 S = lag_0-self.noise
1416 1416
1417 1417 #---------------- Frecuencia Doppler promedio ---------------------
1418 1418 lag_1 = lag_1/(self.n-1)
1419 1419 R1 = numpy.abs(lag_1)
1420 1420
1421 1421 #---------------- Calculo del SNR----------------------------------
1422 1422 data_snrPP = S/self.noise
1423 1423 for i in range(self.__nch):
1424 1424 for j in range(self.__nHeis):
1425 1425 if data_snrPP[i][j] < 1.e-20:
1426 1426 data_snrPP[i][j] = 1.e-20
1427 1427
1428 1428 #----------------- Calculo del ancho espectral ----------------------
1429 1429 L = S/R1
1430 1430 L = numpy.where(L<0,1,L)
1431 1431 L = numpy.log(L)
1432 1432 tmp = numpy.sqrt(numpy.absolute(L))
1433 1433 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1434 1434 n = self.__profIndex
1435 1435
1436 1436 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1437 1437 self.__profIndex = 0
1438 1438 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,data_ccf,n
1439 1439
1440 1440
1441 1441 def pulsePairbyProfiles(self,dataOut):
1442 1442
1443 1443 self.__dataReady = False
1444 1444 data_power = None
1445 1445 data_intensity = None
1446 1446 data_velocity = None
1447 1447 data_specwidth = None
1448 1448 data_snrPP = None
1449 1449 data_ccf = None
1450 1450 self.putData(data=dataOut.data)
1451 1451 if self.__profIndex == self.n:
1452 1452 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, n = self.pushData(dataOut=dataOut)
1453 1453 self.__dataReady = True
1454 1454
1455 1455 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf
1456 1456
1457 1457
1458 1458 def pulsePairOp(self, dataOut, datatime= None):
1459 1459
1460 1460 if self.__initime == None:
1461 1461 self.__initime = datatime
1462 1462 data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf = self.pulsePairbyProfiles(dataOut)
1463 1463 self.__lastdatatime = datatime
1464 1464
1465 1465 if data_power is None:
1466 1466 return None, None, None,None,None,None,None
1467 1467
1468 1468 avgdatatime = self.__initime
1469 1469 deltatime = datatime - self.__lastdatatime
1470 1470 self.__initime = datatime
1471 1471
1472 1472 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime
1473 1473
1474 1474 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1475 1475 #print("hey")
1476 1476 #print(dataOut.data.shape)
1477 1477 #exit(1)
1478 1478 #print(self.__profIndex)
1479 1479 if not self.isConfig:
1480 1480 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1481 1481 self.isConfig = True
1482 1482 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1483 1483 dataOut.flagNoData = True
1484 1484
1485 1485 if self.__dataReady:
1486 1486 ###print("READY ----------------------------------")
1487 1487 dataOut.nCohInt *= self.n
1488 1488 dataOut.dataPP_POW = data_intensity # S
1489 1489 dataOut.dataPP_POWER = data_power # P valor que corresponde a POTENCIA MOMENTO
1490 1490 dataOut.dataPP_DOP = data_velocity
1491 1491 dataOut.dataPP_SNR = data_snrPP
1492 1492 dataOut.dataPP_WIDTH = data_specwidth
1493 1493 dataOut.dataPP_CCF = data_ccf
1494 1494 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1495 1495 dataOut.nProfiles = int(dataOut.nProfiles/n)
1496 1496 dataOut.utctime = avgdatatime
1497 1497 dataOut.flagNoData = False
1498 1498 return dataOut
1499 1499
1500 1500 class PulsePair_vRF(Operation):
1501 1501 '''
1502 1502 Function PulsePair(Signal Power, Velocity)
1503 1503 The real component of Lag[0] provides Intensity Information
1504 1504 The imag component of Lag[1] Phase provides Velocity Information
1505 1505
1506 1506 Configuration Parameters:
1507 1507 nPRF = Number of Several PRF
1508 1508 theta = Degree Azimuth angel Boundaries
1509 1509
1510 1510 Input:
1511 1511 self.dataOut
1512 1512 lag[N]
1513 1513 Affected:
1514 1514 self.dataOut.spc
1515 1515 '''
1516 1516 isConfig = False
1517 1517 __profIndex = 0
1518 1518 __initime = None
1519 1519 __lastdatatime = None
1520 1520 __buffer = None
1521 1521 noise = None
1522 1522 __dataReady = False
1523 1523 n = None
1524 1524 __nch = 0
1525 1525 __nHeis = 0
1526 1526 removeDC = False
1527 1527 ipp = None
1528 1528 lambda_ = 0
1529 1529
1530 1530 def __init__(self,**kwargs):
1531 1531 Operation.__init__(self,**kwargs)
1532 1532
1533 1533 def setup(self, dataOut, n = None, removeDC=False):
1534 1534 '''
1535 1535 n= Numero de PRF's de entrada
1536 1536 '''
1537 1537 self.__initime = None
1538 1538 ####print("[INICIO]-setup del METODO PULSE PAIR")
1539 1539 self.__lastdatatime = 0
1540 1540 self.__dataReady = False
1541 1541 self.__buffer = 0
1542 1542 self.__profIndex = 0
1543 1543 self.noise = None
1544 1544 self.__nch = dataOut.nChannels
1545 1545 self.__nHeis = dataOut.nHeights
1546 1546 self.removeDC = removeDC
1547 1547 self.lambda_ = 3.0e8/(9345.0e6)
1548 1548 self.ippSec = dataOut.ippSeconds
1549 1549 self.nCohInt = dataOut.nCohInt
1550 1550 ####print("IPPseconds",dataOut.ippSeconds)
1551 1551 ####print("ELVALOR DE n es:", n)
1552 1552 if n == None:
1553 1553 raise ValueError("n should be specified.")
1554 1554
1555 1555 if n != None:
1556 1556 if n<2:
1557 1557 raise ValueError("n should be greater than 2")
1558 1558
1559 1559 self.n = n
1560 1560 self.__nProf = n
1561 1561
1562 1562 self.__buffer = numpy.zeros((dataOut.nChannels,
1563 1563 n,
1564 1564 dataOut.nHeights),
1565 1565 dtype='complex')
1566 1566
1567 1567 def putData(self,data):
1568 1568 '''
1569 1569 Add a profile to he __buffer and increase in one the __profiel Index
1570 1570 '''
1571 1571 self.__buffer[:,self.__profIndex,:]= data
1572 1572 self.__profIndex += 1
1573 1573 return
1574 1574
1575 1575 def putDataByBlock(self,data,n):
1576 1576 '''
1577 1577 Add a profile to he __buffer and increase in one the __profiel Index
1578 1578 '''
1579 1579 self.__buffer[:]= data
1580 1580 self.__profIndex = n
1581 1581 return
1582 1582
1583 1583 def pushData(self,dataOut):
1584 1584 '''
1585 1585 Return the PULSEPAIR and the profiles used in the operation
1586 1586 Affected : self.__profileIndex
1587 1587 '''
1588 1588 #----------------- Remove DC-----------------------------------
1589 1589 if self.removeDC==True:
1590 1590 mean = numpy.mean(self.__buffer,1)
1591 1591 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1592 1592 dc= numpy.tile(tmp,[1,self.__nProf,1])
1593 1593 self.__buffer = self.__buffer - dc
1594 1594 #------------------Calculo de Potencia ------------------------
1595 1595 pair0 = self.__buffer*numpy.conj(self.__buffer)
1596 1596 pair0 = pair0.real
1597 1597 lag_0 = numpy.sum(pair0,1)
1598 1598 #-----------------Calculo de Cscp------------------------------ New
1599 1599 cspc_pair01 = self.__buffer[0]*self.__buffer[1]
1600 1600 #------------------Calculo de Ruido x canal--------------------
1601 1601 self.noise = numpy.zeros(self.__nch)
1602 1602 for i in range(self.__nch):
1603 1603 daux = numpy.sort(pair0[i,:,:],axis= None)
1604 1604 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1605 1605
1606 1606 self.noise = self.noise.reshape(self.__nch,1)
1607 1607 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1608 1608 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1609 1609 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1610 1610 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1611 1611 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1612 1612 #-------------------- Power --------------------------------------------------
1613 1613 data_power = lag_0/(self.n*self.nCohInt)
1614 1614 #--------------------CCF------------------------------------------------------
1615 1615 data_ccf =numpy.sum(cspc_pair01,axis=0)/(self.n*self.nCohInt)
1616 1616 #------------------ Senal --------------------------------------------------
1617 1617 data_intensity = pair0 - noise_buffer
1618 1618 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1619 1619 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1620 1620 for i in range(self.__nch):
1621 1621 for j in range(self.__nHeis):
1622 1622 if data_intensity[i][j] < 0:
1623 1623 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1624 1624
1625 1625 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1626 1626 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1627 1627 lag_1 = numpy.sum(pair1,1)
1628 1628 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1629 1629 data_velocity = (self.lambda_/2.0)*data_freq
1630 1630
1631 1631 #---------------- Potencia promedio estimada de la Senal-----------
1632 1632 lag_0 = lag_0/self.n
1633 1633 S = lag_0-self.noise
1634 1634
1635 1635 #---------------- Frecuencia Doppler promedio ---------------------
1636 1636 lag_1 = lag_1/(self.n-1)
1637 1637 R1 = numpy.abs(lag_1)
1638 1638
1639 1639 #---------------- Calculo del SNR----------------------------------
1640 1640 data_snrPP = S/self.noise
1641 1641 for i in range(self.__nch):
1642 1642 for j in range(self.__nHeis):
1643 1643 if data_snrPP[i][j] < 1.e-20:
1644 1644 data_snrPP[i][j] = 1.e-20
1645 1645
1646 1646 #----------------- Calculo del ancho espectral ----------------------
1647 1647 L = S/R1
1648 1648 L = numpy.where(L<0,1,L)
1649 1649 L = numpy.log(L)
1650 1650 tmp = numpy.sqrt(numpy.absolute(L))
1651 1651 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1652 1652 n = self.__profIndex
1653 1653
1654 1654 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1655 1655 self.__profIndex = 0
1656 1656 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,data_ccf,n
1657 1657
1658 1658
1659 1659 def pulsePairbyProfiles(self,dataOut,n):
1660 1660
1661 1661 self.__dataReady = False
1662 1662 data_power = None
1663 1663 data_intensity = None
1664 1664 data_velocity = None
1665 1665 data_specwidth = None
1666 1666 data_snrPP = None
1667 1667 data_ccf = None
1668 1668
1669 1669 if dataOut.flagDataAsBlock:
1670 1670 self.putDataByBlock(data=dataOut.data,n=n)
1671 1671 else:
1672 1672 self.putData(data=dataOut.data)
1673 1673 if self.__profIndex == self.n:
1674 1674 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, n = self.pushData(dataOut=dataOut)
1675 1675 self.__dataReady = True
1676 1676
1677 1677 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf
1678 1678
1679 1679
1680 1680 def pulsePairOp(self, dataOut, n, datatime= None):
1681 1681
1682 1682 if self.__initime == None:
1683 1683 self.__initime = datatime
1684 1684 data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf = self.pulsePairbyProfiles(dataOut,n)
1685 1685 self.__lastdatatime = datatime
1686 1686
1687 1687 if data_power is None:
1688 1688 return None, None, None,None,None,None,None
1689 1689
1690 1690 avgdatatime = self.__initime
1691 1691 deltatime = datatime - self.__lastdatatime
1692 1692 self.__initime = datatime
1693 1693
1694 1694 return data_power, data_intensity, data_velocity, data_snrPP,data_specwidth,data_ccf, avgdatatime
1695 1695
1696 1696 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1697 1697
1698 1698 if dataOut.flagDataAsBlock:
1699 n = dataOut.nProfiles
1699 n = int(dataOut.nProfiles)
1700 #print("n",n)
1700 1701
1701 1702 if not self.isConfig:
1702 1703 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1703 1704 self.isConfig = True
1704 1705
1705 1706
1706 1707 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth,data_ccf, avgdatatime = self.pulsePairOp(dataOut, n, dataOut.utctime)
1707 1708
1708 1709
1709 1710 dataOut.flagNoData = True
1710 1711
1711 1712 if self.__dataReady:
1712 1713 ###print("READY ----------------------------------")
1713 1714 dataOut.nCohInt *= self.n
1714 1715 dataOut.dataPP_POW = data_intensity # S
1715 1716 dataOut.dataPP_POWER = data_power # P valor que corresponde a POTENCIA MOMENTO
1716 1717 dataOut.dataPP_DOP = data_velocity
1717 1718 dataOut.dataPP_SNR = data_snrPP
1718 1719 dataOut.dataPP_WIDTH = data_specwidth
1719 1720 dataOut.dataPP_CCF = data_ccf
1720 1721 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1721 1722 dataOut.nProfiles = int(dataOut.nProfiles/n)
1722 1723 dataOut.utctime = avgdatatime
1723 1724 dataOut.flagNoData = False
1724 1725 return dataOut
1725 1726
1726 1727 # import collections
1727 1728 # from scipy.stats import mode
1728 1729 #
1729 1730 # class Synchronize(Operation):
1730 1731 #
1731 1732 # isConfig = False
1732 1733 # __profIndex = 0
1733 1734 #
1734 1735 # def __init__(self, **kwargs):
1735 1736 #
1736 1737 # Operation.__init__(self, **kwargs)
1737 1738 # # self.isConfig = False
1738 1739 # self.__powBuffer = None
1739 1740 # self.__startIndex = 0
1740 1741 # self.__pulseFound = False
1741 1742 #
1742 1743 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1743 1744 #
1744 1745 # #Read data
1745 1746 #
1746 1747 # powerdB = dataOut.getPower(channel = channel)
1747 1748 # noisedB = dataOut.getNoise(channel = channel)[0]
1748 1749 #
1749 1750 # self.__powBuffer.extend(powerdB.flatten())
1750 1751 #
1751 1752 # dataArray = numpy.array(self.__powBuffer)
1752 1753 #
1753 1754 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1754 1755 #
1755 1756 # maxValue = numpy.nanmax(filteredPower)
1756 1757 #
1757 1758 # if maxValue < noisedB + 10:
1758 1759 # #No se encuentra ningun pulso de transmision
1759 1760 # return None
1760 1761 #
1761 1762 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1762 1763 #
1763 1764 # if len(maxValuesIndex) < 2:
1764 1765 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1765 1766 # return None
1766 1767 #
1767 1768 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1768 1769 #
1769 1770 # #Seleccionar solo valores con un espaciamiento de nSamples
1770 1771 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1771 1772 #
1772 1773 # if len(pulseIndex) < 2:
1773 1774 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1774 1775 # return None
1775 1776 #
1776 1777 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1777 1778 #
1778 1779 # #remover senales que se distancien menos de 10 unidades o muestras
1779 1780 # #(No deberian existir IPP menor a 10 unidades)
1780 1781 #
1781 1782 # realIndex = numpy.where(spacing > 10 )[0]
1782 1783 #
1783 1784 # if len(realIndex) < 2:
1784 1785 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1785 1786 # return None
1786 1787 #
1787 1788 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1788 1789 # realPulseIndex = pulseIndex[realIndex]
1789 1790 #
1790 1791 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1791 1792 #
1792 1793 # print "IPP = %d samples" %period
1793 1794 #
1794 1795 # self.__newNSamples = dataOut.nHeights #int(period)
1795 1796 # self.__startIndex = int(realPulseIndex[0])
1796 1797 #
1797 1798 # return 1
1798 1799 #
1799 1800 #
1800 1801 # def setup(self, nSamples, nChannels, buffer_size = 4):
1801 1802 #
1802 1803 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1803 1804 # maxlen = buffer_size*nSamples)
1804 1805 #
1805 1806 # bufferList = []
1806 1807 #
1807 1808 # for i in range(nChannels):
1808 1809 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1809 1810 # maxlen = buffer_size*nSamples)
1810 1811 #
1811 1812 # bufferList.append(bufferByChannel)
1812 1813 #
1813 1814 # self.__nSamples = nSamples
1814 1815 # self.__nChannels = nChannels
1815 1816 # self.__bufferList = bufferList
1816 1817 #
1817 1818 # def run(self, dataOut, channel = 0):
1818 1819 #
1819 1820 # if not self.isConfig:
1820 1821 # nSamples = dataOut.nHeights
1821 1822 # nChannels = dataOut.nChannels
1822 1823 # self.setup(nSamples, nChannels)
1823 1824 # self.isConfig = True
1824 1825 #
1825 1826 # #Append new data to internal buffer
1826 1827 # for thisChannel in range(self.__nChannels):
1827 1828 # bufferByChannel = self.__bufferList[thisChannel]
1828 1829 # bufferByChannel.extend(dataOut.data[thisChannel])
1829 1830 #
1830 1831 # if self.__pulseFound:
1831 1832 # self.__startIndex -= self.__nSamples
1832 1833 #
1833 1834 # #Finding Tx Pulse
1834 1835 # if not self.__pulseFound:
1835 1836 # indexFound = self.__findTxPulse(dataOut, channel)
1836 1837 #
1837 1838 # if indexFound == None:
1838 1839 # dataOut.flagNoData = True
1839 1840 # return
1840 1841 #
1841 1842 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1842 1843 # self.__pulseFound = True
1843 1844 # self.__startIndex = indexFound
1844 1845 #
1845 1846 # #If pulse was found ...
1846 1847 # for thisChannel in range(self.__nChannels):
1847 1848 # bufferByChannel = self.__bufferList[thisChannel]
1848 1849 # #print self.__startIndex
1849 1850 # x = numpy.array(bufferByChannel)
1850 1851 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1851 1852 #
1852 1853 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1853 1854 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1854 1855 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1855 1856 #
1856 1857 # dataOut.data = self.__arrayBuffer
1857 1858 #
1858 1859 # self.__startIndex += self.__newNSamples
1859 1860 #
1860 1861 # return
General Comments 0
You need to be logged in to leave comments. Login now