##// END OF EJS Templates
Merge branch 'v3.0-devel' of http://jro-dev.igp.gob.pe/rhodecode/schain into v3.0-devel
Juan C. Espinoza -
r1226:df48418de9d9 merge
parent child
Show More
@@ -1,1363 +1,1365
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10 import json
11 11
12 12 from schainpy.utils import log
13 13 from .jroheaderIO import SystemHeader, RadarControllerHeader
14 14
15 15
16 16 def getNumpyDtype(dataTypeCode):
17 17
18 18 if dataTypeCode == 0:
19 19 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
20 20 elif dataTypeCode == 1:
21 21 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
22 22 elif dataTypeCode == 2:
23 23 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
24 24 elif dataTypeCode == 3:
25 25 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
26 26 elif dataTypeCode == 4:
27 27 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
28 28 elif dataTypeCode == 5:
29 29 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
30 30 else:
31 31 raise ValueError('dataTypeCode was not defined')
32 32
33 33 return numpyDtype
34 34
35 35
36 36 def getDataTypeCode(numpyDtype):
37 37
38 38 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
39 39 datatype = 0
40 40 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
41 41 datatype = 1
42 42 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
43 43 datatype = 2
44 44 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
45 45 datatype = 3
46 46 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
47 47 datatype = 4
48 48 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
49 49 datatype = 5
50 50 else:
51 51 datatype = None
52 52
53 53 return datatype
54 54
55 55
56 56 def hildebrand_sekhon(data, navg):
57 57 """
58 58 This method is for the objective determination of the noise level in Doppler spectra. This
59 59 implementation technique is based on the fact that the standard deviation of the spectral
60 60 densities is equal to the mean spectral density for white Gaussian noise
61 61
62 62 Inputs:
63 63 Data : heights
64 64 navg : numbers of averages
65 65
66 66 Return:
67 67 mean : noise's level
68 68 """
69 69
70 70 sortdata = numpy.sort(data, axis=None)
71 71 lenOfData = len(sortdata)
72 72 nums_min = lenOfData*0.2
73 73
74 74 if nums_min <= 5:
75 75
76 76 nums_min = 5
77 77
78 78 sump = 0.
79 79 sumq = 0.
80 80
81 81 j = 0
82 82 cont = 1
83 83
84 84 while((cont == 1)and(j < lenOfData)):
85 85
86 86 sump += sortdata[j]
87 87 sumq += sortdata[j]**2
88 88
89 89 if j > nums_min:
90 90 rtest = float(j)/(j-1) + 1.0/navg
91 91 if ((sumq*j) > (rtest*sump**2)):
92 92 j = j - 1
93 93 sump = sump - sortdata[j]
94 94 sumq = sumq - sortdata[j]**2
95 95 cont = 0
96 96
97 97 j += 1
98 98
99 99 lnoise = sump / j
100 100
101 101 return lnoise
102 102
103 103
104 104 class Beam:
105 105
106 106 def __init__(self):
107 107 self.codeList = []
108 108 self.azimuthList = []
109 109 self.zenithList = []
110 110
111 111
112 112 class GenericData(object):
113 113
114 114 flagNoData = True
115 115
116 116 def copy(self, inputObj=None):
117 117
118 118 if inputObj == None:
119 119 return copy.deepcopy(self)
120 120
121 121 for key in list(inputObj.__dict__.keys()):
122 122
123 123 attribute = inputObj.__dict__[key]
124 124
125 125 # If this attribute is a tuple or list
126 126 if type(inputObj.__dict__[key]) in (tuple, list):
127 127 self.__dict__[key] = attribute[:]
128 128 continue
129 129
130 130 # If this attribute is another object or instance
131 131 if hasattr(attribute, '__dict__'):
132 132 self.__dict__[key] = attribute.copy()
133 133 continue
134 134
135 135 self.__dict__[key] = inputObj.__dict__[key]
136 136
137 137 def deepcopy(self):
138 138
139 139 return copy.deepcopy(self)
140 140
141 141 def isEmpty(self):
142 142
143 143 return self.flagNoData
144 144
145 145
146 146 class JROData(GenericData):
147 147
148 148 # m_BasicHeader = BasicHeader()
149 149 # m_ProcessingHeader = ProcessingHeader()
150 150
151 151 systemHeaderObj = SystemHeader()
152 152 radarControllerHeaderObj = RadarControllerHeader()
153 153 # data = None
154 154 type = None
155 155 datatype = None # dtype but in string
156 156 # dtype = None
157 157 # nChannels = None
158 158 # nHeights = None
159 159 nProfiles = None
160 160 heightList = None
161 161 channelList = None
162 162 flagDiscontinuousBlock = False
163 163 useLocalTime = False
164 164 utctime = None
165 165 timeZone = None
166 166 dstFlag = None
167 167 errorCount = None
168 168 blocksize = None
169 169 # nCode = None
170 170 # nBaud = None
171 171 # code = None
172 172 flagDecodeData = False # asumo q la data no esta decodificada
173 173 flagDeflipData = False # asumo q la data no esta sin flip
174 174 flagShiftFFT = False
175 175 # ippSeconds = None
176 176 # timeInterval = None
177 177 nCohInt = None
178 178 # noise = None
179 179 windowOfFilter = 1
180 180 # Speed of ligth
181 181 C = 3e8
182 182 frequency = 49.92e6
183 183 realtime = False
184 184 beacon_heiIndexList = None
185 185 last_block = None
186 186 blocknow = None
187 187 azimuth = None
188 188 zenith = None
189 189 beam = Beam()
190 190 profileIndex = None
191 191 error = None
192 192 data = None
193 193 nmodes = None
194 194
195 195 def __str__(self):
196 196
197 197 return '{} - {}'.format(self.type, self.getDatatime())
198 198
199 199 def getNoise(self):
200 200
201 201 raise NotImplementedError
202 202
203 203 def getNChannels(self):
204 204
205 205 return len(self.channelList)
206 206
207 207 def getChannelIndexList(self):
208 208
209 209 return list(range(self.nChannels))
210 210
211 211 def getNHeights(self):
212 212
213 213 return len(self.heightList)
214 214
215 215 def getHeiRange(self, extrapoints=0):
216 216
217 217 heis = self.heightList
218 218 # deltah = self.heightList[1] - self.heightList[0]
219 219 #
220 220 # heis.append(self.heightList[-1])
221 221
222 222 return heis
223 223
224 224 def getDeltaH(self):
225 225
226 226 delta = self.heightList[1] - self.heightList[0]
227 227
228 228 return delta
229 229
230 230 def getltctime(self):
231 231
232 232 if self.useLocalTime:
233 233 return self.utctime - self.timeZone * 60
234 234
235 235 return self.utctime
236 236
237 237 def getDatatime(self):
238 238
239 239 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
240 240 return datatimeValue
241 241
242 242 def getTimeRange(self):
243 243
244 244 datatime = []
245 245
246 246 datatime.append(self.ltctime)
247 247 datatime.append(self.ltctime + self.timeInterval + 1)
248 248
249 249 datatime = numpy.array(datatime)
250 250
251 251 return datatime
252 252
253 253 def getFmaxTimeResponse(self):
254 254
255 255 period = (10**-6) * self.getDeltaH() / (0.15)
256 256
257 257 PRF = 1. / (period * self.nCohInt)
258 258
259 259 fmax = PRF
260 260
261 261 return fmax
262 262
263 263 def getFmax(self):
264 264 PRF = 1. / (self.ippSeconds * self.nCohInt)
265 265
266 266 fmax = PRF
267 267 return fmax
268 268
269 269 def getVmax(self):
270 270
271 271 _lambda = self.C / self.frequency
272 272
273 273 vmax = self.getFmax() * _lambda / 2
274 274
275 275 return vmax
276 276
277 277 def get_ippSeconds(self):
278 278 '''
279 279 '''
280 280 return self.radarControllerHeaderObj.ippSeconds
281 281
282 282 def set_ippSeconds(self, ippSeconds):
283 283 '''
284 284 '''
285 285
286 286 self.radarControllerHeaderObj.ippSeconds = ippSeconds
287 287
288 288 return
289 289
290 290 def get_dtype(self):
291 291 '''
292 292 '''
293 293 return getNumpyDtype(self.datatype)
294 294
295 295 def set_dtype(self, numpyDtype):
296 296 '''
297 297 '''
298 298
299 299 self.datatype = getDataTypeCode(numpyDtype)
300 300
301 301 def get_code(self):
302 302 '''
303 303 '''
304 304 return self.radarControllerHeaderObj.code
305 305
306 306 def set_code(self, code):
307 307 '''
308 308 '''
309 309 self.radarControllerHeaderObj.code = code
310 310
311 311 return
312 312
313 313 def get_ncode(self):
314 314 '''
315 315 '''
316 316 return self.radarControllerHeaderObj.nCode
317 317
318 318 def set_ncode(self, nCode):
319 319 '''
320 320 '''
321 321 self.radarControllerHeaderObj.nCode = nCode
322 322
323 323 return
324 324
325 325 def get_nbaud(self):
326 326 '''
327 327 '''
328 328 return self.radarControllerHeaderObj.nBaud
329 329
330 330 def set_nbaud(self, nBaud):
331 331 '''
332 332 '''
333 333 self.radarControllerHeaderObj.nBaud = nBaud
334 334
335 335 return
336 336
337 337 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
338 338 channelIndexList = property(
339 339 getChannelIndexList, "I'm the 'channelIndexList' property.")
340 340 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
341 341 #noise = property(getNoise, "I'm the 'nHeights' property.")
342 342 datatime = property(getDatatime, "I'm the 'datatime' property")
343 343 ltctime = property(getltctime, "I'm the 'ltctime' property")
344 344 ippSeconds = property(get_ippSeconds, set_ippSeconds)
345 345 dtype = property(get_dtype, set_dtype)
346 346 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
347 347 code = property(get_code, set_code)
348 348 nCode = property(get_ncode, set_ncode)
349 349 nBaud = property(get_nbaud, set_nbaud)
350 350
351 351
352 352 class Voltage(JROData):
353 353
354 354 # data es un numpy array de 2 dmensiones (canales, alturas)
355 355 data = None
356 356
357 357 def __init__(self):
358 358 '''
359 359 Constructor
360 360 '''
361 361
362 362 self.useLocalTime = True
363 363 self.radarControllerHeaderObj = RadarControllerHeader()
364 364 self.systemHeaderObj = SystemHeader()
365 365 self.type = "Voltage"
366 366 self.data = None
367 367 # self.dtype = None
368 368 # self.nChannels = 0
369 369 # self.nHeights = 0
370 370 self.nProfiles = None
371 371 self.heightList = None
372 372 self.channelList = None
373 373 # self.channelIndexList = None
374 374 self.flagNoData = True
375 375 self.flagDiscontinuousBlock = False
376 376 self.utctime = None
377 377 self.timeZone = None
378 378 self.dstFlag = None
379 379 self.errorCount = None
380 380 self.nCohInt = None
381 381 self.blocksize = None
382 382 self.flagDecodeData = False # asumo q la data no esta decodificada
383 383 self.flagDeflipData = False # asumo q la data no esta sin flip
384 384 self.flagShiftFFT = False
385 385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
386 386 self.profileIndex = 0
387 387
388 388 def getNoisebyHildebrand(self, channel=None):
389 389 """
390 390 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
391 391
392 392 Return:
393 393 noiselevel
394 394 """
395 395
396 396 if channel != None:
397 397 data = self.data[channel]
398 398 nChannels = 1
399 399 else:
400 400 data = self.data
401 401 nChannels = self.nChannels
402 402
403 403 noise = numpy.zeros(nChannels)
404 404 power = data * numpy.conjugate(data)
405 405
406 406 for thisChannel in range(nChannels):
407 407 if nChannels == 1:
408 408 daux = power[:].real
409 409 else:
410 410 daux = power[thisChannel, :].real
411 411 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
412 412
413 413 return noise
414 414
415 415 def getNoise(self, type=1, channel=None):
416 416
417 417 if type == 1:
418 418 noise = self.getNoisebyHildebrand(channel)
419 419
420 420 return noise
421 421
422 422 def getPower(self, channel=None):
423 423
424 424 if channel != None:
425 425 data = self.data[channel]
426 426 else:
427 427 data = self.data
428 428
429 429 power = data * numpy.conjugate(data)
430 430 powerdB = 10 * numpy.log10(power.real)
431 431 powerdB = numpy.squeeze(powerdB)
432 432
433 433 return powerdB
434 434
435 435 def getTimeInterval(self):
436 436
437 437 timeInterval = self.ippSeconds * self.nCohInt
438 438
439 439 return timeInterval
440 440
441 441 noise = property(getNoise, "I'm the 'nHeights' property.")
442 442 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
443 443
444 444
445 445 class Spectra(JROData):
446 446
447 447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
448 448 data_spc = None
449 449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
450 450 data_cspc = None
451 451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
452 452 data_dc = None
453 453 # data power
454 454 data_pwr = None
455 455 nFFTPoints = None
456 456 # nPairs = None
457 457 pairsList = None
458 458 nIncohInt = None
459 459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
460 460 nCohInt = None # se requiere para determinar el valor de timeInterval
461 461 ippFactor = None
462 462 profileIndex = 0
463 463 plotting = "spectra"
464 464
465 465 def __init__(self):
466 466 '''
467 467 Constructor
468 468 '''
469 469
470 470 self.useLocalTime = True
471 471 self.radarControllerHeaderObj = RadarControllerHeader()
472 472 self.systemHeaderObj = SystemHeader()
473 473 self.type = "Spectra"
474 474 # self.data = None
475 475 # self.dtype = None
476 476 # self.nChannels = 0
477 477 # self.nHeights = 0
478 478 self.nProfiles = None
479 479 self.heightList = None
480 480 self.channelList = None
481 481 # self.channelIndexList = None
482 482 self.pairsList = None
483 483 self.flagNoData = True
484 484 self.flagDiscontinuousBlock = False
485 485 self.utctime = None
486 486 self.nCohInt = None
487 487 self.nIncohInt = None
488 488 self.blocksize = None
489 489 self.nFFTPoints = None
490 490 self.wavelength = None
491 491 self.flagDecodeData = False # asumo q la data no esta decodificada
492 492 self.flagDeflipData = False # asumo q la data no esta sin flip
493 493 self.flagShiftFFT = False
494 494 self.ippFactor = 1
495 495 #self.noise = None
496 496 self.beacon_heiIndexList = []
497 497 self.noise_estimation = None
498 498
499 499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
500 500 """
501 501 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
502 502
503 503 Return:
504 504 noiselevel
505 505 """
506 506
507 507 noise = numpy.zeros(self.nChannels)
508 508
509 509 for channel in range(self.nChannels):
510 510 daux = self.data_spc[channel,
511 511 xmin_index:xmax_index, ymin_index:ymax_index]
512 512 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
513 513
514 514 return noise
515 515
516 516 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
517 517
518 518 if self.noise_estimation is not None:
519 519 # this was estimated by getNoise Operation defined in jroproc_spectra.py
520 520 return self.noise_estimation
521 521 else:
522 522 noise = self.getNoisebyHildebrand(
523 523 xmin_index, xmax_index, ymin_index, ymax_index)
524 524 return noise
525 525
526 526 def getFreqRangeTimeResponse(self, extrapoints=0):
527 527
528 528 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
529 529 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
530 530
531 531 return freqrange
532 532
533 533 def getAcfRange(self, extrapoints=0):
534 534
535 535 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
536 536 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
537 537
538 538 return freqrange
539 539
540 540 def getFreqRange(self, extrapoints=0):
541 541
542 542 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
543 543 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
544 544
545 545 return freqrange
546 546
547 547 def getVelRange(self, extrapoints=0):
548 548
549 549 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
550 550 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
551 551
552 552 if self.nmodes:
553 553 return velrange/self.nmodes
554 554 else:
555 555 return velrange
556 556
557 557 def getNPairs(self):
558 558
559 559 return len(self.pairsList)
560 560
561 561 def getPairsIndexList(self):
562 562
563 563 return list(range(self.nPairs))
564 564
565 565 def getNormFactor(self):
566 566
567 567 pwcode = 1
568 568
569 569 if self.flagDecodeData:
570 570 pwcode = numpy.sum(self.code[0]**2)
571 571 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
572 572 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
573 573
574 574 return normFactor
575 575
576 576 def getFlagCspc(self):
577 577
578 578 if self.data_cspc is None:
579 579 return True
580 580
581 581 return False
582 582
583 583 def getFlagDc(self):
584 584
585 585 if self.data_dc is None:
586 586 return True
587 587
588 588 return False
589 589
590 590 def getTimeInterval(self):
591 591
592 592 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
593 593 if self.nmodes:
594 594 return self.nmodes*timeInterval
595 595 else:
596 596 return timeInterval
597 597
598 598 def getPower(self):
599 599
600 600 factor = self.normFactor
601 601 z = self.data_spc / factor
602 602 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
603 603 avg = numpy.average(z, axis=1)
604 604
605 605 return 10 * numpy.log10(avg)
606 606
607 607 def getCoherence(self, pairsList=None, phase=False):
608 608
609 609 z = []
610 610 if pairsList is None:
611 611 pairsIndexList = self.pairsIndexList
612 612 else:
613 613 pairsIndexList = []
614 614 for pair in pairsList:
615 615 if pair not in self.pairsList:
616 616 raise ValueError("Pair %s is not in dataOut.pairsList" % (
617 617 pair))
618 618 pairsIndexList.append(self.pairsList.index(pair))
619 619 for i in range(len(pairsIndexList)):
620 620 pair = self.pairsList[pairsIndexList[i]]
621 621 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
622 622 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
623 623 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
624 624 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
625 625 if phase:
626 626 data = numpy.arctan2(avgcoherenceComplex.imag,
627 627 avgcoherenceComplex.real) * 180 / numpy.pi
628 628 else:
629 629 data = numpy.abs(avgcoherenceComplex)
630 630
631 631 z.append(data)
632 632
633 633 return numpy.array(z)
634 634
635 635 def setValue(self, value):
636 636
637 637 print("This property should not be initialized")
638 638
639 639 return
640 640
641 641 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
642 642 pairsIndexList = property(
643 643 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
644 644 normFactor = property(getNormFactor, setValue,
645 645 "I'm the 'getNormFactor' property.")
646 646 flag_cspc = property(getFlagCspc, setValue)
647 647 flag_dc = property(getFlagDc, setValue)
648 648 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
649 649 timeInterval = property(getTimeInterval, setValue,
650 650 "I'm the 'timeInterval' property")
651 651
652 652
653 653 class SpectraHeis(Spectra):
654 654
655 655 data_spc = None
656 656 data_cspc = None
657 657 data_dc = None
658 658 nFFTPoints = None
659 659 # nPairs = None
660 660 pairsList = None
661 661 nCohInt = None
662 662 nIncohInt = None
663 663
664 664 def __init__(self):
665 665
666 666 self.radarControllerHeaderObj = RadarControllerHeader()
667 667
668 668 self.systemHeaderObj = SystemHeader()
669 669
670 670 self.type = "SpectraHeis"
671 671
672 672 # self.dtype = None
673 673
674 674 # self.nChannels = 0
675 675
676 676 # self.nHeights = 0
677 677
678 678 self.nProfiles = None
679 679
680 680 self.heightList = None
681 681
682 682 self.channelList = None
683 683
684 684 # self.channelIndexList = None
685 685
686 686 self.flagNoData = True
687 687
688 688 self.flagDiscontinuousBlock = False
689 689
690 690 # self.nPairs = 0
691 691
692 692 self.utctime = None
693 693
694 694 self.blocksize = None
695 695
696 696 self.profileIndex = 0
697 697
698 698 self.nCohInt = 1
699 699
700 700 self.nIncohInt = 1
701 701
702 702 def getNormFactor(self):
703 703 pwcode = 1
704 704 if self.flagDecodeData:
705 705 pwcode = numpy.sum(self.code[0]**2)
706 706
707 707 normFactor = self.nIncohInt * self.nCohInt * pwcode
708 708
709 709 return normFactor
710 710
711 711 def getTimeInterval(self):
712 712
713 713 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
714 714
715 715 return timeInterval
716 716
717 717 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
718 718 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
719 719
720 720
721 721 class Fits(JROData):
722 722
723 723 heightList = None
724 724 channelList = None
725 725 flagNoData = True
726 726 flagDiscontinuousBlock = False
727 727 useLocalTime = False
728 728 utctime = None
729 729 timeZone = None
730 730 # ippSeconds = None
731 731 # timeInterval = None
732 732 nCohInt = None
733 733 nIncohInt = None
734 734 noise = None
735 735 windowOfFilter = 1
736 736 # Speed of ligth
737 737 C = 3e8
738 738 frequency = 49.92e6
739 739 realtime = False
740 740
741 741 def __init__(self):
742 742
743 743 self.type = "Fits"
744 744
745 745 self.nProfiles = None
746 746
747 747 self.heightList = None
748 748
749 749 self.channelList = None
750 750
751 751 # self.channelIndexList = None
752 752
753 753 self.flagNoData = True
754 754
755 755 self.utctime = None
756 756
757 757 self.nCohInt = 1
758 758
759 759 self.nIncohInt = 1
760 760
761 761 self.useLocalTime = True
762 762
763 763 self.profileIndex = 0
764 764
765 765 # self.utctime = None
766 766 # self.timeZone = None
767 767 # self.ltctime = None
768 768 # self.timeInterval = None
769 769 # self.header = None
770 770 # self.data_header = None
771 771 # self.data = None
772 772 # self.datatime = None
773 773 # self.flagNoData = False
774 774 # self.expName = ''
775 775 # self.nChannels = None
776 776 # self.nSamples = None
777 777 # self.dataBlocksPerFile = None
778 778 # self.comments = ''
779 779 #
780 780
781 781 def getltctime(self):
782 782
783 783 if self.useLocalTime:
784 784 return self.utctime - self.timeZone * 60
785 785
786 786 return self.utctime
787 787
788 788 def getDatatime(self):
789 789
790 790 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
791 791 return datatime
792 792
793 793 def getTimeRange(self):
794 794
795 795 datatime = []
796 796
797 797 datatime.append(self.ltctime)
798 798 datatime.append(self.ltctime + self.timeInterval)
799 799
800 800 datatime = numpy.array(datatime)
801 801
802 802 return datatime
803 803
804 804 def getHeiRange(self):
805 805
806 806 heis = self.heightList
807 807
808 808 return heis
809 809
810 810 def getNHeights(self):
811 811
812 812 return len(self.heightList)
813 813
814 814 def getNChannels(self):
815 815
816 816 return len(self.channelList)
817 817
818 818 def getChannelIndexList(self):
819 819
820 820 return list(range(self.nChannels))
821 821
822 822 def getNoise(self, type=1):
823 823
824 824 #noise = numpy.zeros(self.nChannels)
825 825
826 826 if type == 1:
827 827 noise = self.getNoisebyHildebrand()
828 828
829 829 if type == 2:
830 830 noise = self.getNoisebySort()
831 831
832 832 if type == 3:
833 833 noise = self.getNoisebyWindow()
834 834
835 835 return noise
836 836
837 837 def getTimeInterval(self):
838 838
839 839 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
840 840
841 841 return timeInterval
842 842
843 843 def get_ippSeconds(self):
844 844 '''
845 845 '''
846 846 return self.ipp_sec
847 847
848 848
849 849 datatime = property(getDatatime, "I'm the 'datatime' property")
850 850 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
851 851 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
852 852 channelIndexList = property(
853 853 getChannelIndexList, "I'm the 'channelIndexList' property.")
854 854 noise = property(getNoise, "I'm the 'nHeights' property.")
855 855
856 856 ltctime = property(getltctime, "I'm the 'ltctime' property")
857 857 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
858 858 ippSeconds = property(get_ippSeconds, '')
859 859
860 860 class Correlation(JROData):
861 861
862 862 noise = None
863 863 SNR = None
864 864 #--------------------------------------------------
865 865 mode = None
866 866 split = False
867 867 data_cf = None
868 868 lags = None
869 869 lagRange = None
870 870 pairsList = None
871 871 normFactor = None
872 872 #--------------------------------------------------
873 873 # calculateVelocity = None
874 874 nLags = None
875 875 nPairs = None
876 876 nAvg = None
877 877
878 878 def __init__(self):
879 879 '''
880 880 Constructor
881 881 '''
882 882 self.radarControllerHeaderObj = RadarControllerHeader()
883 883
884 884 self.systemHeaderObj = SystemHeader()
885 885
886 886 self.type = "Correlation"
887 887
888 888 self.data = None
889 889
890 890 self.dtype = None
891 891
892 892 self.nProfiles = None
893 893
894 894 self.heightList = None
895 895
896 896 self.channelList = None
897 897
898 898 self.flagNoData = True
899 899
900 900 self.flagDiscontinuousBlock = False
901 901
902 902 self.utctime = None
903 903
904 904 self.timeZone = None
905 905
906 906 self.dstFlag = None
907 907
908 908 self.errorCount = None
909 909
910 910 self.blocksize = None
911 911
912 912 self.flagDecodeData = False # asumo q la data no esta decodificada
913 913
914 914 self.flagDeflipData = False # asumo q la data no esta sin flip
915 915
916 916 self.pairsList = None
917 917
918 918 self.nPoints = None
919 919
920 920 def getPairsList(self):
921 921
922 922 return self.pairsList
923 923
924 924 def getNoise(self, mode=2):
925 925
926 926 indR = numpy.where(self.lagR == 0)[0][0]
927 927 indT = numpy.where(self.lagT == 0)[0][0]
928 928
929 929 jspectra0 = self.data_corr[:, :, indR, :]
930 930 jspectra = copy.copy(jspectra0)
931 931
932 932 num_chan = jspectra.shape[0]
933 933 num_hei = jspectra.shape[2]
934 934
935 935 freq_dc = jspectra.shape[1] / 2
936 936 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
937 937
938 938 if ind_vel[0] < 0:
939 939 ind_vel[list(range(0, 1))] = ind_vel[list(
940 940 range(0, 1))] + self.num_prof
941 941
942 942 if mode == 1:
943 943 jspectra[:, freq_dc, :] = (
944 944 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
945 945
946 946 if mode == 2:
947 947
948 948 vel = numpy.array([-2, -1, 1, 2])
949 949 xx = numpy.zeros([4, 4])
950 950
951 951 for fil in range(4):
952 952 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
953 953
954 954 xx_inv = numpy.linalg.inv(xx)
955 955 xx_aux = xx_inv[0, :]
956 956
957 957 for ich in range(num_chan):
958 958 yy = jspectra[ich, ind_vel, :]
959 959 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
960 960
961 961 junkid = jspectra[ich, freq_dc, :] <= 0
962 962 cjunkid = sum(junkid)
963 963
964 964 if cjunkid.any():
965 965 jspectra[ich, freq_dc, junkid.nonzero()] = (
966 966 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
967 967
968 968 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
969 969
970 970 return noise
971 971
972 972 def getTimeInterval(self):
973 973
974 974 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
975 975
976 976 return timeInterval
977 977
978 978 def splitFunctions(self):
979 979
980 980 pairsList = self.pairsList
981 981 ccf_pairs = []
982 982 acf_pairs = []
983 983 ccf_ind = []
984 984 acf_ind = []
985 985 for l in range(len(pairsList)):
986 986 chan0 = pairsList[l][0]
987 987 chan1 = pairsList[l][1]
988 988
989 989 # Obteniendo pares de Autocorrelacion
990 990 if chan0 == chan1:
991 991 acf_pairs.append(chan0)
992 992 acf_ind.append(l)
993 993 else:
994 994 ccf_pairs.append(pairsList[l])
995 995 ccf_ind.append(l)
996 996
997 997 data_acf = self.data_cf[acf_ind]
998 998 data_ccf = self.data_cf[ccf_ind]
999 999
1000 1000 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1001 1001
1002 1002 def getNormFactor(self):
1003 1003 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1004 1004 acf_pairs = numpy.array(acf_pairs)
1005 1005 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1006 1006
1007 1007 for p in range(self.nPairs):
1008 1008 pair = self.pairsList[p]
1009 1009
1010 1010 ch0 = pair[0]
1011 1011 ch1 = pair[1]
1012 1012
1013 1013 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1014 1014 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1015 1015 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1016 1016
1017 1017 return normFactor
1018 1018
1019 1019 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1020 1020 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1021 1021
1022 1022
1023 1023 class Parameters(Spectra):
1024 1024
1025 1025 experimentInfo = None # Information about the experiment
1026 1026 # Information from previous data
1027 1027 inputUnit = None # Type of data to be processed
1028 1028 operation = None # Type of operation to parametrize
1029 1029 # normFactor = None #Normalization Factor
1030 1030 groupList = None # List of Pairs, Groups, etc
1031 1031 # Parameters
1032 1032 data_param = None # Parameters obtained
1033 1033 data_pre = None # Data Pre Parametrization
1034 1034 data_SNR = None # Signal to Noise Ratio
1035 1035 # heightRange = None #Heights
1036 1036 abscissaList = None # Abscissa, can be velocities, lags or time
1037 1037 # noise = None #Noise Potency
1038 1038 utctimeInit = None # Initial UTC time
1039 1039 paramInterval = None # Time interval to calculate Parameters in seconds
1040 1040 useLocalTime = True
1041 1041 # Fitting
1042 1042 data_error = None # Error of the estimation
1043 1043 constants = None
1044 1044 library = None
1045 1045 # Output signal
1046 1046 outputInterval = None # Time interval to calculate output signal in seconds
1047 1047 data_output = None # Out signal
1048 1048 nAvg = None
1049 1049 noise_estimation = None
1050 1050 GauSPC = None # Fit gaussian SPC
1051 1051
1052 1052 def __init__(self):
1053 1053 '''
1054 1054 Constructor
1055 1055 '''
1056 1056 self.radarControllerHeaderObj = RadarControllerHeader()
1057 1057
1058 1058 self.systemHeaderObj = SystemHeader()
1059 1059
1060 1060 self.type = "Parameters"
1061 1061
1062 1062 def getTimeRange1(self, interval):
1063 1063
1064 1064 datatime = []
1065 1065
1066 1066 if self.useLocalTime:
1067 1067 time1 = self.utctimeInit - self.timeZone * 60
1068 1068 else:
1069 1069 time1 = self.utctimeInit
1070 1070
1071 1071 datatime.append(time1)
1072 1072 datatime.append(time1 + interval)
1073 1073 datatime = numpy.array(datatime)
1074 1074
1075 1075 return datatime
1076 1076
1077 1077 def getTimeInterval(self):
1078 1078
1079 1079 if hasattr(self, 'timeInterval1'):
1080 1080 return self.timeInterval1
1081 1081 else:
1082 1082 return self.paramInterval
1083 1083
1084 1084 def setValue(self, value):
1085 1085
1086 1086 print("This property should not be initialized")
1087 1087
1088 1088 return
1089 1089
1090 1090 def getNoise(self):
1091 1091
1092 1092 return self.spc_noise
1093 1093
1094 1094 timeInterval = property(getTimeInterval)
1095 1095 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1096 1096
1097 1097
1098 1098 class PlotterData(object):
1099 1099 '''
1100 1100 Object to hold data to be plotted
1101 1101 '''
1102 1102
1103 1103 MAXNUMX = 100
1104 1104 MAXNUMY = 100
1105 1105
1106 def __init__(self, code, throttle_value, exp_code, buffering=True):
1106 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1107 1107
1108 self.key = code
1108 1109 self.throttle = throttle_value
1109 1110 self.exp_code = exp_code
1110 1111 self.buffering = buffering
1111 1112 self.ready = False
1112 1113 self.localtime = False
1113 1114 self.data = {}
1114 1115 self.meta = {}
1115 1116 self.__times = []
1116 1117 self.__heights = []
1117 1118
1118 1119 if 'snr' in code:
1119 1120 self.plottypes = ['snr']
1120 1121 elif code == 'spc':
1121 1122 self.plottypes = ['spc', 'noise', 'rti']
1122 1123 elif code == 'rti':
1123 1124 self.plottypes = ['noise', 'rti']
1124 1125 else:
1125 1126 self.plottypes = [code]
1126
1127
1128 if 'snr' not in self.plottypes and snr:
1129 self.plottypes.append('snr')
1130
1127 1131 for plot in self.plottypes:
1128 1132 self.data[plot] = {}
1129 1133
1130 1134 def __str__(self):
1131 1135 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1132 1136 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1133 1137
1134 1138 def __len__(self):
1135 1139 return len(self.__times)
1136 1140
1137 1141 def __getitem__(self, key):
1138 1142
1139 1143 if key not in self.data:
1140 1144 raise KeyError(log.error('Missing key: {}'.format(key)))
1141 1145 if 'spc' in key or not self.buffering:
1142 1146 ret = self.data[key]
1143 1147 elif 'scope' in key:
1144 1148 ret = numpy.array(self.data[key][float(self.tm)])
1145 1149 else:
1146 1150 ret = numpy.array([self.data[key][x] for x in self.times])
1147 1151 if ret.ndim > 1:
1148 1152 ret = numpy.swapaxes(ret, 0, 1)
1149 1153 return ret
1150 1154
1151 1155 def __contains__(self, key):
1152 1156 return key in self.data
1153 1157
1154 1158 def setup(self):
1155 1159 '''
1156 1160 Configure object
1157 1161 '''
1158 1162
1159 1163 self.type = ''
1160 1164 self.ready = False
1161 1165 self.data = {}
1162 1166 self.__times = []
1163 1167 self.__heights = []
1164 1168 self.__all_heights = set()
1165 1169 for plot in self.plottypes:
1166 1170 if 'snr' in plot:
1167 1171 plot = 'snr'
1168 1172 elif 'spc_moments' == plot:
1169 1173 plot = 'moments'
1170 1174 self.data[plot] = {}
1171 1175
1172 1176 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1173 1177 self.data['noise'] = {}
1174 1178 self.data['rti'] = {}
1175 1179 if 'noise' not in self.plottypes:
1176 1180 self.plottypes.append('noise')
1177 1181 if 'rti' not in self.plottypes:
1178 1182 self.plottypes.append('rti')
1179 1183
1180 1184 def shape(self, key):
1181 1185 '''
1182 1186 Get the shape of the one-element data for the given key
1183 1187 '''
1184 1188
1185 1189 if len(self.data[key]):
1186 1190 if 'spc' in key or not self.buffering:
1187 1191 return self.data[key].shape
1188 1192 return self.data[key][self.__times[0]].shape
1189 1193 return (0,)
1190 1194
1191 1195 def update(self, dataOut, tm):
1192 1196 '''
1193 1197 Update data object with new dataOut
1194 1198 '''
1195 1199
1196 1200 if tm in self.__times:
1197 1201 return
1198 1202 self.profileIndex = dataOut.profileIndex
1199 1203 self.tm = tm
1200 1204 self.type = dataOut.type
1201 1205 self.parameters = getattr(dataOut, 'parameters', [])
1202 if hasattr(dataOut, 'pairsList'):
1203 self.pairs = dataOut.pairsList
1204 1206 if hasattr(dataOut, 'meta'):
1205 self.meta = dataOut.meta
1207 self.meta.update(dataOut.meta)
1206 1208 self.channels = dataOut.channelList
1207 1209 self.interval = dataOut.getTimeInterval()
1208 1210 self.localtime = dataOut.useLocalTime
1209 1211 if 'spc' in self.plottypes or 'cspc' in self.plottypes or 'spc_moments' in self.plottypes:
1210 1212 self.xrange = (dataOut.getFreqRange(1)/1000.,
1211 1213 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1212 1214 self.factor = dataOut.normFactor
1213 1215 self.__heights.append(dataOut.heightList)
1214 1216 self.__all_heights.update(dataOut.heightList)
1215 1217 self.__times.append(tm)
1216 1218
1217 1219 for plot in self.plottypes:
1218 1220 if plot in ('spc', 'spc_moments'):
1219 1221 z = dataOut.data_spc/dataOut.normFactor
1220 1222 buffer = 10*numpy.log10(z)
1221 1223 if plot == 'cspc':
1222 1224 z = dataOut.data_spc/dataOut.normFactor
1223 1225 buffer = (dataOut.data_spc, dataOut.data_cspc)
1224 1226 if plot == 'noise':
1225 1227 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1226 1228 if plot == 'rti':
1227 1229 buffer = dataOut.getPower()
1228 1230 if plot == 'snr_db':
1229 1231 buffer = dataOut.data_SNR
1230 1232 if plot == 'snr':
1231 1233 buffer = 10*numpy.log10(dataOut.data_SNR)
1232 1234 if plot == 'dop':
1233 1235 buffer = 10*numpy.log10(dataOut.data_DOP)
1234 1236 if plot == 'mean':
1235 1237 buffer = dataOut.data_MEAN
1236 1238 if plot == 'std':
1237 1239 buffer = dataOut.data_STD
1238 1240 if plot == 'coh':
1239 1241 buffer = dataOut.getCoherence()
1240 1242 if plot == 'phase':
1241 1243 buffer = dataOut.getCoherence(phase=True)
1242 1244 if plot == 'output':
1243 1245 buffer = dataOut.data_output
1244 1246 if plot == 'param':
1245 1247 buffer = dataOut.data_param
1246 1248 if plot == 'scope':
1247 1249 buffer = dataOut.data
1248 1250 self.flagDataAsBlock = dataOut.flagDataAsBlock
1249 1251 self.nProfiles = dataOut.nProfiles
1250 1252
1251 1253 if plot == 'spc':
1252 1254 self.data['spc'] = buffer
1253 1255 elif plot == 'cspc':
1254 1256 self.data['spc'] = buffer[0]
1255 1257 self.data['cspc'] = buffer[1]
1256 1258 elif plot == 'spc_moments':
1257 1259 self.data['spc'] = buffer
1258 1260 self.data['moments'][tm] = dataOut.moments
1259 1261 else:
1260 1262 if self.buffering:
1261 1263 self.data[plot][tm] = buffer
1262 1264 else:
1263 1265 self.data[plot] = buffer
1264 1266
1265 1267 def normalize_heights(self):
1266 1268 '''
1267 1269 Ensure same-dimension of the data for different heighList
1268 1270 '''
1269 1271
1270 1272 H = numpy.array(list(self.__all_heights))
1271 1273 H.sort()
1272 1274 for key in self.data:
1273 1275 shape = self.shape(key)[:-1] + H.shape
1274 1276 for tm, obj in list(self.data[key].items()):
1275 1277 h = self.__heights[self.__times.index(tm)]
1276 1278 if H.size == h.size:
1277 1279 continue
1278 1280 index = numpy.where(numpy.in1d(H, h))[0]
1279 1281 dummy = numpy.zeros(shape) + numpy.nan
1280 1282 if len(shape) == 2:
1281 1283 dummy[:, index] = obj
1282 1284 else:
1283 1285 dummy[index] = obj
1284 1286 self.data[key][tm] = dummy
1285 1287
1286 1288 self.__heights = [H for tm in self.__times]
1287 1289
1288 def jsonify(self, decimate=False):
1290 def jsonify(self, plot_name, plot_type, decimate=False):
1289 1291 '''
1290 1292 Convert data to json
1291 1293 '''
1292 1294
1293 data = {}
1294 1295 tm = self.times[-1]
1295 1296 dy = int(self.heights.size/self.MAXNUMY) + 1
1296 for key in self.data:
1297 if key in ('spc', 'cspc') or not self.buffering:
1298 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1299 data[key] = self.roundFloats(
1300 self.data[key][::, ::dx, ::dy].tolist())
1301 else:
1302 data[key] = self.roundFloats(self.data[key][tm].tolist())
1303
1304 ret = {'data': data}
1305 ret['exp_code'] = self.exp_code
1306 ret['time'] = float(tm)
1307 ret['interval'] = float(self.interval)
1308 ret['localtime'] = self.localtime
1309 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1310 if 'spc' in self.data or 'cspc' in self.data:
1311 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1297 if self.key in ('spc', 'cspc') or not self.buffering:
1298 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1299 data = self.roundFloats(
1300 self.data[self.key][::, ::dx, ::dy].tolist())
1312 1301 else:
1313 ret['xrange'] = []
1314 if hasattr(self, 'pairs'):
1315 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1302 data = self.roundFloats(self.data[self.key][tm].tolist())
1303 if self.key is 'noise':
1304 data = [[x] for x in data]
1305
1306 meta = {}
1307 ret = {
1308 'plot': plot_name,
1309 'code': self.exp_code,
1310 'time': float(tm),
1311 'data': data,
1312 }
1313 meta['type'] = plot_type
1314 meta['interval'] = float(self.interval)
1315 meta['localtime'] = self.localtime
1316 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1317 if 'spc' in self.data or 'cspc' in self.data:
1318 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1316 1319 else:
1317 ret['pairs'] = []
1318
1319 for key, value in list(self.meta.items()):
1320 ret[key] = value
1320 meta['xrange'] = []
1321 1321
1322 meta.update(self.meta)
1323 ret['metadata'] = meta
1322 1324 return json.dumps(ret)
1323 1325
1324 1326 @property
1325 1327 def times(self):
1326 1328 '''
1327 1329 Return the list of times of the current data
1328 1330 '''
1329 1331
1330 1332 ret = numpy.array(self.__times)
1331 1333 ret.sort()
1332 1334 return ret
1333 1335
1334 1336 @property
1335 1337 def min_time(self):
1336 1338 '''
1337 1339 Return the minimun time value
1338 1340 '''
1339 1341
1340 1342 return self.times[0]
1341 1343
1342 1344 @property
1343 1345 def max_time(self):
1344 1346 '''
1345 1347 Return the maximun time value
1346 1348 '''
1347 1349
1348 1350 return self.times[-1]
1349 1351
1350 1352 @property
1351 1353 def heights(self):
1352 1354 '''
1353 1355 Return the list of heights of the current data
1354 1356 '''
1355 1357
1356 1358 return numpy.array(self.__heights[-1])
1357 1359
1358 1360 @staticmethod
1359 1361 def roundFloats(obj):
1360 1362 if isinstance(obj, list):
1361 1363 return list(map(PlotterData.roundFloats, obj))
1362 1364 elif isinstance(obj, float):
1363 1365 return round(obj, 2)
@@ -1,784 +1,808
1 1
2 2 import os
3 3 import sys
4 4 import zmq
5 5 import time
6 6 import numpy
7 7 import datetime
8 8 from functools import wraps
9 9 from threading import Thread
10 10 import matplotlib
11 11
12 12 if 'BACKEND' in os.environ:
13 13 matplotlib.use(os.environ['BACKEND'])
14 14 elif 'linux' in sys.platform:
15 15 matplotlib.use("TkAgg")
16 16 elif 'darwin' in sys.platform:
17 17 matplotlib.use('WxAgg')
18 18 else:
19 19 from schainpy.utils import log
20 20 log.warning('Using default Backend="Agg"', 'INFO')
21 21 matplotlib.use('Agg')
22 22
23 23 import matplotlib.pyplot as plt
24 24 from matplotlib.patches import Polygon
25 25 from mpl_toolkits.axes_grid1 import make_axes_locatable
26 26 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
27 27
28 28 from schainpy.model.data.jrodata import PlotterData
29 29 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
30 30 from schainpy.utils import log
31 31
32 32 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
33 33 blu_values = matplotlib.pyplot.get_cmap(
34 34 'seismic_r', 20)(numpy.arange(20))[10:15]
35 35 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
36 36 'jro', numpy.vstack((blu_values, jet_values)))
37 37 matplotlib.pyplot.register_cmap(cmap=ncmap)
38 38
39 39 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
40 40 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
41 41
42 42 EARTH_RADIUS = 6.3710e3
43 43
44 44 def ll2xy(lat1, lon1, lat2, lon2):
45 45
46 46 p = 0.017453292519943295
47 47 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
48 48 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
49 49 r = 12742 * numpy.arcsin(numpy.sqrt(a))
50 50 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
51 51 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
52 52 theta = -theta + numpy.pi/2
53 53 return r*numpy.cos(theta), r*numpy.sin(theta)
54 54
55 55
56 56 def km2deg(km):
57 57 '''
58 58 Convert distance in km to degrees
59 59 '''
60 60
61 61 return numpy.rad2deg(km/EARTH_RADIUS)
62 62
63 63
64 64 def figpause(interval):
65 65 backend = plt.rcParams['backend']
66 66 if backend in matplotlib.rcsetup.interactive_bk:
67 67 figManager = matplotlib._pylab_helpers.Gcf.get_active()
68 68 if figManager is not None:
69 69 canvas = figManager.canvas
70 70 if canvas.figure.stale:
71 71 canvas.draw()
72 72 try:
73 73 canvas.start_event_loop(interval)
74 74 except:
75 75 pass
76 76 return
77 77
78 78
79 79 def popup(message):
80 80 '''
81 81 '''
82 82
83 83 fig = plt.figure(figsize=(12, 8), facecolor='r')
84 84 text = '\n'.join([s.strip() for s in message.split(':')])
85 85 fig.text(0.01, 0.5, text, ha='left', va='center',
86 86 size='20', weight='heavy', color='w')
87 87 fig.show()
88 88 figpause(1000)
89 89
90 90
91 91 class Throttle(object):
92 92 '''
93 93 Decorator that prevents a function from being called more than once every
94 94 time period.
95 95 To create a function that cannot be called more than once a minute, but
96 96 will sleep until it can be called:
97 97 @Throttle(minutes=1)
98 98 def foo():
99 99 pass
100 100
101 101 for i in range(10):
102 102 foo()
103 103 print "This function has run %s times." % i
104 104 '''
105 105
106 106 def __init__(self, seconds=0, minutes=0, hours=0):
107 107 self.throttle_period = datetime.timedelta(
108 108 seconds=seconds, minutes=minutes, hours=hours
109 109 )
110 110
111 111 self.time_of_last_call = datetime.datetime.min
112 112
113 113 def __call__(self, fn):
114 114 @wraps(fn)
115 115 def wrapper(*args, **kwargs):
116 116 coerce = kwargs.pop('coerce', None)
117 117 if coerce:
118 118 self.time_of_last_call = datetime.datetime.now()
119 119 return fn(*args, **kwargs)
120 120 else:
121 121 now = datetime.datetime.now()
122 122 time_since_last_call = now - self.time_of_last_call
123 123 time_left = self.throttle_period - time_since_last_call
124 124
125 125 if time_left > datetime.timedelta(seconds=0):
126 126 return
127 127
128 128 self.time_of_last_call = datetime.datetime.now()
129 129 return fn(*args, **kwargs)
130 130
131 131 return wrapper
132 132
133 133 def apply_throttle(value):
134 134
135 135 @Throttle(seconds=value)
136 136 def fnThrottled(fn):
137 137 fn()
138 138
139 139 return fnThrottled
140 140
141 141
142 142 @MPDecorator
143 143 class Plot(Operation):
144 144 '''
145 145 Base class for Schain plotting operations
146 146 '''
147 147
148 148 CODE = 'Figure'
149 149 colormap = 'jet'
150 150 bgcolor = 'white'
151 151 __missing = 1E30
152 152
153 153 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
154 154 'zlimits', 'xlabel', 'ylabel', 'xaxis', 'cb_label', 'title',
155 155 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
156 156 'showprofile', 'decimation', 'pause']
157 157
158 158 def __init__(self):
159 159
160 160 Operation.__init__(self)
161 161 self.isConfig = False
162 162 self.isPlotConfig = False
163 163 self.save_counter = 1
164 164 self.sender_counter = 1
165 self.data = None
165 166
166 167 def __fmtTime(self, x, pos):
167 168 '''
168 169 '''
169 170
170 171 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
171 172
172 173 def __setup(self, **kwargs):
173 174 '''
174 175 Initialize variables
175 176 '''
176 177
177 178 self.figures = []
178 179 self.axes = []
179 180 self.cb_axes = []
180 181 self.localtime = kwargs.pop('localtime', True)
181 182 self.show = kwargs.get('show', True)
182 183 self.save = kwargs.get('save', False)
183 184 self.save_period = kwargs.get('save_period', 2)
184 185 self.ftp = kwargs.get('ftp', False)
185 186 self.colormap = kwargs.get('colormap', self.colormap)
186 187 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
187 188 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
188 189 self.colormaps = kwargs.get('colormaps', None)
189 190 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
190 191 self.showprofile = kwargs.get('showprofile', False)
191 192 self.title = kwargs.get('wintitle', self.CODE.upper())
192 193 self.cb_label = kwargs.get('cb_label', None)
193 194 self.cb_labels = kwargs.get('cb_labels', None)
194 195 self.labels = kwargs.get('labels', None)
195 196 self.xaxis = kwargs.get('xaxis', 'frequency')
196 197 self.zmin = kwargs.get('zmin', None)
197 198 self.zmax = kwargs.get('zmax', None)
198 199 self.zlimits = kwargs.get('zlimits', None)
199 200 self.xmin = kwargs.get('xmin', None)
200 201 self.xmax = kwargs.get('xmax', None)
201 202 self.xrange = kwargs.get('xrange', 24)
202 203 self.xscale = kwargs.get('xscale', None)
203 204 self.ymin = kwargs.get('ymin', None)
204 205 self.ymax = kwargs.get('ymax', None)
205 206 self.yscale = kwargs.get('yscale', None)
206 207 self.xlabel = kwargs.get('xlabel', None)
207 208 self.decimation = kwargs.get('decimation', None)
208 209 self.showSNR = kwargs.get('showSNR', False)
209 210 self.oneFigure = kwargs.get('oneFigure', True)
210 211 self.width = kwargs.get('width', None)
211 212 self.height = kwargs.get('height', None)
212 213 self.colorbar = kwargs.get('colorbar', True)
213 214 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
214 215 self.channels = kwargs.get('channels', None)
215 216 self.titles = kwargs.get('titles', [])
216 217 self.polar = False
217 218 self.type = kwargs.get('type', 'iq')
218 219 self.grid = kwargs.get('grid', False)
219 220 self.pause = kwargs.get('pause', False)
220 221 self.save_labels = kwargs.get('save_labels', None)
221 222 self.realtime = kwargs.get('realtime', True)
222 223 self.buffering = kwargs.get('buffering', True)
223 224 self.throttle = kwargs.get('throttle', 2)
224 225 self.exp_code = kwargs.get('exp_code', None)
225 226 self.plot_server = kwargs.get('plot_server', False)
226 227 self.sender_period = kwargs.get('sender_period', 2)
227 228 self.__throttle_plot = apply_throttle(self.throttle)
228 229 self.data = PlotterData(
229 self.CODE, self.throttle, self.exp_code, self.buffering)
230 self.CODE, self.throttle, self.exp_code, self.buffering, snr=self.showSNR)
230 231
231 232 if self.plot_server:
232 233 if not self.plot_server.startswith('tcp://'):
233 234 self.plot_server = 'tcp://{}'.format(self.plot_server)
234 235 log.success(
235 236 'Sending to server: {}'.format(self.plot_server),
236 237 self.name
237 238 )
239 if 'plot_name' in kwargs:
240 self.plot_name = kwargs['plot_name']
238 241
239 242 def __setup_plot(self):
240 243 '''
241 244 Common setup for all figures, here figures and axes are created
242 245 '''
243 246
244 247 self.setup()
245 248
246 self.time_label = 'LT' if self.localtime else 'UTC'
247 if self.data.localtime:
248 self.getDateTime = datetime.datetime.fromtimestamp
249 else:
250 self.getDateTime = datetime.datetime.utcfromtimestamp
249 self.time_label = 'LT' if self.localtime else 'UTC'
251 250
252 251 if self.width is None:
253 252 self.width = 8
254 253
255 254 self.figures = []
256 255 self.axes = []
257 256 self.cb_axes = []
258 257 self.pf_axes = []
259 258 self.cmaps = []
260 259
261 260 size = '15%' if self.ncols == 1 else '30%'
262 261 pad = '4%' if self.ncols == 1 else '8%'
263 262
264 263 if self.oneFigure:
265 264 if self.height is None:
266 265 self.height = 1.4 * self.nrows + 1
267 266 fig = plt.figure(figsize=(self.width, self.height),
268 267 edgecolor='k',
269 268 facecolor='w')
270 269 self.figures.append(fig)
271 270 for n in range(self.nplots):
272 271 ax = fig.add_subplot(self.nrows, self.ncols,
273 272 n + 1, polar=self.polar)
274 273 ax.tick_params(labelsize=8)
275 274 ax.firsttime = True
276 275 ax.index = 0
277 276 ax.press = None
278 277 self.axes.append(ax)
279 278 if self.showprofile:
280 279 cax = self.__add_axes(ax, size=size, pad=pad)
281 280 cax.tick_params(labelsize=8)
282 281 self.pf_axes.append(cax)
283 282 else:
284 283 if self.height is None:
285 284 self.height = 3
286 285 for n in range(self.nplots):
287 286 fig = plt.figure(figsize=(self.width, self.height),
288 287 edgecolor='k',
289 288 facecolor='w')
290 289 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
291 290 ax.tick_params(labelsize=8)
292 291 ax.firsttime = True
293 292 ax.index = 0
294 293 ax.press = None
295 294 self.figures.append(fig)
296 295 self.axes.append(ax)
297 296 if self.showprofile:
298 297 cax = self.__add_axes(ax, size=size, pad=pad)
299 298 cax.tick_params(labelsize=8)
300 299 self.pf_axes.append(cax)
301 300
302 301 for n in range(self.nrows):
303 302 if self.colormaps is not None:
304 303 cmap = plt.get_cmap(self.colormaps[n])
305 304 else:
306 305 cmap = plt.get_cmap(self.colormap)
307 306 cmap.set_bad(self.bgcolor, 1.)
308 307 self.cmaps.append(cmap)
309
308
310 309 for fig in self.figures:
311 310 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
312 311 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
313 312 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
314 313 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
315 314 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
316 if self.show:
317 fig.show()
318 315
319 316 def OnKeyPress(self, event):
320 317 '''
321 318 Event for pressing keys (up, down) change colormap
322 319 '''
323 320 ax = event.inaxes
324 321 if ax in self.axes:
325 322 if event.key == 'down':
326 323 ax.index += 1
327 324 elif event.key == 'up':
328 325 ax.index -= 1
329 326 if ax.index < 0:
330 327 ax.index = len(CMAPS) - 1
331 328 elif ax.index == len(CMAPS):
332 329 ax.index = 0
333 330 cmap = CMAPS[ax.index]
334 331 ax.cbar.set_cmap(cmap)
335 332 ax.cbar.draw_all()
336 333 ax.plt.set_cmap(cmap)
337 334 ax.cbar.patch.figure.canvas.draw()
338 335 self.colormap = cmap.name
339 336
340 337 def OnBtnScroll(self, event):
341 338 '''
342 339 Event for scrolling, scale figure
343 340 '''
344 341 cb_ax = event.inaxes
345 342 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
346 343 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
347 344 pt = ax.cbar.ax.bbox.get_points()[:, 1]
348 345 nrm = ax.cbar.norm
349 346 vmin, vmax, p0, p1, pS = (
350 347 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
351 348 scale = 2 if event.step == 1 else 0.5
352 349 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
353 350 ax.cbar.norm.vmin = point - scale * (point - vmin)
354 351 ax.cbar.norm.vmax = point - scale * (point - vmax)
355 352 ax.plt.set_norm(ax.cbar.norm)
356 353 ax.cbar.draw_all()
357 354 ax.cbar.patch.figure.canvas.draw()
358 355
359 356 def onBtnPress(self, event):
360 357 '''
361 358 Event for mouse button press
362 359 '''
363 360 cb_ax = event.inaxes
364 361 if cb_ax is None:
365 362 return
366 363
367 364 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
368 365 cb_ax.press = event.x, event.y
369 366 else:
370 367 cb_ax.press = None
371 368
372 369 def onMotion(self, event):
373 370 '''
374 371 Event for move inside colorbar
375 372 '''
376 373 cb_ax = event.inaxes
377 374 if cb_ax is None:
378 375 return
379 376 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
380 377 return
381 378 if cb_ax.press is None:
382 379 return
383 380
384 381 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
385 382 xprev, yprev = cb_ax.press
386 383 dx = event.x - xprev
387 384 dy = event.y - yprev
388 385 cb_ax.press = event.x, event.y
389 386 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
390 387 perc = 0.03
391 388
392 389 if event.button == 1:
393 390 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
394 391 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
395 392 elif event.button == 3:
396 393 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
397 394 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
398 395
399 396 ax.cbar.draw_all()
400 397 ax.plt.set_norm(ax.cbar.norm)
401 398 ax.cbar.patch.figure.canvas.draw()
402 399
403 400 def onBtnRelease(self, event):
404 401 '''
405 402 Event for mouse button release
406 403 '''
407 404 cb_ax = event.inaxes
408 405 if cb_ax is not None:
409 406 cb_ax.press = None
410 407
411 408 def __add_axes(self, ax, size='30%', pad='8%'):
412 409 '''
413 410 Add new axes to the given figure
414 411 '''
415 412 divider = make_axes_locatable(ax)
416 413 nax = divider.new_horizontal(size=size, pad=pad)
417 414 ax.figure.add_axes(nax)
418 415 return nax
419 416
420 417 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
421 418 '''
422 419 Create a masked array for missing data
423 420 '''
424 421 if x_buffer.shape[0] < 2:
425 422 return x_buffer, y_buffer, z_buffer
426 423
427 424 deltas = x_buffer[1:] - x_buffer[0:-1]
428 425 x_median = numpy.median(deltas)
429 426
430 427 index = numpy.where(deltas > 5 * x_median)
431 428
432 429 if len(index[0]) != 0:
433 430 z_buffer[::, index[0], ::] = self.__missing
434 431 z_buffer = numpy.ma.masked_inside(z_buffer,
435 432 0.99 * self.__missing,
436 433 1.01 * self.__missing)
437 434
438 435 return x_buffer, y_buffer, z_buffer
439 436
440 437 def decimate(self):
441 438
442 439 # dx = int(len(self.x)/self.__MAXNUMX) + 1
443 440 dy = int(len(self.y) / self.decimation) + 1
444 441
445 442 # x = self.x[::dx]
446 443 x = self.x
447 444 y = self.y[::dy]
448 445 z = self.z[::, ::, ::dy]
449 446
450 447 return x, y, z
451 448
452 449 def format(self):
453 450 '''
454 451 Set min and max values, labels, ticks and titles
455 452 '''
456 453
457 454 if self.xmin is None:
458 455 xmin = self.data.min_time
459 456 else:
460 457 if self.xaxis is 'time':
461 458 dt = self.getDateTime(self.data.min_time)
462 459 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
463 460 datetime.datetime(1970, 1, 1)).total_seconds()
464 461 if self.data.localtime:
465 462 xmin += time.timezone
466 self.tmin = xmin
467 463 else:
468 464 xmin = self.xmin
469 465
470 466 if self.xmax is None:
471 467 xmax = xmin + self.xrange * 60 * 60
472 468 else:
473 469 if self.xaxis is 'time':
474 470 dt = self.getDateTime(self.data.max_time)
475 471 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
476 472 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
477 473 if self.data.localtime:
478 474 xmax += time.timezone
479 475 else:
480 476 xmax = self.xmax
481 477
482 478 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
483 479 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
484 480 #Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000])
485 481
486 482 #i = 1 if numpy.where(
487 483 # abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
488 484 #ystep = Y[i] / 10.
489 485 dig = int(numpy.log10(ymax))
490 486 if dig == 0:
491 487 digD = len(str(ymax)) - 2
492 488 ydec = ymax*(10**digD)
493 489
494 490 dig = int(numpy.log10(ydec))
495 491 ystep = ((ydec + (10**(dig)))//10**(dig))*(10**(dig))
496 492 ystep = ystep/5
497 493 ystep = ystep/(10**digD)
498 494
499 495 else:
500 496 ystep = ((ymax + (10**(dig)))//10**(dig))*(10**(dig))
501 497 ystep = ystep/5
502 498
503 499 if self.xaxis is not 'time':
504 500
505 501 dig = int(numpy.log10(xmax))
506 502
507 503 if dig <= 0:
508 504 digD = len(str(xmax)) - 2
509 505 xdec = xmax*(10**digD)
510 506
511 507 dig = int(numpy.log10(xdec))
512 508 xstep = ((xdec + (10**(dig)))//10**(dig))*(10**(dig))
513 509 xstep = xstep*0.5
514 510 xstep = xstep/(10**digD)
515 511
516 512 else:
517 513 xstep = ((xmax + (10**(dig)))//10**(dig))*(10**(dig))
518 514 xstep = xstep/5
519 515
520 516 for n, ax in enumerate(self.axes):
521 517 if ax.firsttime:
522 518 ax.set_facecolor(self.bgcolor)
523 519 ax.yaxis.set_major_locator(MultipleLocator(ystep))
524 520 if self.xscale:
525 521 ax.xaxis.set_major_formatter(FuncFormatter(
526 522 lambda x, pos: '{0:g}'.format(x*self.xscale)))
527 523 if self.xscale:
528 524 ax.yaxis.set_major_formatter(FuncFormatter(
529 525 lambda x, pos: '{0:g}'.format(x*self.yscale)))
530 526 if self.xaxis is 'time':
531 527 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
532 528 ax.xaxis.set_major_locator(LinearLocator(9))
533 529 else:
534 530 ax.xaxis.set_major_locator(MultipleLocator(xstep))
535 531 if self.xlabel is not None:
536 532 ax.set_xlabel(self.xlabel)
537 533 ax.set_ylabel(self.ylabel)
538 534 ax.firsttime = False
539 535 if self.showprofile:
540 536 self.pf_axes[n].set_ylim(ymin, ymax)
541 537 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
542 538 self.pf_axes[n].set_xlabel('dB')
543 539 self.pf_axes[n].grid(b=True, axis='x')
544 540 [tick.set_visible(False)
545 541 for tick in self.pf_axes[n].get_yticklabels()]
546 542 if self.colorbar:
547 543 ax.cbar = plt.colorbar(
548 544 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
549 545 ax.cbar.ax.tick_params(labelsize=8)
550 546 ax.cbar.ax.press = None
551 547 if self.cb_label:
552 548 ax.cbar.set_label(self.cb_label, size=8)
553 549 elif self.cb_labels:
554 550 ax.cbar.set_label(self.cb_labels[n], size=8)
555 551 else:
556 552 ax.cbar = None
557 553 if self.grid:
558 554 ax.grid(True)
559 555
560 556 if not self.polar:
561 557 ax.set_xlim(xmin, xmax)
562 558 ax.set_ylim(ymin, ymax)
563 559 ax.set_title('{} {} {}'.format(
564 560 self.titles[n],
565 561 self.getDateTime(self.data.max_time).strftime(
566 '%H:%M:%S'),
562 '%Y-%m-%d %H:%M:%S'),
567 563 self.time_label),
568 564 size=8)
569 565 else:
570 566 ax.set_title('{}'.format(self.titles[n]), size=8)
571 567 ax.set_ylim(0, 90)
572 568 ax.set_yticks(numpy.arange(0, 90, 20))
573 569 ax.yaxis.labelpad = 40
574 570
575 571 def clear_figures(self):
576 572 '''
577 573 Reset axes for redraw plots
578 574 '''
579 575
580 576 for ax in self.axes:
581 577 ax.clear()
582 578 ax.firsttime = True
583 579 if ax.cbar:
584 580 ax.cbar.remove()
585 581
586 582 def __plot(self):
587 583 '''
588 584 Main function to plot, format and save figures
589 585 '''
590 586
591 587 try:
592 588 self.plot()
593 589 self.format()
594 590 except Exception as e:
595 591 log.warning('{} Plot could not be updated... check data'.format(
596 592 self.CODE), self.name)
597 593 log.error(str(e), '')
598 594 return
599 595
600 596 for n, fig in enumerate(self.figures):
601 597 if self.nrows == 0 or self.nplots == 0:
602 598 log.warning('No data', self.name)
603 599 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
604 600 fig.canvas.manager.set_window_title(self.CODE)
605 601 continue
606 602
607 603 fig.tight_layout()
608 604 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
609 605 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
610 606 fig.canvas.draw()
607 if self.show:
608 fig.show()
609 figpause(0.1)
611 610
612 611 if self.save:
613 612 self.save_figure(n)
614 613
615 614 if self.plot_server:
616 615 self.send_to_server()
617 616 # t = Thread(target=self.send_to_server)
618 617 # t.start()
619 618
620 619 def save_figure(self, n):
621 620 '''
622 621 '''
623 622
624 623 if self.save_counter < self.save_period:
625 624 self.save_counter += 1
626 625 return
627 626
628 627 self.save_counter = 1
629 628
630 629 fig = self.figures[n]
631 630
632 631 if self.save_labels:
633 632 labels = self.save_labels
634 633 else:
635 634 labels = list(range(self.nrows))
636 635
637 636 if self.oneFigure:
638 637 label = ''
639 638 else:
640 639 label = '-{}'.format(labels[n])
641 640 figname = os.path.join(
642 641 self.save,
643 642 self.CODE,
644 643 '{}{}_{}.png'.format(
645 644 self.CODE,
646 645 label,
647 646 self.getDateTime(self.data.max_time).strftime(
648 647 '%Y%m%d_%H%M%S'
649 648 ),
650 649 )
651 650 )
652 651 log.log('Saving figure: {}'.format(figname), self.name)
653 652 if not os.path.isdir(os.path.dirname(figname)):
654 653 os.makedirs(os.path.dirname(figname))
655 654 fig.savefig(figname)
656 655
657 656 if self.realtime:
658 657 figname = os.path.join(
659 658 self.save,
660 659 '{}{}_{}.png'.format(
661 660 self.CODE,
662 661 label,
663 662 self.getDateTime(self.data.min_time).strftime(
664 663 '%Y%m%d'
665 664 ),
666 665 )
667 666 )
668 667 fig.savefig(figname)
669 668
670 669 def send_to_server(self):
671 670 '''
672 671 '''
673 672
674 673 if self.sender_counter < self.sender_period:
675 674 self.sender_counter += 1
676 675
677 676 self.sender_counter = 1
678
677 self.data.meta['titles'] = self.titles
679 678 retries = 2
680 679 while True:
681 self.socket.send_string(self.data.jsonify())
680 self.socket.send_string(self.data.jsonify(self.plot_name, self.plot_type))
682 681 socks = dict(self.poll.poll(5000))
683 682 if socks.get(self.socket) == zmq.POLLIN:
684 683 reply = self.socket.recv_string()
685 684 if reply == 'ok':
686 685 log.log("Response from server ok", self.name)
687 686 break
688 687 else:
689 688 log.warning(
690 689 "Malformed reply from server: {}".format(reply), self.name)
691 690
692 691 else:
693 692 log.warning(
694 693 "No response from server, retrying...", self.name)
695 694 self.socket.setsockopt(zmq.LINGER, 0)
696 695 self.socket.close()
697 696 self.poll.unregister(self.socket)
698 697 retries -= 1
699 698 if retries == 0:
700 699 log.error(
701 700 "Server seems to be offline, abandoning", self.name)
702 701 self.socket = self.context.socket(zmq.REQ)
703 702 self.socket.connect(self.plot_server)
704 703 self.poll.register(self.socket, zmq.POLLIN)
705 704 time.sleep(1)
706 705 break
707 706 self.socket = self.context.socket(zmq.REQ)
708 707 self.socket.connect(self.plot_server)
709 708 self.poll.register(self.socket, zmq.POLLIN)
710 709 time.sleep(0.5)
711 710
712 711 def setup(self):
713 712 '''
714 713 This method should be implemented in the child class, the following
715 714 attributes should be set:
716 715
717 716 self.nrows: number of rows
718 717 self.ncols: number of cols
719 718 self.nplots: number of plots (channels or pairs)
720 719 self.ylabel: label for Y axes
721 720 self.titles: list of axes title
722 721
723 722 '''
724 723 raise NotImplementedError
725 724
726 725 def plot(self):
727 726 '''
728 727 Must be defined in the child class
729 728 '''
730 729 raise NotImplementedError
731 730
732 731 def run(self, dataOut, **kwargs):
733
734 if dataOut.error:
735 coerce = True
736 else:
737 coerce = False
732 '''
733 Main plotting routine
734 '''
738 735
739 736 if self.isConfig is False:
740 737 self.__setup(**kwargs)
738 if dataOut.type == 'Parameters':
739 t = dataOut.utctimeInit
740 else:
741 t = dataOut.utctime
742
743 if dataOut.useLocalTime:
744 self.getDateTime = datetime.datetime.fromtimestamp
745 if not self.localtime:
746 t += time.timezone
747 else:
748 self.getDateTime = datetime.datetime.utcfromtimestamp
749 if self.localtime:
750 t -= time.timezone
751
752 if self.xmin is None:
753 self.tmin = t
754 else:
755 self.tmin = (
756 self.getDateTime(t).replace(
757 hour=self.xmin,
758 minute=0,
759 second=0) - self.getDateTime(0)).total_seconds()
760
741 761 self.data.setup()
742 762 self.isConfig = True
743 763 if self.plot_server:
744 764 self.context = zmq.Context()
745 765 self.socket = self.context.socket(zmq.REQ)
746 766 self.socket.connect(self.plot_server)
747 767 self.poll = zmq.Poller()
748 768 self.poll.register(self.socket, zmq.POLLIN)
749 769
750 770 if dataOut.type == 'Parameters':
751 771 tm = dataOut.utctimeInit
752 772 else:
753 773 tm = dataOut.utctime
754
755 if dataOut.useLocalTime:
756 if not self.localtime:
757 tm += time.timezone
758 else:
759 if self.localtime:
760 tm -= time.timezone
761 774
762 if self.xaxis is 'time' and self.data and (tm - self.tmin) >= self.xrange*60*60:
775 if not dataOut.useLocalTime and self.localtime:
776 tm -= time.timezone
777 if dataOut.useLocalTime and not self.localtime:
778 tm += time.timezone
779
780 if self.xaxis is 'time' and self.data and (tm - self.tmin) >= self.xrange*60*60:
781 self.save_counter = self.save_period
763 782 self.__plot()
783 self.xmin += self.xrange
784 if self.xmin >= 24:
785 self.xmin -= 24
786 self.tmin += self.xrange*60*60
764 787 self.data.setup()
765 788 self.clear_figures()
766 789
767 790 self.data.update(dataOut, tm)
768 791
769 792 if self.isPlotConfig is False:
770 793 self.__setup_plot()
771 794 self.isPlotConfig = True
772 795
773 796 if self.realtime:
774 797 self.__plot()
775 798 else:
776 self.__throttle_plot(self.__plot, coerce=coerce)
777
778 figpause(0.001)
799 self.__throttle_plot(self.__plot)#, coerce=coerce)
779 800
780 801 def close(self):
781 802
803 if self.data:
804 self.save_counter = self.save_period
805 self.__plot()
782 806 if self.data and self.pause:
783 807 figpause(10)
784 808
@@ -1,748 +1,759
1 1 '''
2 2 New Plots Operations
3 3
4 4 @author: juan.espinoza@jro.igp.gob.pe
5 5 '''
6 6
7 7
8 8 import time
9 9 import datetime
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt
13 13 from schainpy.utils import log
14 14
15 15 EARTH_RADIUS = 6.3710e3
16 16
17 17
18 18 def ll2xy(lat1, lon1, lat2, lon2):
19 19
20 20 p = 0.017453292519943295
21 21 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
22 22 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
23 23 r = 12742 * numpy.arcsin(numpy.sqrt(a))
24 24 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
25 25 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
26 26 theta = -theta + numpy.pi/2
27 27 return r*numpy.cos(theta), r*numpy.sin(theta)
28 28
29 29
30 30 def km2deg(km):
31 31 '''
32 32 Convert distance in km to degrees
33 33 '''
34 34
35 35 return numpy.rad2deg(km/EARTH_RADIUS)
36 36
37 37
38 38 class SpectraPlot(Plot):
39 39 '''
40 40 Plot for Spectra data
41 41 '''
42 42
43 43 CODE = 'spc'
44 44 colormap = 'jro'
45 plot_name = 'Spectra'
46 plot_type = 'pcolor'
45 47
46 48 def setup(self):
47 49 self.nplots = len(self.data.channels)
48 50 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
49 51 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
50 52 self.width = 3.4 * self.ncols
51 53 self.height = 3 * self.nrows
52 54 self.cb_label = 'dB'
53 55 if self.showprofile:
54 56 self.width += 0.8 * self.ncols
55 57
56 58 self.ylabel = 'Range [km]'
57 59
58 60 def plot(self):
59 61 if self.xaxis == "frequency":
60 62 x = self.data.xrange[0]
61 63 self.xlabel = "Frequency (kHz)"
62 64 elif self.xaxis == "time":
63 65 x = self.data.xrange[1]
64 66 self.xlabel = "Time (ms)"
65 67 else:
66 68 x = self.data.xrange[2]
67 69 self.xlabel = "Velocity (m/s)"
68 70
69 71 if self.CODE == 'spc_moments':
70 72 x = self.data.xrange[2]
71 73 self.xlabel = "Velocity (m/s)"
72 74
73 75 self.titles = []
74 76
75 77 y = self.data.heights
76 78 self.y = y
77 79 z = self.data['spc']
78 80
79 81 for n, ax in enumerate(self.axes):
80 82 noise = self.data['noise'][n][-1]
81 83 if self.CODE == 'spc_moments':
82 84 mean = self.data['moments'][n, :, 1, :][-1]
83 85 if ax.firsttime:
84 86 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
85 87 self.xmin = self.xmin if self.xmin else -self.xmax
86 88 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
87 89 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
88 90 ax.plt = ax.pcolormesh(x, y, z[n].T,
89 91 vmin=self.zmin,
90 92 vmax=self.zmax,
91 93 cmap=plt.get_cmap(self.colormap)
92 94 )
93 95
94 96 if self.showprofile:
95 97 ax.plt_profile = self.pf_axes[n].plot(
96 98 self.data['rti'][n][-1], y)[0]
97 99 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
98 100 color="k", linestyle="dashed", lw=1)[0]
99 101 if self.CODE == 'spc_moments':
100 102 ax.plt_mean = ax.plot(mean, y, color='k')[0]
101 103 else:
102 104 ax.plt.set_array(z[n].T.ravel())
103 105 if self.showprofile:
104 106 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
105 107 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
106 108 if self.CODE == 'spc_moments':
107 109 ax.plt_mean.set_data(mean, y)
108 110 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
109 111
110 112
111 113 class CrossSpectraPlot(Plot):
112 114
113 115 CODE = 'cspc'
114 116 colormap = 'jet'
117 plot_name = 'CrossSpectra'
118 plot_type = 'pcolor'
115 119 zmin_coh = None
116 120 zmax_coh = None
117 121 zmin_phase = None
118 122 zmax_phase = None
119 123
120 124 def setup(self):
121 125
122 126 self.ncols = 4
123 127 self.nrows = len(self.data.pairs)
124 128 self.nplots = self.nrows * 4
125 129 self.width = 3.4 * self.ncols
126 130 self.height = 3 * self.nrows
127 131 self.ylabel = 'Range [km]'
128 132 self.showprofile = False
129 133
130 134 def plot(self):
131 135
132 136 if self.xaxis == "frequency":
133 137 x = self.data.xrange[0]
134 138 self.xlabel = "Frequency (kHz)"
135 139 elif self.xaxis == "time":
136 140 x = self.data.xrange[1]
137 141 self.xlabel = "Time (ms)"
138 142 else:
139 143 x = self.data.xrange[2]
140 144 self.xlabel = "Velocity (m/s)"
141 145
142 146 self.titles = []
143 147
144 148 y = self.data.heights
145 149 self.y = y
146 150 spc = self.data['spc']
147 151 cspc = self.data['cspc']
148 152
149 153 for n in range(self.nrows):
150 154 noise = self.data['noise'][n][-1]
151 155 pair = self.data.pairs[n]
152 156 ax = self.axes[4 * n]
153 157 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
154 158 if ax.firsttime:
155 159 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
156 160 self.xmin = self.xmin if self.xmin else -self.xmax
157 161 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
158 162 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
159 163 ax.plt = ax.pcolormesh(x , y , spc0.T,
160 164 vmin=self.zmin,
161 165 vmax=self.zmax,
162 166 cmap=plt.get_cmap(self.colormap)
163 167 )
164 168 else:
165 169 ax.plt.set_array(spc0.T.ravel())
166 170 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
167 171
168 172 ax = self.axes[4 * n + 1]
169 173 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
170 174 if ax.firsttime:
171 175 ax.plt = ax.pcolormesh(x , y, spc1.T,
172 176 vmin=self.zmin,
173 177 vmax=self.zmax,
174 178 cmap=plt.get_cmap(self.colormap)
175 179 )
176 180 else:
177 181 ax.plt.set_array(spc1.T.ravel())
178 182 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
179 183
180 184 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
181 185 coh = numpy.abs(out)
182 186 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
183 187
184 188 ax = self.axes[4 * n + 2]
185 189 if ax.firsttime:
186 190 ax.plt = ax.pcolormesh(x, y, coh.T,
187 191 vmin=0,
188 192 vmax=1,
189 193 cmap=plt.get_cmap(self.colormap_coh)
190 194 )
191 195 else:
192 196 ax.plt.set_array(coh.T.ravel())
193 197 self.titles.append(
194 198 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
195 199
196 200 ax = self.axes[4 * n + 3]
197 201 if ax.firsttime:
198 202 ax.plt = ax.pcolormesh(x, y, phase.T,
199 203 vmin=-180,
200 204 vmax=180,
201 205 cmap=plt.get_cmap(self.colormap_phase)
202 206 )
203 207 else:
204 208 ax.plt.set_array(phase.T.ravel())
205 209 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
206 210
207 211
208 212 class SpectralMomentsPlot(SpectraPlot):
209 213 '''
210 214 Plot for Spectral Moments
211 215 '''
212 216 CODE = 'spc_moments'
213 217 colormap = 'jro'
218 plot_name = 'SpectralMoments'
219 plot_type = 'pcolor'
214 220
215 221
216 222 class RTIPlot(Plot):
217 223 '''
218 224 Plot for RTI data
219 225 '''
220 226
221 227 CODE = 'rti'
222 228 colormap = 'jro'
229 plot_name = 'RTI'
230 plot_type = 'pcolorbuffer'
223 231
224 232 def setup(self):
225 233 self.xaxis = 'time'
226 234 self.ncols = 1
227 235 self.nrows = len(self.data.channels)
228 236 self.nplots = len(self.data.channels)
229 237 self.ylabel = 'Range [km]'
230 238 self.cb_label = 'dB'
231 239 self.titles = ['{} Channel {}'.format(
232 240 self.CODE.upper(), x) for x in range(self.nrows)]
233 241
234 242 def plot(self):
235 243 self.x = self.data.times
236 244 self.y = self.data.heights
237 245 self.z = self.data[self.CODE]
238 246 self.z = numpy.ma.masked_invalid(self.z)
239 247
240 248 if self.decimation is None:
241 249 x, y, z = self.fill_gaps(self.x, self.y, self.z)
242 250 else:
243 251 x, y, z = self.fill_gaps(*self.decimate())
244 252
245 253 for n, ax in enumerate(self.axes):
246 254 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
247 255 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
248 256 if ax.firsttime:
249 257 ax.plt = ax.pcolormesh(x, y, z[n].T,
250 258 vmin=self.zmin,
251 259 vmax=self.zmax,
252 260 cmap=plt.get_cmap(self.colormap)
253 261 )
254 262 if self.showprofile:
255 263 ax.plot_profile = self.pf_axes[n].plot(
256 264 self.data['rti'][n][-1], self.y)[0]
257 265 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
258 266 color="k", linestyle="dashed", lw=1)[0]
259 267 else:
260 268 ax.collections.remove(ax.collections[0])
261 269 ax.plt = ax.pcolormesh(x, y, z[n].T,
262 270 vmin=self.zmin,
263 271 vmax=self.zmax,
264 272 cmap=plt.get_cmap(self.colormap)
265 273 )
266 274 if self.showprofile:
267 275 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
268 276 ax.plot_noise.set_data(numpy.repeat(
269 277 self.data['noise'][n][-1], len(self.y)), self.y)
270 278
271 279
272 280 class CoherencePlot(RTIPlot):
273 281 '''
274 282 Plot for Coherence data
275 283 '''
276 284
277 285 CODE = 'coh'
286 plot_name = 'Coherence'
278 287
279 288 def setup(self):
280 289 self.xaxis = 'time'
281 290 self.ncols = 1
282 291 self.nrows = len(self.data.pairs)
283 292 self.nplots = len(self.data.pairs)
284 293 self.ylabel = 'Range [km]'
285 294 if self.CODE == 'coh':
286 295 self.cb_label = ''
287 296 self.titles = [
288 297 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
289 298 else:
290 299 self.cb_label = 'Degrees'
291 300 self.titles = [
292 301 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
293 302
294 303
295 304 class PhasePlot(CoherencePlot):
296 305 '''
297 306 Plot for Phase map data
298 307 '''
299 308
300 309 CODE = 'phase'
301 310 colormap = 'seismic'
311 plot_name = 'Phase'
302 312
303 313
304 314 class NoisePlot(Plot):
305 315 '''
306 316 Plot for noise
307 317 '''
308 318
309 319 CODE = 'noise'
320 plot_name = 'Noise'
321 plot_type = 'scatterbuffer'
322
310 323
311 324 def setup(self):
312 325 self.xaxis = 'time'
313 326 self.ncols = 1
314 327 self.nrows = 1
315 328 self.nplots = 1
316 329 self.ylabel = 'Intensity [dB]'
317 330 self.titles = ['Noise']
318 331 self.colorbar = False
319 332
320 333 def plot(self):
321 334
322 335 x = self.data.times
323 336 xmin = self.data.min_time
324 337 xmax = xmin + self.xrange * 60 * 60
325 338 Y = self.data[self.CODE]
326 339
327 340 if self.axes[0].firsttime:
328 341 for ch in self.data.channels:
329 342 y = Y[ch]
330 343 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
331 344 plt.legend()
332 345 else:
333 346 for ch in self.data.channels:
334 347 y = Y[ch]
335 348 self.axes[0].lines[ch].set_data(x, y)
336 349
337 350 self.ymin = numpy.nanmin(Y) - 5
338 351 self.ymax = numpy.nanmax(Y) + 5
339 352
340 353
341 354 class SnrPlot(RTIPlot):
342 355 '''
343 356 Plot for SNR Data
344 357 '''
345 358
346 359 CODE = 'snr'
347 360 colormap = 'jet'
361 plot_name = 'SNR'
348 362
349 363
350 364 class DopplerPlot(RTIPlot):
351 365 '''
352 366 Plot for DOPPLER Data
353 367 '''
354 368
355 369 CODE = 'dop'
356 370 colormap = 'jet'
371 plot_name = 'Doppler'
357 372
358 373
359 374 class SkyMapPlot(Plot):
360 375 '''
361 376 Plot for meteors detection data
362 377 '''
363 378
364 379 CODE = 'param'
365 380
366 381 def setup(self):
367 382
368 383 self.ncols = 1
369 384 self.nrows = 1
370 385 self.width = 7.2
371 386 self.height = 7.2
372 387 self.nplots = 1
373 388 self.xlabel = 'Zonal Zenith Angle (deg)'
374 389 self.ylabel = 'Meridional Zenith Angle (deg)'
375 390 self.polar = True
376 391 self.ymin = -180
377 392 self.ymax = 180
378 393 self.colorbar = False
379 394
380 395 def plot(self):
381 396
382 397 arrayParameters = numpy.concatenate(self.data['param'])
383 398 error = arrayParameters[:, -1]
384 399 indValid = numpy.where(error == 0)[0]
385 400 finalMeteor = arrayParameters[indValid, :]
386 401 finalAzimuth = finalMeteor[:, 3]
387 402 finalZenith = finalMeteor[:, 4]
388 403
389 404 x = finalAzimuth * numpy.pi / 180
390 405 y = finalZenith
391 406
392 407 ax = self.axes[0]
393 408
394 409 if ax.firsttime:
395 410 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
396 411 else:
397 412 ax.plot.set_data(x, y)
398 413
399 414 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
400 415 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
401 416 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
402 417 dt2,
403 418 len(x))
404 419 self.titles[0] = title
405 420
406 421
407 422 class ParametersPlot(RTIPlot):
408 423 '''
409 424 Plot for data_param object
410 425 '''
411 426
412 427 CODE = 'param'
413 428 colormap = 'seismic'
429 plot_name = 'Parameters'
414 430
415 431 def setup(self):
416 432 self.xaxis = 'time'
417 433 self.ncols = 1
418 434 self.nrows = self.data.shape(self.CODE)[0]
419 435 self.nplots = self.nrows
420 436 if self.showSNR:
421 437 self.nrows += 1
422 438 self.nplots += 1
423 439
424 440 self.ylabel = 'Height [km]'
425 441 if not self.titles:
426 442 self.titles = self.data.parameters \
427 443 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
428 444 if self.showSNR:
429 445 self.titles.append('SNR')
430 446
431 447 def plot(self):
432 448 self.data.normalize_heights()
433 449 self.x = self.data.times
434 450 self.y = self.data.heights
435 451 if self.showSNR:
436 452 self.z = numpy.concatenate(
437 453 (self.data[self.CODE], self.data['snr'])
438 454 )
439 455 else:
440 456 self.z = self.data[self.CODE]
441 457
442 458 self.z = numpy.ma.masked_invalid(self.z)
443 459
444 460 if self.decimation is None:
445 461 x, y, z = self.fill_gaps(self.x, self.y, self.z)
446 462 else:
447 463 x, y, z = self.fill_gaps(*self.decimate())
448 464
449 465 for n, ax in enumerate(self.axes):
450 466
451 467 self.zmax = self.zmax if self.zmax is not None else numpy.max(
452 468 self.z[n])
453 469 self.zmin = self.zmin if self.zmin is not None else numpy.min(
454 470 self.z[n])
455 471
456 472 if ax.firsttime:
457 473 if self.zlimits is not None:
458 474 self.zmin, self.zmax = self.zlimits[n]
459 475
460 476 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
461 477 vmin=self.zmin,
462 478 vmax=self.zmax,
463 479 cmap=self.cmaps[n]
464 480 )
465 481 else:
466 482 if self.zlimits is not None:
467 483 self.zmin, self.zmax = self.zlimits[n]
468 484 ax.collections.remove(ax.collections[0])
469 485 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
470 486 vmin=self.zmin,
471 487 vmax=self.zmax,
472 488 cmap=self.cmaps[n]
473 489 )
474 490
475 491
476 492 class OutputPlot(ParametersPlot):
477 493 '''
478 494 Plot data_output object
479 495 '''
480 496
481 497 CODE = 'output'
482 498 colormap = 'seismic'
499 plot_name = 'Output'
483 500
484 501
485 502 class PolarMapPlot(Plot):
486 503 '''
487 504 Plot for weather radar
488 505 '''
489 506
490 507 CODE = 'param'
491 508 colormap = 'seismic'
492 509
493 510 def setup(self):
494 511 self.ncols = 1
495 512 self.nrows = 1
496 513 self.width = 9
497 514 self.height = 8
498 515 self.mode = self.data.meta['mode']
499 516 if self.channels is not None:
500 517 self.nplots = len(self.channels)
501 518 self.nrows = len(self.channels)
502 519 else:
503 520 self.nplots = self.data.shape(self.CODE)[0]
504 521 self.nrows = self.nplots
505 522 self.channels = list(range(self.nplots))
506 523 if self.mode == 'E':
507 524 self.xlabel = 'Longitude'
508 525 self.ylabel = 'Latitude'
509 526 else:
510 527 self.xlabel = 'Range (km)'
511 528 self.ylabel = 'Height (km)'
512 529 self.bgcolor = 'white'
513 530 self.cb_labels = self.data.meta['units']
514 531 self.lat = self.data.meta['latitude']
515 532 self.lon = self.data.meta['longitude']
516 533 self.xmin, self.xmax = float(
517 534 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
518 535 self.ymin, self.ymax = float(
519 536 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
520 537 # self.polar = True
521 538
522 539 def plot(self):
523 540
524 541 for n, ax in enumerate(self.axes):
525 542 data = self.data['param'][self.channels[n]]
526 543
527 544 zeniths = numpy.linspace(
528 545 0, self.data.meta['max_range'], data.shape[1])
529 546 if self.mode == 'E':
530 547 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
531 548 r, theta = numpy.meshgrid(zeniths, azimuths)
532 549 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
533 550 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
534 551 x = km2deg(x) + self.lon
535 552 y = km2deg(y) + self.lat
536 553 else:
537 554 azimuths = numpy.radians(self.data.heights)
538 555 r, theta = numpy.meshgrid(zeniths, azimuths)
539 556 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
540 557 self.y = zeniths
541 558
542 559 if ax.firsttime:
543 560 if self.zlimits is not None:
544 561 self.zmin, self.zmax = self.zlimits[n]
545 562 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
546 563 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
547 564 vmin=self.zmin,
548 565 vmax=self.zmax,
549 566 cmap=self.cmaps[n])
550 567 else:
551 568 if self.zlimits is not None:
552 569 self.zmin, self.zmax = self.zlimits[n]
553 570 ax.collections.remove(ax.collections[0])
554 571 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
555 572 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
556 573 vmin=self.zmin,
557 574 vmax=self.zmax,
558 575 cmap=self.cmaps[n])
559 576
560 577 if self.mode == 'A':
561 578 continue
562 579
563 580 # plot district names
564 581 f = open('/data/workspace/schain_scripts/distrito.csv')
565 582 for line in f:
566 583 label, lon, lat = [s.strip() for s in line.split(',') if s]
567 584 lat = float(lat)
568 585 lon = float(lon)
569 586 # ax.plot(lon, lat, '.b', ms=2)
570 587 ax.text(lon, lat, label.decode('utf8'), ha='center',
571 588 va='bottom', size='8', color='black')
572 589
573 590 # plot limites
574 591 limites = []
575 592 tmp = []
576 593 for line in open('/data/workspace/schain_scripts/lima.csv'):
577 594 if '#' in line:
578 595 if tmp:
579 596 limites.append(tmp)
580 597 tmp = []
581 598 continue
582 599 values = line.strip().split(',')
583 600 tmp.append((float(values[0]), float(values[1])))
584 601 for points in limites:
585 602 ax.add_patch(
586 603 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
587 604
588 605 # plot Cuencas
589 606 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
590 607 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
591 608 values = [line.strip().split(',') for line in f]
592 609 points = [(float(s[0]), float(s[1])) for s in values]
593 610 ax.add_patch(Polygon(points, ec='b', fc='none'))
594 611
595 612 # plot grid
596 613 for r in (15, 30, 45, 60):
597 614 ax.add_artist(plt.Circle((self.lon, self.lat),
598 615 km2deg(r), color='0.6', fill=False, lw=0.2))
599 616 ax.text(
600 617 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
601 618 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
602 619 '{}km'.format(r),
603 620 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
604 621
605 622 if self.mode == 'E':
606 623 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
607 624 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
608 625 else:
609 626 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
610 627 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
611 628
612 629 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
613 630 self.titles = ['{} {}'.format(
614 631 self.data.parameters[x], title) for x in self.channels]
615 632
616 633
617 634 class ScopePlot(Plot):
618 635
619 636 '''
620 637 Plot for Scope
621 638 '''
622 639
623 640 CODE = 'scope'
641 plot_name = 'Scope'
642 plot_type = 'scatter'
624 643
625 644 def setup(self):
626 645
627 646 self.xaxis = 'Range (Km)'
628 647 self.ncols = 1
629 648 self.nrows = 1
630 649 self.nplots = 1
631 650 self.ylabel = 'Intensity [dB]'
632 651 self.titles = ['Scope']
633 652 self.colorbar = False
634 653 colspan = 3
635 654 rowspan = 1
636 655
637 656 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
638 657
639 658 yreal = y[channelIndexList,:].real
640 659 yimag = y[channelIndexList,:].imag
641 660 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
642 661 self.xlabel = "Range (Km)"
643 662 self.ylabel = "Intensity - IQ"
644 663
645 664 self.y = yreal
646 665 self.x = x
647 666 self.xmin = min(x)
648 667 self.xmax = max(x)
649 668
650 669
651 670 self.titles[0] = title
652 671
653 672 for i,ax in enumerate(self.axes):
654 673 title = "Channel %d" %(i)
655 674 if ax.firsttime:
656 675 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
657 676 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
658 677 else:
659 678 #pass
660 679 ax.plt_r.set_data(x, yreal[i,:])
661 680 ax.plt_i.set_data(x, yimag[i,:])
662 681
663 682 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
664 683 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
665 684 yreal = y.real
666 685 self.y = yreal
667 686 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
668 687 self.xlabel = "Range (Km)"
669 688 self.ylabel = "Intensity"
670 689 self.xmin = min(x)
671 690 self.xmax = max(x)
672 691
673 692
674 693 self.titles[0] = title
675 694
676 695 for i,ax in enumerate(self.axes):
677 696 title = "Channel %d" %(i)
678 697
679 698 ychannel = yreal[i,:]
680 699
681 700 if ax.firsttime:
682 701 ax.plt_r = ax.plot(x, ychannel)[0]
683 702 else:
684 703 #pass
685 704 ax.plt_r.set_data(x, ychannel)
686 705
687 706
688 707 def plot(self):
689 708
690 709 if self.channels:
691 710 channels = self.channels
692 711 else:
693 712 channels = self.data.channels
694 713
695 714
696 715
697 716 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
698 717
699 718 scope = self.data['scope']
700 719
701 720
702 721 if self.data.flagDataAsBlock:
703 722
704 723 for i in range(self.data.nProfiles):
705 724
706 725 wintitle1 = " [Profile = %d] " %i
707 726
708 727 if self.type == "power":
709 728 self.plot_power(self.data.heights,
710 729 scope[:,i,:],
711 730 channels,
712 731 thisDatetime,
713 732 wintitle1
714 733 )
715 734
716 735 if self.type == "iq":
717 736 self.plot_iq(self.data.heights,
718 737 scope[:,i,:],
719 738 channels,
720 739 thisDatetime,
721 740 wintitle1
722 741 )
723
724
725
726
727
728 742 else:
729 743 wintitle = " [Profile = %d] " %self.data.profileIndex
730 744
731 745 if self.type == "power":
732 746 self.plot_power(self.data.heights,
733 747 scope,
734 748 channels,
735 749 thisDatetime,
736 750 wintitle
737 751 )
738 752
739 753 if self.type == "iq":
740 754 self.plot_iq(self.data.heights,
741 755 scope,
742 756 channels,
743 757 thisDatetime,
744 758 wintitle
745 759 )
746
747
748 No newline at end of file
@@ -1,368 +1,408
1 1 '''
2 2 Created on Nov 9, 2016
3 3
4 4 @author: roj- LouVD
5 5 '''
6 6
7 7
8 8 import os
9 9 import sys
10 10 import time
11 11 import glob
12 12 import datetime
13 13
14 14 import numpy
15 15
16 16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator
17 17 from schainpy.model.data.jrodata import Parameters
18 18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
19 19 from schainpy.utils import log
20 20
21 21 FILE_HEADER_STRUCTURE = numpy.dtype([
22 22 ('FMN', '<u4'),
23 23 ('nrec', '<u4'),
24 24 ('fr_offset', '<u4'),
25 25 ('id', '<u4'),
26 26 ('site', 'u1', (32,))
27 27 ])
28 28
29 29 REC_HEADER_STRUCTURE = numpy.dtype([
30 30 ('rmn', '<u4'),
31 31 ('rcounter', '<u4'),
32 32 ('nr_offset', '<u4'),
33 33 ('tr_offset', '<u4'),
34 34 ('time', '<u4'),
35 35 ('time_msec', '<u4'),
36 36 ('tag', 'u1', (32,)),
37 37 ('comments', 'u1', (32,)),
38 38 ('lat', '<f4'),
39 39 ('lon', '<f4'),
40 40 ('gps_status', '<u4'),
41 41 ('freq', '<u4'),
42 42 ('freq0', '<u4'),
43 43 ('nchan', '<u4'),
44 44 ('delta_r', '<u4'),
45 45 ('nranges', '<u4'),
46 46 ('r0', '<u4'),
47 47 ('prf', '<u4'),
48 48 ('ncoh', '<u4'),
49 49 ('npoints', '<u4'),
50 50 ('polarization', '<i4'),
51 51 ('rx_filter', '<u4'),
52 52 ('nmodes', '<u4'),
53 53 ('dmode_index', '<u4'),
54 54 ('dmode_rngcorr', '<u4'),
55 55 ('nrxs', '<u4'),
56 56 ('acf_length', '<u4'),
57 57 ('acf_lags', '<u4'),
58 58 ('sea_to_atmos', '<f4'),
59 59 ('sea_notch', '<u4'),
60 60 ('lh_sea', '<u4'),
61 61 ('hh_sea', '<u4'),
62 62 ('nbins_sea', '<u4'),
63 63 ('min_snr', '<f4'),
64 64 ('min_cc', '<f4'),
65 65 ('max_time_diff', '<f4')
66 66 ])
67 67
68 68 DATA_STRUCTURE = numpy.dtype([
69 69 ('range', '<u4'),
70 70 ('status', '<u4'),
71 71 ('zonal', '<f4'),
72 72 ('meridional', '<f4'),
73 73 ('vertical', '<f4'),
74 74 ('zonal_a', '<f4'),
75 75 ('meridional_a', '<f4'),
76 76 ('corrected_fading', '<f4'), # seconds
77 77 ('uncorrected_fading', '<f4'), # seconds
78 78 ('time_diff', '<f4'),
79 79 ('major_axis', '<f4'),
80 80 ('axial_ratio', '<f4'),
81 81 ('orientation', '<f4'),
82 82 ('sea_power', '<u4'),
83 83 ('sea_algorithm', '<u4')
84 84 ])
85 85
86 86 @MPDecorator
87 87 class BLTRParamReader(JRODataReader, ProcessingUnit):
88 88 '''
89 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
89 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR
90 from *.sswma files
90 91 '''
91 92
92 93 ext = '.sswma'
93 94
94 95 def __init__(self):
95 96
96 97 ProcessingUnit.__init__(self)
97 98
98 99 self.dataOut = Parameters()
99 100 self.counter_records = 0
100 101 self.flagNoMoreFiles = 0
101 102 self.isConfig = False
102 103 self.filename = None
103 104
104 105 def setup(self,
105 106 path=None,
106 107 startDate=None,
107 108 endDate=None,
108 109 ext=None,
109 110 startTime=datetime.time(0, 0, 0),
110 111 endTime=datetime.time(23, 59, 59),
111 112 timezone=0,
112 113 status_value=0,
113 114 **kwargs):
114 115 self.path = path
115 116 self.startDate = startDate
116 117 self.endDate = endDate
117 118 self.startTime = startTime
118 119 self.endTime = endTime
119 120 self.status_value = status_value
120 121 self.datatime = datetime.datetime(1900,1,1)
122 self.delay = kwargs.get('delay', 10)
123 self.online = kwargs.get('online', False)
124 self.nTries = kwargs.get('nTries', 3)
121 125
122 126 if self.path is None:
123 127 raise ValueError("The path is not valid")
124 128
125 129 if ext is None:
126 130 ext = self.ext
127 131
128 self.search_files(self.path, startDate, endDate, ext)
132 self.fileList = self.search_files(self.path, startDate, endDate, ext)
129 133 self.timezone = timezone
130 134 self.fileIndex = 0
131 135
132 136 if not self.fileList:
133 137 raise Warning("There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' " % (
134 138 path))
135 139
136 140 self.setNextFile()
137 141
142 def search_last_file(self):
143 '''
144 Get last file and add it to the list
145 '''
146
147 for n in range(self.nTries+1):
148 if n>0:
149 log.warning(
150 "Waiting %0.2f seconds for the next file, try %03d ..." % (self.delay, n+1),
151 self.name
152 )
153 time.sleep(self.delay)
154 file_list = os.listdir(self.path)
155 file_list.sort()
156 if file_list:
157 if self.filename:
158 if file_list[-1] not in self.filename:
159 return file_list[-1]
160 else:
161 continue
162 return file_list[-1]
163 return 0
164
138 165 def search_files(self, path, startDate, endDate, ext):
139 166 '''
140 167 Searching for BLTR rawdata file in path
141 168 Creating a list of file to proces included in [startDate,endDate]
142 169
143 170 Input:
144 171 path - Path to find BLTR rawdata files
145 172 startDate - Select file from this date
146 173 enDate - Select file until this date
147 174 ext - Extension of the file to read
148 175 '''
149 176
150 177 log.success('Searching files in {} '.format(path), 'BLTRParamReader')
151 178 foldercounter = 0
152 179 fileList0 = glob.glob1(path, "*%s" % ext)
153 180 fileList0.sort()
154 181
155 self.fileList = []
156 self.dateFileList = []
157
158 182 for thisFile in fileList0:
159 183 year = thisFile[-14:-10]
160 184 if not isNumber(year):
161 185 continue
162 186
163 187 month = thisFile[-10:-8]
164 188 if not isNumber(month):
165 189 continue
166 190
167 191 day = thisFile[-8:-6]
168 192 if not isNumber(day):
169 193 continue
170 194
171 195 year, month, day = int(year), int(month), int(day)
172 196 dateFile = datetime.date(year, month, day)
173 197
174 198 if (startDate > dateFile) or (endDate < dateFile):
175 199 continue
176 200
177 self.fileList.append(thisFile)
178 self.dateFileList.append(dateFile)
201 yield thisFile
179 202
180 203 return
181 204
182 205 def setNextFile(self):
183 206
184 file_id = self.fileIndex
185
186 if file_id == len(self.fileList):
187 self.flagNoMoreFiles = 1
188 return 0
189
190 log.success('Opening {}'.format(self.fileList[file_id]), 'BLTRParamReader')
191 filename = os.path.join(self.path, self.fileList[file_id])
207 if self.online:
208 filename = self.search_last_file()
209 if not filename:
210 self.flagNoMoreFiles = 1
211 return 0
212 else:
213 try:
214 filename = next(self.fileList)
215 except StopIteration:
216 self.flagNoMoreFiles = 1
217 return 0
218
219 log.success('Opening {}'.format(filename), 'BLTRParamReader')
192 220
193 221 dirname, name = os.path.split(filename)
194 222 # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
195 self.siteFile = name.split('.')[0]
223 self.siteFile = filename.split('.')[0]
196 224 if self.filename is not None:
197 225 self.fp.close()
198 self.filename = filename
226 self.filename = os.path.join(self.path, filename)
199 227 self.fp = open(self.filename, 'rb')
200 228 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
201 229 self.nrecords = self.header_file['nrec'][0]
202 230 self.sizeOfFile = os.path.getsize(self.filename)
203 231 self.counter_records = 0
204 232 self.flagIsNewFile = 0
205 233 self.fileIndex += 1
234 time.sleep(2)
206 235
207 236 return 1
208 237
209 238 def readNextBlock(self):
210 239
211 240 while True:
212 if self.counter_records == self.nrecords:
241 if not self.online and self.counter_records == self.nrecords:
213 242 self.flagIsNewFile = 1
214 243 if not self.setNextFile():
215 244 return 0
216 245
217 self.readBlock()
246 try:
247 pointer = self.fp.tell()
248 self.readBlock()
249 except:
250 if self.online and self.waitDataBlock(pointer, 38512) == 1:
251 continue
252 else:
253 if not self.setNextFile():
254 return 0
218 255
219 256 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
220 257 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
221 258 log.warning(
222 259 'Reading Record No. {}/{} -> {} [Skipping]'.format(
223 260 self.counter_records,
224 261 self.nrecords,
225 262 self.datatime.ctime()),
226 263 'BLTRParamReader')
227 264 continue
228 265 break
229 266
230 log.log('Reading Record No. {}/{} -> {}'.format(
267 log.log('Reading Record No. {} -> {}'.format(
231 268 self.counter_records,
232 self.nrecords,
269 # self.nrecords,
233 270 self.datatime.ctime()), 'BLTRParamReader')
234 271
235 272 return 1
236 273
237 274 def readBlock(self):
238 275
239 276 pointer = self.fp.tell()
240 277 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
241 278 self.nchannels = int(header_rec['nchan'][0] / 2)
242 279 self.kchan = header_rec['nrxs'][0]
243 280 self.nmodes = header_rec['nmodes'][0]
244 281 self.nranges = header_rec['nranges'][0]
245 282 self.fp.seek(pointer)
246 283 self.height = numpy.empty((self.nmodes, self.nranges))
247 284 self.snr = numpy.empty((self.nmodes, int(self.nchannels), self.nranges))
248 285 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
249 286 self.flagDiscontinuousBlock = 0
250 287
251 288 for mode in range(self.nmodes):
252 289 self.readHeader()
253 290 data = self.readData()
254 291 self.height[mode] = (data[0] - self.correction) / 1000.
255 292 self.buffer[mode] = data[1]
256 293 self.snr[mode] = data[2]
257 294
258 295 self.counter_records = self.counter_records + self.nmodes
259 296
260 297 return
261 298
262 299 def readHeader(self):
263 300 '''
264 301 RecordHeader of BLTR rawdata file
265 302 '''
266 303
267 304 header_structure = numpy.dtype(
268 305 REC_HEADER_STRUCTURE.descr + [
269 306 ('antenna_coord', 'f4', (2, int(self.nchannels))),
270 307 ('rx_gains', 'u4', (int(self.nchannels),)),
271 308 ('rx_analysis', 'u4', (int(self.nchannels),))
272 309 ]
273 310 )
274 311
275 312 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
276 313 self.lat = self.header_rec['lat'][0]
277 314 self.lon = self.header_rec['lon'][0]
278 315 self.delta = self.header_rec['delta_r'][0]
279 316 self.correction = self.header_rec['dmode_rngcorr'][0]
280 317 self.imode = self.header_rec['dmode_index'][0]
281 318 self.antenna = self.header_rec['antenna_coord']
282 319 self.rx_gains = self.header_rec['rx_gains']
283 320 self.time = self.header_rec['time'][0]
284 321 dt = datetime.datetime.utcfromtimestamp(self.time)
285 322 if dt.date()>self.datatime.date():
286 323 self.flagDiscontinuousBlock = 1
287 324 self.datatime = dt
288 325
289 326 def readData(self):
290 327 '''
291 Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
328 Reading and filtering data block record of BLTR rawdata file,
329 filtering is according to status_value.
292 330
293 331 Input:
294 status_value - Array data is set to NAN for values that are not equal to status_value
332 status_value - Array data is set to NAN for values that are not
333 equal to status_value
295 334
296 335 '''
297 336 self.nchannels = int(self.nchannels)
298 337
299 338 data_structure = numpy.dtype(
300 339 DATA_STRUCTURE.descr + [
301 340 ('rx_saturation', 'u4', (self.nchannels,)),
302 341 ('chan_offset', 'u4', (2 * self.nchannels,)),
303 342 ('rx_amp', 'u4', (self.nchannels,)),
304 343 ('rx_snr', 'f4', (self.nchannels,)),
305 344 ('cross_snr', 'f4', (self.kchan,)),
306 345 ('sea_power_relative', 'f4', (self.kchan,))]
307 346 )
308 347
309 348 data = numpy.fromfile(self.fp, data_structure, self.nranges)
310 349
311 350 height = data['range']
312 351 winds = numpy.array(
313 352 (data['zonal'], data['meridional'], data['vertical']))
314 353 snr = data['rx_snr'].T
315 354
316 355 winds[numpy.where(winds == -9999.)] = numpy.nan
317 356 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
318 357 snr[numpy.where(snr == -9999.)] = numpy.nan
319 358 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
320 359 snr = numpy.power(10, snr / 10)
321 360
322 361 return height, winds, snr
323 362
324 363 def set_output(self):
325 364 '''
326 365 Storing data from databuffer to dataOut object
327 366 '''
328 367
329 368 self.dataOut.data_SNR = self.snr
330 369 self.dataOut.height = self.height
331 370 self.dataOut.data = self.buffer
332 371 self.dataOut.utctimeInit = self.time
333 372 self.dataOut.utctime = self.dataOut.utctimeInit
334 373 self.dataOut.useLocalTime = False
335 374 self.dataOut.paramInterval = 157
336 375 self.dataOut.timezone = self.timezone
337 376 self.dataOut.site = self.siteFile
338 377 self.dataOut.nrecords = self.nrecords / self.nmodes
339 378 self.dataOut.sizeOfFile = self.sizeOfFile
340 379 self.dataOut.lat = self.lat
341 380 self.dataOut.lon = self.lon
342 381 self.dataOut.channelList = list(range(self.nchannels))
343 382 self.dataOut.kchan = self.kchan
344 383 self.dataOut.delta = self.delta
345 384 self.dataOut.correction = self.correction
346 385 self.dataOut.nmodes = self.nmodes
347 386 self.dataOut.imode = self.imode
348 387 self.dataOut.antenna = self.antenna
349 388 self.dataOut.rx_gains = self.rx_gains
350 389 self.dataOut.flagNoData = False
351 390 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
352 391
353 392 def getData(self):
354 393 '''
355 394 Storing data from databuffer to dataOut object
356 395 '''
357 396 if self.flagNoMoreFiles:
358 397 self.dataOut.flagNoData = True
359 398 self.dataOut.error = 'No More files to read'
399 return
360 400
361 401 if not self.readNextBlock():
362 402 self.dataOut.flagNoData = True
363 return 0
403 self.dataOut.error = 'Time for wait new file reach!!!'
364 404
365 405 self.set_output()
366 406
367 407 return 1
368 408 No newline at end of file
@@ -1,1828 +1,1833
1 1 '''
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import sys
8 8 import glob
9 9 import time
10 10 import numpy
11 11 import fnmatch
12 12 import inspect
13 13 import time
14 14 import datetime
15 15 import traceback
16 16 import zmq
17 17
18 18 try:
19 19 from gevent import sleep
20 20 except:
21 21 from time import sleep
22 22
23 23 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
24 24 from schainpy.model.data.jroheaderIO import get_dtype_index, get_numpy_dtype, get_procflag_dtype, get_dtype_width
25 25 from schainpy.utils import log
26 26 import schainpy.admin
27 27
28 28 LOCALTIME = True
29 29
30 30
31 31 def isNumber(cad):
32 32 """
33 33 Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero.
34 34
35 35 Excepciones:
36 36 Si un determinado string no puede ser convertido a numero
37 37 Input:
38 38 str, string al cual se le analiza para determinar si convertible a un numero o no
39 39
40 40 Return:
41 41 True : si el string es uno numerico
42 42 False : no es un string numerico
43 43 """
44 44 try:
45 45 float(cad)
46 46 return True
47 47 except:
48 48 return False
49 49
50 50
51 51 def isFileInEpoch(filename, startUTSeconds, endUTSeconds):
52 52 """
53 53 Esta funcion determina si un archivo de datos se encuentra o no dentro del rango de fecha especificado.
54 54
55 55 Inputs:
56 56 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
57 57
58 58 startUTSeconds : fecha inicial del rango seleccionado. La fecha esta dada en
59 59 segundos contados desde 01/01/1970.
60 60 endUTSeconds : fecha final del rango seleccionado. La fecha esta dada en
61 61 segundos contados desde 01/01/1970.
62 62
63 63 Return:
64 64 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
65 65 fecha especificado, de lo contrario retorna False.
66 66
67 67 Excepciones:
68 68 Si el archivo no existe o no puede ser abierto
69 69 Si la cabecera no puede ser leida.
70 70
71 71 """
72 72 basicHeaderObj = BasicHeader(LOCALTIME)
73 73
74 74 try:
75 75 fp = open(filename, 'rb')
76 76 except IOError:
77 77 print("The file %s can't be opened" % (filename))
78 78 return 0
79 79
80 80 sts = basicHeaderObj.read(fp)
81 81 fp.close()
82 82
83 83 if not(sts):
84 84 print("Skipping the file %s because it has not a valid header" % (filename))
85 85 return 0
86 86
87 87 if not ((startUTSeconds <= basicHeaderObj.utc) and (endUTSeconds > basicHeaderObj.utc)):
88 88 return 0
89 89
90 90 return 1
91 91
92 92
93 93 def isTimeInRange(thisTime, startTime, endTime):
94 94 if endTime >= startTime:
95 95 if (thisTime < startTime) or (thisTime > endTime):
96 96 return 0
97 97 return 1
98 98 else:
99 99 if (thisTime < startTime) and (thisTime > endTime):
100 100 return 0
101 101 return 1
102 102
103 103
104 104 def isFileInTimeRange(filename, startDate, endDate, startTime, endTime):
105 105 """
106 106 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
107 107
108 108 Inputs:
109 109 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
110 110
111 111 startDate : fecha inicial del rango seleccionado en formato datetime.date
112 112
113 113 endDate : fecha final del rango seleccionado en formato datetime.date
114 114
115 115 startTime : tiempo inicial del rango seleccionado en formato datetime.time
116 116
117 117 endTime : tiempo final del rango seleccionado en formato datetime.time
118 118
119 119 Return:
120 120 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
121 121 fecha especificado, de lo contrario retorna False.
122 122
123 123 Excepciones:
124 124 Si el archivo no existe o no puede ser abierto
125 125 Si la cabecera no puede ser leida.
126 126
127 127 """
128 128
129 129 try:
130 130 fp = open(filename, 'rb')
131 131 except IOError:
132 132 print("The file %s can't be opened" % (filename))
133 133 return None
134 134
135 135 firstBasicHeaderObj = BasicHeader(LOCALTIME)
136 136 systemHeaderObj = SystemHeader()
137 137 radarControllerHeaderObj = RadarControllerHeader()
138 138 processingHeaderObj = ProcessingHeader()
139 139
140 140 lastBasicHeaderObj = BasicHeader(LOCALTIME)
141 141
142 142 sts = firstBasicHeaderObj.read(fp)
143 143
144 144 if not(sts):
145 145 print("[Reading] Skipping the file %s because it has not a valid header" % (filename))
146 146 return None
147 147
148 148 if not systemHeaderObj.read(fp):
149 149 return None
150 150
151 151 if not radarControllerHeaderObj.read(fp):
152 152 return None
153 153
154 154 if not processingHeaderObj.read(fp):
155 155 return None
156 156
157 157 filesize = os.path.getsize(filename)
158 158
159 159 offset = processingHeaderObj.blockSize + 24 # header size
160 160
161 161 if filesize <= offset:
162 162 print("[Reading] %s: This file has not enough data" % filename)
163 163 return None
164 164
165 165 fp.seek(-offset, 2)
166 166
167 167 sts = lastBasicHeaderObj.read(fp)
168 168
169 169 fp.close()
170 170
171 171 thisDatetime = lastBasicHeaderObj.datatime
172 172 thisTime_last_block = thisDatetime.time()
173 173
174 174 thisDatetime = firstBasicHeaderObj.datatime
175 175 thisDate = thisDatetime.date()
176 176 thisTime_first_block = thisDatetime.time()
177 177
178 178 # General case
179 179 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
180 180 #-----------o----------------------------o-----------
181 181 # startTime endTime
182 182
183 183 if endTime >= startTime:
184 184 if (thisTime_last_block < startTime) or (thisTime_first_block > endTime):
185 185 return None
186 186
187 187 return thisDatetime
188 188
189 189 # If endTime < startTime then endTime belongs to the next day
190 190
191 191 #<<<<<<<<<<<o o>>>>>>>>>>>
192 192 #-----------o----------------------------o-----------
193 193 # endTime startTime
194 194
195 195 if (thisDate == startDate) and (thisTime_last_block < startTime):
196 196 return None
197 197
198 198 if (thisDate == endDate) and (thisTime_first_block > endTime):
199 199 return None
200 200
201 201 if (thisTime_last_block < startTime) and (thisTime_first_block > endTime):
202 202 return None
203 203
204 204 return thisDatetime
205 205
206 206
207 207 def isFolderInDateRange(folder, startDate=None, endDate=None):
208 208 """
209 209 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
210 210
211 211 Inputs:
212 212 folder : nombre completo del directorio.
213 213 Su formato deberia ser "/path_root/?YYYYDDD"
214 214
215 215 siendo:
216 216 YYYY : Anio (ejemplo 2015)
217 217 DDD : Dia del anio (ejemplo 305)
218 218
219 219 startDate : fecha inicial del rango seleccionado en formato datetime.date
220 220
221 221 endDate : fecha final del rango seleccionado en formato datetime.date
222 222
223 223 Return:
224 224 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
225 225 fecha especificado, de lo contrario retorna False.
226 226 Excepciones:
227 227 Si el directorio no tiene el formato adecuado
228 228 """
229 229
230 230 basename = os.path.basename(folder)
231 231
232 232 if not isRadarFolder(basename):
233 233 print("The folder %s has not the rigth format" % folder)
234 234 return 0
235 235
236 236 if startDate and endDate:
237 237 thisDate = getDateFromRadarFolder(basename)
238 238
239 239 if thisDate < startDate:
240 240 return 0
241 241
242 242 if thisDate > endDate:
243 243 return 0
244 244
245 245 return 1
246 246
247 247
248 248 def isFileInDateRange(filename, startDate=None, endDate=None):
249 249 """
250 250 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
251 251
252 252 Inputs:
253 253 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
254 254
255 255 Su formato deberia ser "?YYYYDDDsss"
256 256
257 257 siendo:
258 258 YYYY : Anio (ejemplo 2015)
259 259 DDD : Dia del anio (ejemplo 305)
260 260 sss : set
261 261
262 262 startDate : fecha inicial del rango seleccionado en formato datetime.date
263 263
264 264 endDate : fecha final del rango seleccionado en formato datetime.date
265 265
266 266 Return:
267 267 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
268 268 fecha especificado, de lo contrario retorna False.
269 269 Excepciones:
270 270 Si el archivo no tiene el formato adecuado
271 271 """
272 272
273 273 basename = os.path.basename(filename)
274 274
275 275 if not isRadarFile(basename):
276 276 print("The filename %s has not the rigth format" % filename)
277 277 return 0
278 278
279 279 if startDate and endDate:
280 280 thisDate = getDateFromRadarFile(basename)
281 281
282 282 if thisDate < startDate:
283 283 return 0
284 284
285 285 if thisDate > endDate:
286 286 return 0
287 287
288 288 return 1
289 289
290 290
291 291 def getFileFromSet(path, ext, set):
292 292 validFilelist = []
293 293 fileList = os.listdir(path)
294 294
295 295 # 0 1234 567 89A BCDE
296 296 # H YYYY DDD SSS .ext
297 297
298 298 for thisFile in fileList:
299 299 try:
300 300 year = int(thisFile[1:5])
301 301 doy = int(thisFile[5:8])
302 302 except:
303 303 continue
304 304
305 305 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
306 306 continue
307 307
308 308 validFilelist.append(thisFile)
309 309
310 310 myfile = fnmatch.filter(
311 311 validFilelist, '*%4.4d%3.3d%3.3d*' % (year, doy, set))
312 312
313 313 if len(myfile) != 0:
314 314 return myfile[0]
315 315 else:
316 316 filename = '*%4.4d%3.3d%3.3d%s' % (year, doy, set, ext.lower())
317 317 print('the filename %s does not exist' % filename)
318 318 print('...going to the last file: ')
319 319
320 320 if validFilelist:
321 321 validFilelist = sorted(validFilelist, key=str.lower)
322 322 return validFilelist[-1]
323 323
324 324 return None
325 325
326 326
327 327 def getlastFileFromPath(path, ext):
328 328 """
329 329 Depura el fileList dejando solo los que cumplan el formato de "PYYYYDDDSSS.ext"
330 330 al final de la depuracion devuelve el ultimo file de la lista que quedo.
331 331
332 332 Input:
333 333 fileList : lista conteniendo todos los files (sin path) que componen una determinada carpeta
334 334 ext : extension de los files contenidos en una carpeta
335 335
336 336 Return:
337 337 El ultimo file de una determinada carpeta, no se considera el path.
338 338 """
339 339 validFilelist = []
340 340 fileList = os.listdir(path)
341 341
342 342 # 0 1234 567 89A BCDE
343 343 # H YYYY DDD SSS .ext
344 344
345 345 for thisFile in fileList:
346 346
347 347 year = thisFile[1:5]
348 348 if not isNumber(year):
349 349 continue
350 350
351 351 doy = thisFile[5:8]
352 352 if not isNumber(doy):
353 353 continue
354 354
355 355 year = int(year)
356 356 doy = int(doy)
357 357
358 358 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
359 359 continue
360 360
361 361 validFilelist.append(thisFile)
362 362
363 363 if validFilelist:
364 364 validFilelist = sorted(validFilelist, key=str.lower)
365 365 return validFilelist[-1]
366 366
367 367 return None
368 368
369 369
370 370 def checkForRealPath(path, foldercounter, year, doy, set, ext):
371 371 """
372 372 Por ser Linux Case Sensitive entonces checkForRealPath encuentra el nombre correcto de un path,
373 373 Prueba por varias combinaciones de nombres entre mayusculas y minusculas para determinar
374 374 el path exacto de un determinado file.
375 375
376 376 Example :
377 377 nombre correcto del file es .../.../D2009307/P2009307367.ext
378 378
379 379 Entonces la funcion prueba con las siguientes combinaciones
380 380 .../.../y2009307367.ext
381 381 .../.../Y2009307367.ext
382 382 .../.../x2009307/y2009307367.ext
383 383 .../.../x2009307/Y2009307367.ext
384 384 .../.../X2009307/y2009307367.ext
385 385 .../.../X2009307/Y2009307367.ext
386 386 siendo para este caso, la ultima combinacion de letras, identica al file buscado
387 387
388 388 Return:
389 389 Si encuentra la cobinacion adecuada devuelve el path completo y el nombre del file
390 390 caso contrario devuelve None como path y el la ultima combinacion de nombre en mayusculas
391 391 para el filename
392 392 """
393 393 fullfilename = None
394 394 find_flag = False
395 395 filename = None
396 396
397 397 prefixDirList = [None, 'd', 'D']
398 398 if ext.lower() == ".r": # voltage
399 399 prefixFileList = ['d', 'D']
400 400 elif ext.lower() == ".pdata": # spectra
401 401 prefixFileList = ['p', 'P']
402 402 else:
403 403 return None, filename
404 404
405 405 # barrido por las combinaciones posibles
406 406 for prefixDir in prefixDirList:
407 407 thispath = path
408 408 if prefixDir != None:
409 409 # formo el nombre del directorio xYYYYDDD (x=d o x=D)
410 410 if foldercounter == 0:
411 411 thispath = os.path.join(path, "%s%04d%03d" %
412 412 (prefixDir, year, doy))
413 413 else:
414 414 thispath = os.path.join(path, "%s%04d%03d_%02d" % (
415 415 prefixDir, year, doy, foldercounter))
416 416 for prefixFile in prefixFileList: # barrido por las dos combinaciones posibles de "D"
417 417 # formo el nombre del file xYYYYDDDSSS.ext
418 418 filename = "%s%04d%03d%03d%s" % (prefixFile, year, doy, set, ext)
419 419 fullfilename = os.path.join(
420 420 thispath, filename) # formo el path completo
421 421
422 422 if os.path.exists(fullfilename): # verifico que exista
423 423 find_flag = True
424 424 break
425 425 if find_flag:
426 426 break
427 427
428 428 if not(find_flag):
429 429 return None, filename
430 430
431 431 return fullfilename, filename
432 432
433 433
434 434 def isRadarFolder(folder):
435 435 try:
436 436 year = int(folder[1:5])
437 437 doy = int(folder[5:8])
438 438 except:
439 439 return 0
440 440
441 441 return 1
442 442
443 443
444 444 def isRadarFile(file):
445 445 try:
446 446 year = int(file[1:5])
447 447 doy = int(file[5:8])
448 448 set = int(file[8:11])
449 449 except:
450 450 return 0
451 451
452 452 return 1
453 453
454 454
455 455 def getDateFromRadarFile(file):
456 456 try:
457 457 year = int(file[1:5])
458 458 doy = int(file[5:8])
459 459 set = int(file[8:11])
460 460 except:
461 461 return None
462 462
463 463 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
464 464 return thisDate
465 465
466 466
467 467 def getDateFromRadarFolder(folder):
468 468 try:
469 469 year = int(folder[1:5])
470 470 doy = int(folder[5:8])
471 471 except:
472 472 return None
473 473
474 474 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
475 475 return thisDate
476 476
477 477
478 478 class JRODataIO:
479 479
480 480 c = 3E8
481 481
482 482 isConfig = False
483 483
484 484 basicHeaderObj = None
485 485
486 486 systemHeaderObj = None
487 487
488 488 radarControllerHeaderObj = None
489 489
490 490 processingHeaderObj = None
491 491
492 492 dtype = None
493 493
494 494 pathList = []
495 495
496 496 filenameList = []
497 497
498 498 filename = None
499 499
500 500 ext = None
501 501
502 502 flagIsNewFile = 1
503 503
504 504 flagDiscontinuousBlock = 0
505 505
506 506 flagIsNewBlock = 0
507 507
508 508 fp = None
509 509
510 510 firstHeaderSize = 0
511 511
512 512 basicHeaderSize = 24
513 513
514 514 versionFile = 1103
515 515
516 516 fileSize = None
517 517
518 518 # ippSeconds = None
519 519
520 520 fileSizeByHeader = None
521 521
522 522 fileIndex = None
523 523
524 524 profileIndex = None
525 525
526 526 blockIndex = None
527 527
528 528 nTotalBlocks = None
529 529
530 530 maxTimeStep = 30
531 531
532 532 lastUTTime = None
533 533
534 534 datablock = None
535 535
536 536 dataOut = None
537 537
538 538 blocksize = None
539 539
540 540 getByBlock = False
541 541
542 542 def __init__(self):
543 543
544 544 raise NotImplementedError
545 545
546 546 def run(self):
547 547
548 548 raise NotImplementedError
549 549
550 550 def getDtypeWidth(self):
551 551
552 552 dtype_index = get_dtype_index(self.dtype)
553 553 dtype_width = get_dtype_width(dtype_index)
554 554
555 555 return dtype_width
556 556
557 557 def getAllowedArgs(self):
558 558 if hasattr(self, '__attrs__'):
559 559 return self.__attrs__
560 560 else:
561 561 return inspect.getargspec(self.run).args
562 562
563 563
564 564 class JRODataReader(JRODataIO):
565 565
566 566 online = 0
567 567
568 568 realtime = 0
569 569
570 570 nReadBlocks = 0
571 571
572 572 delay = 10 # number of seconds waiting a new file
573 573
574 574 nTries = 3 # quantity tries
575 575
576 576 nFiles = 3 # number of files for searching
577 577
578 578 path = None
579 579
580 580 foldercounter = 0
581 581
582 582 flagNoMoreFiles = 0
583 583
584 584 datetimeList = []
585 585
586 586 __isFirstTimeOnline = 1
587 587
588 588 __printInfo = True
589 589
590 590 profileIndex = None
591 591
592 592 nTxs = 1
593 593
594 594 txIndex = None
595 595
596 596 # Added--------------------
597 597
598 598 selBlocksize = None
599 599
600 600 selBlocktime = None
601 601
602 602 def __init__(self):
603 603 """
604 604 This class is used to find data files
605 605
606 606 Example:
607 607 reader = JRODataReader()
608 608 fileList = reader.findDataFiles()
609 609
610 610 """
611 611 pass
612 612
613 613 def createObjByDefault(self):
614 614 """
615 615
616 616 """
617 617 raise NotImplementedError
618 618
619 619 def getBlockDimension(self):
620 620
621 621 raise NotImplementedError
622 622
623 623 def searchFilesOffLine(self,
624 624 path,
625 625 startDate=None,
626 626 endDate=None,
627 627 startTime=datetime.time(0, 0, 0),
628 628 endTime=datetime.time(23, 59, 59),
629 629 set=None,
630 630 expLabel='',
631 631 ext='.r',
632 632 cursor=None,
633 633 skip=None,
634 634 walk=True):
635 635
636 636 self.filenameList = []
637 637 self.datetimeList = []
638 638
639 639 pathList = []
640 640
641 641 dateList, pathList = self.findDatafiles(
642 642 path, startDate, endDate, expLabel, ext, walk, include_path=True)
643 643
644 644 if dateList == []:
645 645 return [], []
646 646
647 647 if len(dateList) > 1:
648 648 print("[Reading] Data found for date range [%s - %s]: total days = %d" % (startDate, endDate, len(dateList)))
649 649 else:
650 650 print("[Reading] Data found for date range [%s - %s]: date = %s" % (startDate, endDate, dateList[0]))
651 651
652 652 filenameList = []
653 653 datetimeList = []
654 654
655 655 for thisPath in pathList:
656 656
657 657 fileList = glob.glob1(thisPath, "*%s" % ext)
658 658 fileList.sort()
659 659
660 660 for file in fileList:
661 661
662 662 filename = os.path.join(thisPath, file)
663 663
664 664 if not isFileInDateRange(filename, startDate, endDate):
665 665 continue
666 666
667 667 thisDatetime = isFileInTimeRange(
668 668 filename, startDate, endDate, startTime, endTime)
669 669
670 670 if not(thisDatetime):
671 671 continue
672 672
673 673 filenameList.append(filename)
674 674 datetimeList.append(thisDatetime)
675 675
676 676 if cursor is not None and skip is not None:
677 677 filenameList = filenameList[cursor * skip:cursor * skip + skip]
678 678 datetimeList = datetimeList[cursor * skip:cursor * skip + skip]
679 679
680 680 if not(filenameList):
681 681 print("[Reading] Time range selected invalid [%s - %s]: No *%s files in %s)" % (startTime, endTime, ext, path))
682 682 return [], []
683 683
684 684 print("[Reading] %d file(s) was(were) found in time range: %s - %s" % (len(filenameList), startTime, endTime))
685 685
686 686 # for i in range(len(filenameList)):
687 687 # print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
688 688
689 689 self.filenameList = filenameList
690 690 self.datetimeList = datetimeList
691 691
692 692 return pathList, filenameList
693 693
694 694 def __searchFilesOnLine(self, path, expLabel="", ext=None, walk=True, set=None):
695 695 """
696 696 Busca el ultimo archivo de la ultima carpeta (determinada o no por startDateTime) y
697 697 devuelve el archivo encontrado ademas de otros datos.
698 698
699 699 Input:
700 700 path : carpeta donde estan contenidos los files que contiene data
701 701
702 702 expLabel : Nombre del subexperimento (subfolder)
703 703
704 704 ext : extension de los files
705 705
706 706 walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath)
707 707
708 708 Return:
709 709 directory : eL directorio donde esta el file encontrado
710 710 filename : el ultimo file de una determinada carpeta
711 711 year : el anho
712 712 doy : el numero de dia del anho
713 713 set : el set del archivo
714 714
715 715
716 716 """
717 717 if not os.path.isdir(path):
718 718 return None, None, None, None, None, None
719 719
720 720 dirList = []
721 721
722 722 if not walk:
723 723 fullpath = path
724 724 foldercounter = 0
725 725 else:
726 726 # Filtra solo los directorios
727 727 for thisPath in os.listdir(path):
728 728 if not os.path.isdir(os.path.join(path, thisPath)):
729 729 continue
730 730 if not isRadarFolder(thisPath):
731 731 continue
732 732
733 733 dirList.append(thisPath)
734 734
735 735 if not(dirList):
736 736 return None, None, None, None, None, None
737 737
738 738 dirList = sorted(dirList, key=str.lower)
739 739
740 740 doypath = dirList[-1]
741 741 foldercounter = int(doypath.split('_')[1]) if len(
742 742 doypath.split('_')) > 1 else 0
743 743 fullpath = os.path.join(path, doypath, expLabel)
744 744
745 745 print("[Reading] %s folder was found: " % (fullpath))
746 746
747 747 if set == None:
748 748 filename = getlastFileFromPath(fullpath, ext)
749 749 else:
750 750 filename = getFileFromSet(fullpath, ext, set)
751 751
752 752 if not(filename):
753 753 return None, None, None, None, None, None
754 754
755 755 print("[Reading] %s file was found" % (filename))
756 756
757 757 if not(self.__verifyFile(os.path.join(fullpath, filename))):
758 758 return None, None, None, None, None, None
759 759
760 760 year = int(filename[1:5])
761 761 doy = int(filename[5:8])
762 762 set = int(filename[8:11])
763 763
764 764 return fullpath, foldercounter, filename, year, doy, set
765 765
766 766 def __setNextFileOffline(self):
767 767
768 768 idFile = self.fileIndex
769 769
770 770 while (True):
771 771 idFile += 1
772 772 if not(idFile < len(self.filenameList)):
773 773 self.flagNoMoreFiles = 1
774 774 # print "[Reading] No more Files"
775 775 return 0
776 776
777 777 filename = self.filenameList[idFile]
778 778
779 779 if not(self.__verifyFile(filename)):
780 780 continue
781 781
782 782 fileSize = os.path.getsize(filename)
783 783 fp = open(filename, 'rb')
784 784 break
785 785
786 786 self.flagIsNewFile = 1
787 787 self.fileIndex = idFile
788 788 self.filename = filename
789 789 self.fileSize = fileSize
790 790 self.fp = fp
791 791
792 792 # print "[Reading] Setting the file: %s"%self.filename
793 793
794 794 return 1
795 795
796 796 def __setNextFileOnline(self):
797 797 """
798 798 Busca el siguiente file que tenga suficiente data para ser leida, dentro de un folder especifico, si
799 799 no encuentra un file valido espera un tiempo determinado y luego busca en los posibles n files
800 800 siguientes.
801 801
802 802 Affected:
803 803 self.flagIsNewFile
804 804 self.filename
805 805 self.fileSize
806 806 self.fp
807 807 self.set
808 808 self.flagNoMoreFiles
809 809
810 810 Return:
811 811 0 : si luego de una busqueda del siguiente file valido este no pudo ser encontrado
812 812 1 : si el file fue abierto con exito y esta listo a ser leido
813 813
814 814 Excepciones:
815 815 Si un determinado file no puede ser abierto
816 816 """
817 817 nFiles = 0
818 818 fileOk_flag = False
819 819 firstTime_flag = True
820 820
821 821 self.set += 1
822 822
823 823 if self.set > 999:
824 824 self.set = 0
825 825 self.foldercounter += 1
826 826
827 827 # busca el 1er file disponible
828 828 fullfilename, filename = checkForRealPath(
829 829 self.path, self.foldercounter, self.year, self.doy, self.set, self.ext)
830 830 if fullfilename:
831 831 if self.__verifyFile(fullfilename, False):
832 832 fileOk_flag = True
833 833
834 834 # si no encuentra un file entonces espera y vuelve a buscar
835 835 if not(fileOk_flag):
836 836 # busco en los siguientes self.nFiles+1 files posibles
837 837 for nFiles in range(self.nFiles + 1):
838 838
839 839 if firstTime_flag: # si es la 1era vez entonces hace el for self.nTries veces
840 840 tries = self.nTries
841 841 else:
842 842 tries = 1 # si no es la 1era vez entonces solo lo hace una vez
843 843
844 844 for nTries in range(tries):
845 845 if firstTime_flag:
846 846 print("\t[Reading] Waiting %0.2f sec for the next file: \"%s\" , try %03d ..." % (self.delay, filename, nTries + 1))
847 847 sleep(self.delay)
848 848 else:
849 849 print("\t[Reading] Searching the next \"%s%04d%03d%03d%s\" file ..." % (self.optchar, self.year, self.doy, self.set, self.ext))
850 850
851 851 fullfilename, filename = checkForRealPath(
852 852 self.path, self.foldercounter, self.year, self.doy, self.set, self.ext)
853 853 if fullfilename:
854 854 if self.__verifyFile(fullfilename):
855 855 fileOk_flag = True
856 856 break
857 857
858 858 if fileOk_flag:
859 859 break
860 860
861 861 firstTime_flag = False
862 862
863 863 log.warning('Skipping the file {} due to this file doesn\'t exist'.format(filename))
864 864 self.set += 1
865 865
866 866 # si no encuentro el file buscado cambio de carpeta y busco en la siguiente carpeta
867 867 if nFiles == (self.nFiles - 1):
868 868 self.set = 0
869 869 self.doy += 1
870 870 self.foldercounter = 0
871 871
872 872 if fileOk_flag:
873 873 self.fileSize = os.path.getsize(fullfilename)
874 874 self.filename = fullfilename
875 875 self.flagIsNewFile = 1
876 876 if self.fp != None:
877 877 self.fp.close()
878 878 self.fp = open(fullfilename, 'rb')
879 879 self.flagNoMoreFiles = 0
880 880 # print '[Reading] Setting the file: %s' % fullfilename
881 881 else:
882 882 self.fileSize = 0
883 883 self.filename = None
884 884 self.flagIsNewFile = 0
885 885 self.fp = None
886 886 self.flagNoMoreFiles = 1
887 887 # print '[Reading] No more files to read'
888 888
889 889 return fileOk_flag
890 890
891 891 def setNextFile(self):
892 892 if self.fp != None:
893 893 self.fp.close()
894 894
895 895 if self.online:
896 896 newFile = self.__setNextFileOnline()
897 897 else:
898 898 newFile = self.__setNextFileOffline()
899 899
900 900 if not(newFile):
901 901 self.dataOut.error = 'No more files to read'
902 902 return 0
903 903
904 904 if self.verbose:
905 905 print('[Reading] Setting the file: %s' % self.filename)
906 906
907 907 self.__readFirstHeader()
908 908 self.nReadBlocks = 0
909 909 return 1
910 910
911 911 def __waitNewBlock(self):
912 912 """
913 913 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
914 914
915 915 Si el modo de lectura es OffLine siempre retorn 0
916 916 """
917 917 if not self.online:
918 918 return 0
919 919
920 920 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
921 921 return 0
922 922
923 923 currentPointer = self.fp.tell()
924 924
925 925 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
926 926
927 927 for nTries in range(self.nTries):
928 928
929 929 self.fp.close()
930 930 self.fp = open(self.filename, 'rb')
931 931 self.fp.seek(currentPointer)
932 932
933 933 self.fileSize = os.path.getsize(self.filename)
934 934 currentSize = self.fileSize - currentPointer
935 935
936 936 if (currentSize >= neededSize):
937 937 self.basicHeaderObj.read(self.fp)
938 938 return 1
939 939
940 940 if self.fileSize == self.fileSizeByHeader:
941 941 # self.flagEoF = True
942 942 return 0
943 943
944 944 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
945 945 sleep(self.delay)
946 946
947 947 return 0
948 948
949 def waitDataBlock(self, pointer_location):
949 def waitDataBlock(self, pointer_location, blocksize=None):
950 950
951 951 currentPointer = pointer_location
952
953 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
952 if blocksize is None:
953 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
954 else:
955 neededSize = blocksize
954 956
955 957 for nTries in range(self.nTries):
956 958 self.fp.close()
957 959 self.fp = open(self.filename, 'rb')
958 960 self.fp.seek(currentPointer)
959 961
960 962 self.fileSize = os.path.getsize(self.filename)
961 963 currentSize = self.fileSize - currentPointer
962 964
963 965 if (currentSize >= neededSize):
964 966 return 1
965 967
966 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
968 log.warning(
969 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
970 self.name
971 )
967 972 sleep(self.delay)
968 973
969 974 return 0
970 975
971 976 def __jumpToLastBlock(self):
972 977
973 978 if not(self.__isFirstTimeOnline):
974 979 return
975 980
976 981 csize = self.fileSize - self.fp.tell()
977 982 blocksize = self.processingHeaderObj.blockSize
978 983
979 984 # salta el primer bloque de datos
980 985 if csize > self.processingHeaderObj.blockSize:
981 986 self.fp.seek(self.fp.tell() + blocksize)
982 987 else:
983 988 return
984 989
985 990 csize = self.fileSize - self.fp.tell()
986 991 neededsize = self.processingHeaderObj.blockSize + self.basicHeaderSize
987 992 while True:
988 993
989 994 if self.fp.tell() < self.fileSize:
990 995 self.fp.seek(self.fp.tell() + neededsize)
991 996 else:
992 997 self.fp.seek(self.fp.tell() - neededsize)
993 998 break
994 999
995 1000 # csize = self.fileSize - self.fp.tell()
996 1001 # neededsize = self.processingHeaderObj.blockSize + self.basicHeaderSize
997 1002 # factor = int(csize/neededsize)
998 1003 # if factor > 0:
999 1004 # self.fp.seek(self.fp.tell() + factor*neededsize)
1000 1005
1001 1006 self.flagIsNewFile = 0
1002 1007 self.__isFirstTimeOnline = 0
1003 1008
1004 1009 def __setNewBlock(self):
1005 1010 # if self.server is None:
1006 1011 if self.fp == None:
1007 1012 return 0
1008 1013
1009 1014 # if self.online:
1010 1015 # self.__jumpToLastBlock()
1011 1016
1012 1017 if self.flagIsNewFile:
1013 1018 self.lastUTTime = self.basicHeaderObj.utc
1014 1019 return 1
1015 1020
1016 1021 if self.realtime:
1017 1022 self.flagDiscontinuousBlock = 1
1018 1023 if not(self.setNextFile()):
1019 1024 return 0
1020 1025 else:
1021 1026 return 1
1022 1027 # if self.server is None:
1023 1028 currentSize = self.fileSize - self.fp.tell()
1024 1029 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
1025 1030 if (currentSize >= neededSize):
1026 1031 self.basicHeaderObj.read(self.fp)
1027 1032 self.lastUTTime = self.basicHeaderObj.utc
1028 1033 return 1
1029 1034 # else:
1030 1035 # self.basicHeaderObj.read(self.zHeader)
1031 1036 # self.lastUTTime = self.basicHeaderObj.utc
1032 1037 # return 1
1033 1038 if self.__waitNewBlock():
1034 1039 self.lastUTTime = self.basicHeaderObj.utc
1035 1040 return 1
1036 1041 # if self.server is None:
1037 1042 if not(self.setNextFile()):
1038 1043 return 0
1039 1044
1040 1045 deltaTime = self.basicHeaderObj.utc - self.lastUTTime
1041 1046 self.lastUTTime = self.basicHeaderObj.utc
1042 1047
1043 1048 self.flagDiscontinuousBlock = 0
1044 1049
1045 1050 if deltaTime > self.maxTimeStep:
1046 1051 self.flagDiscontinuousBlock = 1
1047 1052
1048 1053 return 1
1049 1054
1050 1055 def readNextBlock(self):
1051 1056
1052 1057 # Skip block out of startTime and endTime
1053 1058 while True:
1054 1059 if not(self.__setNewBlock()):
1055 1060 self.dataOut.error = 'No more files to read'
1056 1061 return 0
1057 1062
1058 1063 if not(self.readBlock()):
1059 1064 return 0
1060 1065
1061 1066 self.getBasicHeader()
1062 1067 if (self.dataOut.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or (self.dataOut.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
1063 1068 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.nReadBlocks,
1064 1069 self.processingHeaderObj.dataBlocksPerFile,
1065 1070 self.dataOut.datatime.ctime()))
1066 1071 continue
1067 1072
1068 1073 break
1069 1074
1070 1075 if self.verbose:
1071 1076 print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
1072 1077 self.processingHeaderObj.dataBlocksPerFile,
1073 1078 self.dataOut.datatime.ctime()))
1074 1079 return 1
1075 1080
1076 1081 def __readFirstHeader(self):
1077 1082
1078 1083 self.basicHeaderObj.read(self.fp)
1079 1084 self.systemHeaderObj.read(self.fp)
1080 1085 self.radarControllerHeaderObj.read(self.fp)
1081 1086 self.processingHeaderObj.read(self.fp)
1082 1087
1083 1088 self.firstHeaderSize = self.basicHeaderObj.size
1084 1089
1085 1090 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
1086 1091 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
1087 1092 if datatype == 0:
1088 1093 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
1089 1094 elif datatype == 1:
1090 1095 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
1091 1096 elif datatype == 2:
1092 1097 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
1093 1098 elif datatype == 3:
1094 1099 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
1095 1100 elif datatype == 4:
1096 1101 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
1097 1102 elif datatype == 5:
1098 1103 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
1099 1104 else:
1100 1105 raise ValueError('Data type was not defined')
1101 1106
1102 1107 self.dtype = datatype_str
1103 1108 #self.ippSeconds = 2 * 1000 * self.radarControllerHeaderObj.ipp / self.c
1104 1109 self.fileSizeByHeader = self.processingHeaderObj.dataBlocksPerFile * self.processingHeaderObj.blockSize + \
1105 1110 self.firstHeaderSize + self.basicHeaderSize * \
1106 1111 (self.processingHeaderObj.dataBlocksPerFile - 1)
1107 1112 # self.dataOut.channelList = numpy.arange(self.systemHeaderObj.numChannels)
1108 1113 # self.dataOut.channelIndexList = numpy.arange(self.systemHeaderObj.numChannels)
1109 1114 self.getBlockDimension()
1110 1115
1111 1116 def __verifyFile(self, filename, msgFlag=True):
1112 1117
1113 1118 msg = None
1114 1119
1115 1120 try:
1116 1121 fp = open(filename, 'rb')
1117 1122 except IOError:
1118 1123
1119 1124 if msgFlag:
1120 1125 print("[Reading] File %s can't be opened" % (filename))
1121 1126
1122 1127 return False
1123 1128
1124 1129 currentPosition = fp.tell()
1125 1130 neededSize = self.processingHeaderObj.blockSize + self.firstHeaderSize
1126 1131
1127 1132 if neededSize == 0:
1128 1133 basicHeaderObj = BasicHeader(LOCALTIME)
1129 1134 systemHeaderObj = SystemHeader()
1130 1135 radarControllerHeaderObj = RadarControllerHeader()
1131 1136 processingHeaderObj = ProcessingHeader()
1132 1137
1133 1138 if not(basicHeaderObj.read(fp)):
1134 1139 fp.close()
1135 1140 return False
1136 1141
1137 1142 if not(systemHeaderObj.read(fp)):
1138 1143 fp.close()
1139 1144 return False
1140 1145
1141 1146 if not(radarControllerHeaderObj.read(fp)):
1142 1147 fp.close()
1143 1148 return False
1144 1149
1145 1150 if not(processingHeaderObj.read(fp)):
1146 1151 fp.close()
1147 1152 return False
1148 1153
1149 1154 neededSize = processingHeaderObj.blockSize + basicHeaderObj.size
1150 1155 else:
1151 1156 msg = "[Reading] Skipping the file %s due to it hasn't enough data" % filename
1152 1157
1153 1158 fp.close()
1154 1159
1155 1160 fileSize = os.path.getsize(filename)
1156 1161 currentSize = fileSize - currentPosition
1157 1162
1158 1163 if currentSize < neededSize:
1159 1164 if msgFlag and (msg != None):
1160 1165 print(msg)
1161 1166 return False
1162 1167
1163 1168 return True
1164 1169
1165 1170 def findDatafiles(self, path, startDate=None, endDate=None, expLabel='', ext='.r', walk=True, include_path=False):
1166 1171
1167 1172 path_empty = True
1168 1173
1169 1174 dateList = []
1170 1175 pathList = []
1171 1176
1172 1177 multi_path = path.split(',')
1173 1178
1174 1179 if not walk:
1175 1180
1176 1181 for single_path in multi_path:
1177 1182
1178 1183 if not os.path.isdir(single_path):
1179 1184 continue
1180 1185
1181 1186 fileList = glob.glob1(single_path, "*" + ext)
1182 1187
1183 1188 if not fileList:
1184 1189 continue
1185 1190
1186 1191 path_empty = False
1187 1192
1188 1193 fileList.sort()
1189 1194
1190 1195 for thisFile in fileList:
1191 1196
1192 1197 if not os.path.isfile(os.path.join(single_path, thisFile)):
1193 1198 continue
1194 1199
1195 1200 if not isRadarFile(thisFile):
1196 1201 continue
1197 1202
1198 1203 if not isFileInDateRange(thisFile, startDate, endDate):
1199 1204 continue
1200 1205
1201 1206 thisDate = getDateFromRadarFile(thisFile)
1202 1207
1203 1208 if thisDate in dateList:
1204 1209 continue
1205 1210
1206 1211 dateList.append(thisDate)
1207 1212 pathList.append(single_path)
1208 1213
1209 1214 else:
1210 1215 for single_path in multi_path:
1211 1216
1212 1217 if not os.path.isdir(single_path):
1213 1218 continue
1214 1219
1215 1220 dirList = []
1216 1221
1217 1222 for thisPath in os.listdir(single_path):
1218 1223
1219 1224 if not os.path.isdir(os.path.join(single_path, thisPath)):
1220 1225 continue
1221 1226
1222 1227 if not isRadarFolder(thisPath):
1223 1228 continue
1224 1229
1225 1230 if not isFolderInDateRange(thisPath, startDate, endDate):
1226 1231 continue
1227 1232
1228 1233 dirList.append(thisPath)
1229 1234
1230 1235 if not dirList:
1231 1236 continue
1232 1237
1233 1238 dirList.sort()
1234 1239
1235 1240 for thisDir in dirList:
1236 1241
1237 1242 datapath = os.path.join(single_path, thisDir, expLabel)
1238 1243 fileList = glob.glob1(datapath, "*" + ext)
1239 1244
1240 1245 if not fileList:
1241 1246 continue
1242 1247
1243 1248 path_empty = False
1244 1249
1245 1250 thisDate = getDateFromRadarFolder(thisDir)
1246 1251
1247 1252 pathList.append(datapath)
1248 1253 dateList.append(thisDate)
1249 1254
1250 1255 dateList.sort()
1251 1256
1252 1257 if walk:
1253 1258 pattern_path = os.path.join(multi_path[0], "[dYYYYDDD]", expLabel)
1254 1259 else:
1255 1260 pattern_path = multi_path[0]
1256 1261
1257 1262 if path_empty:
1258 1263 print("[Reading] No *%s files in %s for %s to %s" % (ext, pattern_path, startDate, endDate))
1259 1264 else:
1260 1265 if not dateList:
1261 1266 print("[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" % (startDate, endDate, ext, path))
1262 1267
1263 1268 if include_path:
1264 1269 return dateList, pathList
1265 1270
1266 1271 return dateList
1267 1272
1268 1273 def setup(self,
1269 1274 path=None,
1270 1275 startDate=None,
1271 1276 endDate=None,
1272 1277 startTime=datetime.time(0, 0, 0),
1273 1278 endTime=datetime.time(23, 59, 59),
1274 1279 set=None,
1275 1280 expLabel="",
1276 1281 ext=None,
1277 1282 online=False,
1278 1283 delay=60,
1279 1284 walk=True,
1280 1285 getblock=False,
1281 1286 nTxs=1,
1282 1287 realtime=False,
1283 1288 blocksize=None,
1284 1289 blocktime=None,
1285 1290 skip=None,
1286 1291 cursor=None,
1287 1292 warnings=True,
1288 1293 verbose=True,
1289 1294 server=None,
1290 1295 format=None,
1291 1296 oneDDict=None,
1292 1297 twoDDict=None,
1293 1298 independentParam=None):
1294 1299 if server is not None:
1295 1300 if 'tcp://' in server:
1296 1301 address = server
1297 1302 else:
1298 1303 address = 'ipc:///tmp/%s' % server
1299 1304 self.server = address
1300 1305 self.context = zmq.Context()
1301 1306 self.receiver = self.context.socket(zmq.PULL)
1302 1307 self.receiver.connect(self.server)
1303 1308 time.sleep(0.5)
1304 1309 print('[Starting] ReceiverData from {}'.format(self.server))
1305 1310 else:
1306 1311 self.server = None
1307 1312 if path == None:
1308 1313 raise ValueError("[Reading] The path is not valid")
1309 1314
1310 1315 if ext == None:
1311 1316 ext = self.ext
1312 1317
1313 1318 if online:
1314 1319 print("[Reading] Searching files in online mode...")
1315 1320
1316 1321 for nTries in range(self.nTries):
1317 1322 fullpath, foldercounter, file, year, doy, set = self.__searchFilesOnLine(
1318 1323 path=path, expLabel=expLabel, ext=ext, walk=walk, set=set)
1319 1324
1320 1325 if fullpath:
1321 1326 break
1322 1327
1323 1328 print('[Reading] Waiting %0.2f sec for an valid file in %s: try %02d ...' % (self.delay, path, nTries + 1))
1324 1329 sleep(self.delay)
1325 1330
1326 1331 if not(fullpath):
1327 1332 self.dataOut.error = 'There isn\'t any valid file in {}'.format(path)
1328 1333 return
1329 1334
1330 1335 self.year = year
1331 1336 self.doy = doy
1332 1337 self.set = set - 1
1333 1338 self.path = path
1334 1339 self.foldercounter = foldercounter
1335 1340 last_set = None
1336 1341 else:
1337 1342 print("[Reading] Searching files in offline mode ...")
1338 1343 pathList, filenameList = self.searchFilesOffLine(path, startDate=startDate, endDate=endDate,
1339 1344 startTime=startTime, endTime=endTime,
1340 1345 set=set, expLabel=expLabel, ext=ext,
1341 1346 walk=walk, cursor=cursor,
1342 1347 skip=skip)
1343 1348
1344 1349 if not(pathList):
1345 1350 self.fileIndex = -1
1346 1351 self.pathList = []
1347 1352 self.filenameList = []
1348 1353 return
1349 1354
1350 1355 self.fileIndex = -1
1351 1356 self.pathList = pathList
1352 1357 self.filenameList = filenameList
1353 1358 file_name = os.path.basename(filenameList[-1])
1354 1359 basename, ext = os.path.splitext(file_name)
1355 1360 last_set = int(basename[-3:])
1356 1361
1357 1362 self.online = online
1358 1363 self.realtime = realtime
1359 1364 self.delay = delay
1360 1365 ext = ext.lower()
1361 1366 self.ext = ext
1362 1367 self.getByBlock = getblock
1363 1368 self.nTxs = nTxs
1364 1369 self.startTime = startTime
1365 1370 self.endTime = endTime
1366 1371 self.endDate = endDate
1367 1372 self.startDate = startDate
1368 1373 # Added-----------------
1369 1374 self.selBlocksize = blocksize
1370 1375 self.selBlocktime = blocktime
1371 1376
1372 1377 # Verbose-----------
1373 1378 self.verbose = verbose
1374 1379 self.warnings = warnings
1375 1380
1376 1381 if not(self.setNextFile()):
1377 1382 if (startDate != None) and (endDate != None):
1378 1383 print("[Reading] No files in range: %s - %s" % (datetime.datetime.combine(startDate, startTime).ctime(), datetime.datetime.combine(endDate, endTime).ctime()))
1379 1384 elif startDate != None:
1380 1385 print("[Reading] No files in range: %s" % (datetime.datetime.combine(startDate, startTime).ctime()))
1381 1386 else:
1382 1387 print("[Reading] No files")
1383 1388
1384 1389 self.fileIndex = -1
1385 1390 self.pathList = []
1386 1391 self.filenameList = []
1387 1392 return
1388 1393
1389 1394 # self.getBasicHeader()
1390 1395
1391 1396 if last_set != None:
1392 1397 self.dataOut.last_block = last_set * \
1393 1398 self.processingHeaderObj.dataBlocksPerFile + self.basicHeaderObj.dataBlock
1394 1399 return
1395 1400
1396 1401 def getBasicHeader(self):
1397 1402
1398 1403 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
1399 1404 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
1400 1405
1401 1406 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
1402 1407
1403 1408 self.dataOut.timeZone = self.basicHeaderObj.timeZone
1404 1409
1405 1410 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
1406 1411
1407 1412 self.dataOut.errorCount = self.basicHeaderObj.errorCount
1408 1413
1409 1414 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1410 1415
1411 1416 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1412 1417
1413 1418 # self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock*self.nTxs
1414 1419
1415 1420 def getFirstHeader(self):
1416 1421
1417 1422 raise NotImplementedError
1418 1423
1419 1424 def getData(self):
1420 1425
1421 1426 raise NotImplementedError
1422 1427
1423 1428 def hasNotDataInBuffer(self):
1424 1429
1425 1430 raise NotImplementedError
1426 1431
1427 1432 def readBlock(self):
1428 1433
1429 1434 raise NotImplementedError
1430 1435
1431 1436 def isEndProcess(self):
1432 1437
1433 1438 return self.flagNoMoreFiles
1434 1439
1435 1440 def printReadBlocks(self):
1436 1441
1437 1442 print("[Reading] Number of read blocks per file %04d" % self.nReadBlocks)
1438 1443
1439 1444 def printTotalBlocks(self):
1440 1445
1441 1446 print("[Reading] Number of read blocks %04d" % self.nTotalBlocks)
1442 1447
1443 1448 def printNumberOfBlock(self):
1444 1449 'SPAM!'
1445 1450
1446 1451 # if self.flagIsNewBlock:
1447 1452 # print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
1448 1453 # self.processingHeaderObj.dataBlocksPerFile,
1449 1454 # self.dataOut.datatime.ctime())
1450 1455
1451 1456 def printInfo(self):
1452 1457
1453 1458 if self.__printInfo == False:
1454 1459 return
1455 1460
1456 1461 self.basicHeaderObj.printInfo()
1457 1462 self.systemHeaderObj.printInfo()
1458 1463 self.radarControllerHeaderObj.printInfo()
1459 1464 self.processingHeaderObj.printInfo()
1460 1465
1461 1466 self.__printInfo = False
1462 1467
1463 1468 def run(self,
1464 1469 path=None,
1465 1470 startDate=None,
1466 1471 endDate=None,
1467 1472 startTime=datetime.time(0, 0, 0),
1468 1473 endTime=datetime.time(23, 59, 59),
1469 1474 set=None,
1470 1475 expLabel="",
1471 1476 ext=None,
1472 1477 online=False,
1473 1478 delay=60,
1474 1479 walk=True,
1475 1480 getblock=False,
1476 1481 nTxs=1,
1477 1482 realtime=False,
1478 1483 blocksize=None,
1479 1484 blocktime=None,
1480 1485 skip=None,
1481 1486 cursor=None,
1482 1487 warnings=True,
1483 1488 server=None,
1484 1489 verbose=True,
1485 1490 format=None,
1486 1491 oneDDict=None,
1487 1492 twoDDict=None,
1488 1493 independentParam=None, **kwargs):
1489 1494
1490 1495 if not(self.isConfig):
1491 1496 self.setup(path=path,
1492 1497 startDate=startDate,
1493 1498 endDate=endDate,
1494 1499 startTime=startTime,
1495 1500 endTime=endTime,
1496 1501 set=set,
1497 1502 expLabel=expLabel,
1498 1503 ext=ext,
1499 1504 online=online,
1500 1505 delay=delay,
1501 1506 walk=walk,
1502 1507 getblock=getblock,
1503 1508 nTxs=nTxs,
1504 1509 realtime=realtime,
1505 1510 blocksize=blocksize,
1506 1511 blocktime=blocktime,
1507 1512 skip=skip,
1508 1513 cursor=cursor,
1509 1514 warnings=warnings,
1510 1515 server=server,
1511 1516 verbose=verbose,
1512 1517 format=format,
1513 1518 oneDDict=oneDDict,
1514 1519 twoDDict=twoDDict,
1515 1520 independentParam=independentParam)
1516 1521 self.isConfig = True
1517 1522 if server is None:
1518 1523 self.getData()
1519 1524 else:
1520 1525 self.getFromServer()
1521 1526
1522 1527
1523 1528 class JRODataWriter(JRODataIO):
1524 1529
1525 1530 """
1526 1531 Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura
1527 1532 de los datos siempre se realiza por bloques.
1528 1533 """
1529 1534
1530 1535 blockIndex = 0
1531 1536
1532 1537 path = None
1533 1538
1534 1539 setFile = None
1535 1540
1536 1541 profilesPerBlock = None
1537 1542
1538 1543 blocksPerFile = None
1539 1544
1540 1545 nWriteBlocks = 0
1541 1546
1542 1547 fileDate = None
1543 1548
1544 1549 def __init__(self, dataOut=None):
1545 1550 raise NotImplementedError
1546 1551
1547 1552 def hasAllDataInBuffer(self):
1548 1553 raise NotImplementedError
1549 1554
1550 1555 def setBlockDimension(self):
1551 1556 raise NotImplementedError
1552 1557
1553 1558 def writeBlock(self):
1554 1559 raise NotImplementedError
1555 1560
1556 1561 def putData(self):
1557 1562 raise NotImplementedError
1558 1563
1559 1564 def getProcessFlags(self):
1560 1565
1561 1566 processFlags = 0
1562 1567
1563 1568 dtype_index = get_dtype_index(self.dtype)
1564 1569 procflag_dtype = get_procflag_dtype(dtype_index)
1565 1570
1566 1571 processFlags += procflag_dtype
1567 1572
1568 1573 if self.dataOut.flagDecodeData:
1569 1574 processFlags += PROCFLAG.DECODE_DATA
1570 1575
1571 1576 if self.dataOut.flagDeflipData:
1572 1577 processFlags += PROCFLAG.DEFLIP_DATA
1573 1578
1574 1579 if self.dataOut.code is not None:
1575 1580 processFlags += PROCFLAG.DEFINE_PROCESS_CODE
1576 1581
1577 1582 if self.dataOut.nCohInt > 1:
1578 1583 processFlags += PROCFLAG.COHERENT_INTEGRATION
1579 1584
1580 1585 if self.dataOut.type == "Spectra":
1581 1586 if self.dataOut.nIncohInt > 1:
1582 1587 processFlags += PROCFLAG.INCOHERENT_INTEGRATION
1583 1588
1584 1589 if self.dataOut.data_dc is not None:
1585 1590 processFlags += PROCFLAG.SAVE_CHANNELS_DC
1586 1591
1587 1592 if self.dataOut.flagShiftFFT:
1588 1593 processFlags += PROCFLAG.SHIFT_FFT_DATA
1589 1594
1590 1595 return processFlags
1591 1596
1592 1597 def setBasicHeader(self):
1593 1598
1594 1599 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1595 1600 self.basicHeaderObj.version = self.versionFile
1596 1601 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1597 1602
1598 1603 utc = numpy.floor(self.dataOut.utctime)
1599 1604 milisecond = (self.dataOut.utctime - utc) * 1000.0
1600 1605
1601 1606 self.basicHeaderObj.utc = utc
1602 1607 self.basicHeaderObj.miliSecond = milisecond
1603 1608 self.basicHeaderObj.timeZone = self.dataOut.timeZone
1604 1609 self.basicHeaderObj.dstFlag = self.dataOut.dstFlag
1605 1610 self.basicHeaderObj.errorCount = self.dataOut.errorCount
1606 1611
1607 1612 def setFirstHeader(self):
1608 1613 """
1609 1614 Obtiene una copia del First Header
1610 1615
1611 1616 Affected:
1612 1617
1613 1618 self.basicHeaderObj
1614 1619 self.systemHeaderObj
1615 1620 self.radarControllerHeaderObj
1616 1621 self.processingHeaderObj self.
1617 1622
1618 1623 Return:
1619 1624 None
1620 1625 """
1621 1626
1622 1627 raise NotImplementedError
1623 1628
1624 1629 def __writeFirstHeader(self):
1625 1630 """
1626 1631 Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader)
1627 1632
1628 1633 Affected:
1629 1634 __dataType
1630 1635
1631 1636 Return:
1632 1637 None
1633 1638 """
1634 1639
1635 1640 # CALCULAR PARAMETROS
1636 1641
1637 1642 sizeLongHeader = self.systemHeaderObj.size + \
1638 1643 self.radarControllerHeaderObj.size + self.processingHeaderObj.size
1639 1644 self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader
1640 1645
1641 1646 self.basicHeaderObj.write(self.fp)
1642 1647 self.systemHeaderObj.write(self.fp)
1643 1648 self.radarControllerHeaderObj.write(self.fp)
1644 1649 self.processingHeaderObj.write(self.fp)
1645 1650
1646 1651 def __setNewBlock(self):
1647 1652 """
1648 1653 Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header
1649 1654
1650 1655 Return:
1651 1656 0 : si no pudo escribir nada
1652 1657 1 : Si escribio el Basic el First Header
1653 1658 """
1654 1659 if self.fp == None:
1655 1660 self.setNextFile()
1656 1661
1657 1662 if self.flagIsNewFile:
1658 1663 return 1
1659 1664
1660 1665 if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile:
1661 1666 self.basicHeaderObj.write(self.fp)
1662 1667 return 1
1663 1668
1664 1669 if not(self.setNextFile()):
1665 1670 return 0
1666 1671
1667 1672 return 1
1668 1673
1669 1674 def writeNextBlock(self):
1670 1675 """
1671 1676 Selecciona el bloque siguiente de datos y los escribe en un file
1672 1677
1673 1678 Return:
1674 1679 0 : Si no hizo pudo escribir el bloque de datos
1675 1680 1 : Si no pudo escribir el bloque de datos
1676 1681 """
1677 1682 if not(self.__setNewBlock()):
1678 1683 return 0
1679 1684
1680 1685 self.writeBlock()
1681 1686
1682 1687 print("[Writing] Block No. %d/%d" % (self.blockIndex,
1683 1688 self.processingHeaderObj.dataBlocksPerFile))
1684 1689
1685 1690 return 1
1686 1691
1687 1692 def setNextFile(self):
1688 1693 """
1689 1694 Determina el siguiente file que sera escrito
1690 1695
1691 1696 Affected:
1692 1697 self.filename
1693 1698 self.subfolder
1694 1699 self.fp
1695 1700 self.setFile
1696 1701 self.flagIsNewFile
1697 1702
1698 1703 Return:
1699 1704 0 : Si el archivo no puede ser escrito
1700 1705 1 : Si el archivo esta listo para ser escrito
1701 1706 """
1702 1707 ext = self.ext
1703 1708 path = self.path
1704 1709
1705 1710 if self.fp != None:
1706 1711 self.fp.close()
1707 1712
1708 1713 timeTuple = time.localtime(self.dataOut.utctime)
1709 1714 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
1710 1715
1711 1716 fullpath = os.path.join(path, subfolder)
1712 1717 setFile = self.setFile
1713 1718
1714 1719 if not(os.path.exists(fullpath)):
1715 1720 os.mkdir(fullpath)
1716 1721 setFile = -1 # inicializo mi contador de seteo
1717 1722 else:
1718 1723 filesList = os.listdir(fullpath)
1719 1724 if len(filesList) > 0:
1720 1725 filesList = sorted(filesList, key=str.lower)
1721 1726 filen = filesList[-1]
1722 1727 # el filename debera tener el siguiente formato
1723 1728 # 0 1234 567 89A BCDE (hex)
1724 1729 # x YYYY DDD SSS .ext
1725 1730 if isNumber(filen[8:11]):
1726 1731 # inicializo mi contador de seteo al seteo del ultimo file
1727 1732 setFile = int(filen[8:11])
1728 1733 else:
1729 1734 setFile = -1
1730 1735 else:
1731 1736 setFile = -1 # inicializo mi contador de seteo
1732 1737
1733 1738 setFile += 1
1734 1739
1735 1740 # If this is a new day it resets some values
1736 1741 if self.dataOut.datatime.date() > self.fileDate:
1737 1742 setFile = 0
1738 1743 self.nTotalBlocks = 0
1739 1744
1740 1745 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1741 1746 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1742 1747
1743 1748 filename = os.path.join(path, subfolder, filen)
1744 1749
1745 1750 fp = open(filename, 'wb')
1746 1751
1747 1752 self.blockIndex = 0
1748 1753
1749 1754 # guardando atributos
1750 1755 self.filename = filename
1751 1756 self.subfolder = subfolder
1752 1757 self.fp = fp
1753 1758 self.setFile = setFile
1754 1759 self.flagIsNewFile = 1
1755 1760 self.fileDate = self.dataOut.datatime.date()
1756 1761
1757 1762 self.setFirstHeader()
1758 1763
1759 1764 print('[Writing] Opening file: %s' % self.filename)
1760 1765
1761 1766 self.__writeFirstHeader()
1762 1767
1763 1768 return 1
1764 1769
1765 1770 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1766 1771 """
1767 1772 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1768 1773
1769 1774 Inputs:
1770 1775 path : directory where data will be saved
1771 1776 profilesPerBlock : number of profiles per block
1772 1777 set : initial file set
1773 1778 datatype : An integer number that defines data type:
1774 1779 0 : int8 (1 byte)
1775 1780 1 : int16 (2 bytes)
1776 1781 2 : int32 (4 bytes)
1777 1782 3 : int64 (8 bytes)
1778 1783 4 : float32 (4 bytes)
1779 1784 5 : double64 (8 bytes)
1780 1785
1781 1786 Return:
1782 1787 0 : Si no realizo un buen seteo
1783 1788 1 : Si realizo un buen seteo
1784 1789 """
1785 1790
1786 1791 if ext == None:
1787 1792 ext = self.ext
1788 1793
1789 1794 self.ext = ext.lower()
1790 1795
1791 1796 self.path = path
1792 1797
1793 1798 if set is None:
1794 1799 self.setFile = -1
1795 1800 else:
1796 1801 self.setFile = set - 1
1797 1802
1798 1803 self.blocksPerFile = blocksPerFile
1799 1804
1800 1805 self.profilesPerBlock = profilesPerBlock
1801 1806
1802 1807 self.dataOut = dataOut
1803 1808 self.fileDate = self.dataOut.datatime.date()
1804 1809 # By default
1805 1810 self.dtype = self.dataOut.dtype
1806 1811
1807 1812 if datatype is not None:
1808 1813 self.dtype = get_numpy_dtype(datatype)
1809 1814
1810 1815 if not(self.setNextFile()):
1811 1816 print("[Writing] There isn't a next file")
1812 1817 return 0
1813 1818
1814 1819 self.setBlockDimension()
1815 1820
1816 1821 return 1
1817 1822
1818 1823 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1819 1824
1820 1825 if not(self.isConfig):
1821 1826
1822 1827 self.setup(dataOut, path, blocksPerFile, profilesPerBlock=profilesPerBlock,
1823 1828 set=set, ext=ext, datatype=datatype, **kwargs)
1824 1829 self.isConfig = True
1825 1830
1826 1831 self.dataOut = dataOut
1827 1832 self.putData()
1828 1833 return self.dataOut No newline at end of file
@@ -1,390 +1,404
1 1 '''
2 2 Updated for multiprocessing
3 3 Author : Sergio Cortez
4 4 Jan 2018
5 5 Abstract:
6 6 Base class for processing units and operations. A decorator provides multiprocessing features and interconnect the processes created.
7 7 The argument (kwargs) sent from the controller is parsed and filtered via the decorator for each processing unit or operation instantiated.
8 8 The decorator handle also the methods inside the processing unit to be called from the main script (not as operations) (OPERATION -> type ='self').
9 9
10 10 Based on:
11 11 $Author: murco $
12 12 $Id: jroproc_base.py 1 2012-11-12 18:56:07Z murco $
13 13 '''
14 14
15 15 import os
16 16 import inspect
17 17 import zmq
18 18 import time
19 19 import pickle
20 20 from multiprocessing import Process
21 21 from zmq.utils.monitor import recv_monitor_message
22 22
23 23 from schainpy.utils import log
24 24
25 25
26 26 class ProcessingUnit(object):
27 27
28 28 """
29 29 Update - Jan 2018 - MULTIPROCESSING
30 30 All the "call" methods present in the previous base were removed.
31 31 The majority of operations are independant processes, thus
32 32 the decorator is in charge of communicate the operation processes
33 33 with the proccessing unit via IPC.
34 34
35 35 The constructor does not receive any argument. The remaining methods
36 36 are related with the operations to execute.
37 37
38 38
39 39 """
40 40
41 41 def __init__(self):
42 42
43 43 self.dataIn = None
44 44 self.dataOut = None
45 45 self.isConfig = False
46 46 self.operations = []
47 47 self.plots = []
48 48
49 49 def getAllowedArgs(self):
50 50 if hasattr(self, '__attrs__'):
51 51 return self.__attrs__
52 52 else:
53 53 return inspect.getargspec(self.run).args
54 54
55 55 def addOperation(self, conf, operation):
56 56 """
57 57 This method is used in the controller, and update the dictionary containing the operations to execute. The dict
58 58 posses the id of the operation process (IPC purposes)
59 59
60 60 Agrega un objeto del tipo "Operation" (opObj) a la lista de objetos "self.objectList" y retorna el
61 61 identificador asociado a este objeto.
62 62
63 63 Input:
64 64
65 65 object : objeto de la clase "Operation"
66 66
67 67 Return:
68 68
69 69 objId : identificador del objeto, necesario para comunicar con master(procUnit)
70 70 """
71 71
72 72 self.operations.append(
73 73 (operation, conf.type, conf.id, conf.getKwargs()))
74 74
75 75 if 'plot' in self.name.lower():
76 76 self.plots.append(operation.CODE)
77 77
78 78 def getOperationObj(self, objId):
79 79
80 80 if objId not in list(self.operations.keys()):
81 81 return None
82 82
83 83 return self.operations[objId]
84 84
85 85 def operation(self, **kwargs):
86 86 """
87 87 Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
88 88 atributos del objeto dataOut
89 89
90 90 Input:
91 91
92 92 **kwargs : Diccionario de argumentos de la funcion a ejecutar
93 93 """
94 94
95 95 raise NotImplementedError
96 96
97 97 def setup(self):
98 98
99 99 raise NotImplementedError
100 100
101 101 def run(self):
102 102
103 103 raise NotImplementedError
104 104
105 105 def close(self):
106 106
107 107 return
108 108
109 109
110 110 class Operation(object):
111 111
112 112 """
113 113 Update - Jan 2018 - MULTIPROCESSING
114 114
115 115 Most of the methods remained the same. The decorator parse the arguments and executed the run() method for each process.
116 116 The constructor doe snot receive any argument, neither the baseclass.
117 117
118 118
119 119 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
120 120 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
121 121 acumulacion dentro de esta clase
122 122
123 123 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
124 124
125 125 """
126 126
127 127 def __init__(self):
128 128
129 129 self.id = None
130 130 self.isConfig = False
131 131
132 132 if not hasattr(self, 'name'):
133 133 self.name = self.__class__.__name__
134 134
135 135 def getAllowedArgs(self):
136 136 if hasattr(self, '__attrs__'):
137 137 return self.__attrs__
138 138 else:
139 139 return inspect.getargspec(self.run).args
140 140
141 141 def setup(self):
142 142
143 143 self.isConfig = True
144 144
145 145 raise NotImplementedError
146 146
147 147 def run(self, dataIn, **kwargs):
148 148 """
149 149 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
150 150 atributos del objeto dataIn.
151 151
152 152 Input:
153 153
154 154 dataIn : objeto del tipo JROData
155 155
156 156 Return:
157 157
158 158 None
159 159
160 160 Affected:
161 161 __buffer : buffer de recepcion de datos.
162 162
163 163 """
164 164 if not self.isConfig:
165 165 self.setup(**kwargs)
166 166
167 167 raise NotImplementedError
168 168
169 169 def close(self):
170 170
171 171 return
172 172
173 173
174 174 def MPDecorator(BaseClass):
175 175 """
176 176 Multiprocessing class decorator
177 177
178 178 This function add multiprocessing features to a BaseClass. Also, it handle
179 179 the communication beetween processes (readers, procUnits and operations).
180 180 """
181 181
182 182 class MPClass(BaseClass, Process):
183 183
184 184 def __init__(self, *args, **kwargs):
185 185 super(MPClass, self).__init__()
186 186 Process.__init__(self)
187 187 self.operationKwargs = {}
188 188 self.args = args
189 189 self.kwargs = kwargs
190 190 self.sender = None
191 191 self.receiver = None
192 self.i = 0
192 193 self.name = BaseClass.__name__
193 194 if 'plot' in self.name.lower() and not self.name.endswith('_'):
194 195 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
195 196 self.start_time = time.time()
196 197
197 198 if len(self.args) is 3:
198 199 self.typeProc = "ProcUnit"
199 200 self.id = args[0]
200 201 self.inputId = args[1]
201 202 self.project_id = args[2]
202 203 elif len(self.args) is 2:
203 204 self.id = args[0]
204 205 self.inputId = args[0]
205 206 self.project_id = args[1]
206 207 self.typeProc = "Operation"
208
209 def fix_publish(self,valor,multiple1):
210 return True if valor%multiple1 ==0 else False
207 211
208 212 def subscribe(self):
209 213 '''
210 214 This function create a socket to receive objects from the
211 215 topic `inputId`.
212 216 '''
213 217
214 218 c = zmq.Context()
215 219 self.receiver = c.socket(zmq.SUB)
216 220 self.receiver.connect(
217 221 'ipc:///tmp/schain/{}_pub'.format(self.project_id))
218 222 self.receiver.setsockopt(zmq.SUBSCRIBE, self.inputId.encode())
219 223
220 224 def listen(self):
221 225 '''
222 226 This function waits for objects and deserialize using pickle
223 227 '''
224
225 data = pickle.loads(self.receiver.recv_multipart()[1])
228 try:
229 data = pickle.loads(self.receiver.recv_multipart()[1])
230 except zmq.ZMQError as e:
231 if e.errno == zmq.ETERM:
232 print (e.errno)
226 233
227 234 return data
228 235
229 236 def set_publisher(self):
230 237 '''
231 238 This function create a socket for publishing purposes.
232 239 '''
233 240
234 241 time.sleep(1)
235 242 c = zmq.Context()
236 243 self.sender = c.socket(zmq.PUB)
237 244 self.sender.connect(
238 245 'ipc:///tmp/schain/{}_sub'.format(self.project_id))
239 246
240 247 def publish(self, data, id):
241 248 '''
242 249 This function publish an object, to a specific topic.
250 The fix method only affect inputId None which is Read Unit
251 Use value between 64 80, you should notice a little retard in processing
243 252 '''
253 if self.inputId is None:
254 self.i+=1
255 if self.fix_publish(self.i,80) == True:# value n
256 time.sleep(0.01)
257
244 258 self.sender.send_multipart([str(id).encode(), pickle.dumps(data)])
245 259
246 260 def runReader(self):
247 261 '''
248 262 Run fuction for read units
249 263 '''
250 264 while True:
251 265
252 266 BaseClass.run(self, **self.kwargs)
253 267
254 268 for op, optype, opId, kwargs in self.operations:
255 269 if optype == 'self' and not self.dataOut.flagNoData:
256 270 op(**kwargs)
257 271 elif optype == 'other' and not self.dataOut.flagNoData:
258 272 self.dataOut = op.run(self.dataOut, **self.kwargs)
259 273 elif optype == 'external':
260 274 self.publish(self.dataOut, opId)
261 275
262 276 if self.dataOut.flagNoData and not self.dataOut.error:
263 277 continue
264 278
265 279 self.publish(self.dataOut, self.id)
266 280
267 281 if self.dataOut.error:
268 282 log.error(self.dataOut.error, self.name)
269 283 # self.sender.send_multipart([str(self.project_id).encode(), 'end'.encode()])
270 284 break
271 285
272 286 time.sleep(1)
273 287
274 288 def runProc(self):
275 289 '''
276 290 Run function for proccessing units
277 291 '''
278 292
279 293 while True:
280 294 self.dataIn = self.listen()
281 295
282 296 if self.dataIn.flagNoData and self.dataIn.error is None:
283 297 continue
284 298
285 299 BaseClass.run(self, **self.kwargs)
286 300
287 301 if self.dataIn.error:
288 302 self.dataOut.error = self.dataIn.error
289 303 self.dataOut.flagNoData = True
290 304
291 305 for op, optype, opId, kwargs in self.operations:
292 306 if optype == 'self' and not self.dataOut.flagNoData:
293 307 op(**kwargs)
294 308 elif optype == 'other' and not self.dataOut.flagNoData:
295 309 self.dataOut = op.run(self.dataOut, **kwargs)
296 310 elif optype == 'external' and not self.dataOut.flagNoData:
297 311 self.publish(self.dataOut, opId)
298 312
299 313 if not self.dataOut.flagNoData or self.dataOut.error:
300 314 self.publish(self.dataOut, self.id)
301 315 for op, optype, opId, kwargs in self.operations:
302 316 if optype == 'self' and self.dataOut.error:
303 317 op(**kwargs)
304 318 elif optype == 'other' and self.dataOut.error:
305 319 self.dataOut = op.run(self.dataOut, **kwargs)
306 320 elif optype == 'external' and self.dataOut.error:
307 321 self.publish(self.dataOut, opId)
308 322
309 323 if self.dataIn.error:
310 324 break
311 325
312 326 time.sleep(1)
313 327
314 328 def runOp(self):
315 329 '''
316 330 Run function for external operations (this operations just receive data
317 331 ex: plots, writers, publishers)
318 332 '''
319 333
320 334 while True:
321 335
322 336 dataOut = self.listen()
323 337
324 338 BaseClass.run(self, dataOut, **self.kwargs)
325 339
326 340 if dataOut.error:
327 341 break
328 342
329 343 time.sleep(1)
330 344
331 345 def run(self):
332 346 if self.typeProc is "ProcUnit":
333 347
334 348 if self.inputId is not None:
335 349
336 350 self.subscribe()
337 351
338 352 self.set_publisher()
339 353
340 354 if 'Reader' not in BaseClass.__name__:
341 355 self.runProc()
342 356 else:
343 357 self.runReader()
344 358
345 359 elif self.typeProc is "Operation":
346 360
347 361 self.subscribe()
348 362 self.runOp()
349 363
350 364 else:
351 365 raise ValueError("Unknown type")
352 366
353 367 self.close()
354 368
355 369 def event_monitor(self, monitor):
356 370
357 371 events = {}
358 372
359 373 for name in dir(zmq):
360 374 if name.startswith('EVENT_'):
361 375 value = getattr(zmq, name)
362 376 events[value] = name
363 377
364 378 while monitor.poll():
365 379 evt = recv_monitor_message(monitor)
366 380 if evt['event'] == 32:
367 381 self.connections += 1
368 382 if evt['event'] == 512:
369 383 pass
370 384
371 385 evt.update({'description': events[evt['event']]})
372 386
373 387 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
374 388 break
375 389 monitor.close()
376 390 print('event monitor thread done!')
377 391
378 392 def close(self):
379 393
380 394 BaseClass.close(self)
381 395
382 396 if self.sender:
383 397 self.sender.close()
384 398
385 399 if self.receiver:
386 400 self.receiver.close()
387 401
388 402 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
389 403
390 404 return MPClass
General Comments 0
You need to be logged in to leave comments. Login now