##// END OF EJS Templates
Fix ParametersPlot & BLTRParamReader
Juan C. Espinoza -
r1322:f141c0751dba
parent child
Show More
@@ -1,1397 +1,1397
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 import schainpy.admin
13 13 from schainpy.utils import log
14 14 from .jroheaderIO import SystemHeader, RadarControllerHeader
15 15 from schainpy.model.data import _noise
16 16
17 17
18 18 def getNumpyDtype(dataTypeCode):
19 19
20 20 if dataTypeCode == 0:
21 21 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
22 22 elif dataTypeCode == 1:
23 23 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
24 24 elif dataTypeCode == 2:
25 25 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
26 26 elif dataTypeCode == 3:
27 27 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
28 28 elif dataTypeCode == 4:
29 29 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
30 30 elif dataTypeCode == 5:
31 31 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
32 32 else:
33 33 raise ValueError('dataTypeCode was not defined')
34 34
35 35 return numpyDtype
36 36
37 37
38 38 def getDataTypeCode(numpyDtype):
39 39
40 40 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
41 41 datatype = 0
42 42 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
43 43 datatype = 1
44 44 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
45 45 datatype = 2
46 46 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
47 47 datatype = 3
48 48 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
49 49 datatype = 4
50 50 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
51 51 datatype = 5
52 52 else:
53 53 datatype = None
54 54
55 55 return datatype
56 56
57 57
58 58 def hildebrand_sekhon(data, navg):
59 59 """
60 60 This method is for the objective determination of the noise level in Doppler spectra. This
61 61 implementation technique is based on the fact that the standard deviation of the spectral
62 62 densities is equal to the mean spectral density for white Gaussian noise
63 63
64 64 Inputs:
65 65 Data : heights
66 66 navg : numbers of averages
67 67
68 68 Return:
69 69 mean : noise's level
70 70 """
71 71
72 72 sortdata = numpy.sort(data, axis=None)
73 73 '''
74 74 lenOfData = len(sortdata)
75 75 nums_min = lenOfData*0.2
76 76
77 77 if nums_min <= 5:
78 78
79 79 nums_min = 5
80 80
81 81 sump = 0.
82 82 sumq = 0.
83 83
84 84 j = 0
85 85 cont = 1
86 86
87 87 while((cont == 1)and(j < lenOfData)):
88 88
89 89 sump += sortdata[j]
90 90 sumq += sortdata[j]**2
91 91
92 92 if j > nums_min:
93 93 rtest = float(j)/(j-1) + 1.0/navg
94 94 if ((sumq*j) > (rtest*sump**2)):
95 95 j = j - 1
96 96 sump = sump - sortdata[j]
97 97 sumq = sumq - sortdata[j]**2
98 98 cont = 0
99 99
100 100 j += 1
101 101
102 102 lnoise = sump / j
103 103 '''
104 104 return _noise.hildebrand_sekhon(sortdata, navg)
105 105
106 106
107 107 class Beam:
108 108
109 109 def __init__(self):
110 110 self.codeList = []
111 111 self.azimuthList = []
112 112 self.zenithList = []
113 113
114 114
115 115 class GenericData(object):
116 116
117 117 flagNoData = True
118 118
119 119 def copy(self, inputObj=None):
120 120
121 121 if inputObj == None:
122 122 return copy.deepcopy(self)
123 123
124 124 for key in list(inputObj.__dict__.keys()):
125 125
126 126 attribute = inputObj.__dict__[key]
127 127
128 128 # If this attribute is a tuple or list
129 129 if type(inputObj.__dict__[key]) in (tuple, list):
130 130 self.__dict__[key] = attribute[:]
131 131 continue
132 132
133 133 # If this attribute is another object or instance
134 134 if hasattr(attribute, '__dict__'):
135 135 self.__dict__[key] = attribute.copy()
136 136 continue
137 137
138 138 self.__dict__[key] = inputObj.__dict__[key]
139 139
140 140 def deepcopy(self):
141 141
142 142 return copy.deepcopy(self)
143 143
144 144 def isEmpty(self):
145 145
146 146 return self.flagNoData
147 147
148 148 def isReady(self):
149 149
150 150 return not self.flagNoData
151 151
152 152
153 153 class JROData(GenericData):
154 154
155 155 # m_BasicHeader = BasicHeader()
156 156 # m_ProcessingHeader = ProcessingHeader()
157 157
158 158 systemHeaderObj = SystemHeader()
159 159 radarControllerHeaderObj = RadarControllerHeader()
160 160 # data = None
161 161 type = None
162 162 datatype = None # dtype but in string
163 163 # dtype = None
164 164 # nChannels = None
165 165 # nHeights = None
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 # nCode = None
177 177 # nBaud = None
178 178 # code = None
179 179 flagDecodeData = False # asumo q la data no esta decodificada
180 180 flagDeflipData = False # asumo q la data no esta sin flip
181 181 flagShiftFFT = False
182 182 # ippSeconds = None
183 183 # timeInterval = None
184 184 nCohInt = None
185 185 # noise = None
186 186 windowOfFilter = 1
187 187 # Speed of ligth
188 188 C = 3e8
189 189 frequency = 49.92e6
190 190 realtime = False
191 191 beacon_heiIndexList = None
192 192 last_block = None
193 193 blocknow = None
194 194 azimuth = None
195 195 zenith = None
196 196 beam = Beam()
197 197 profileIndex = None
198 198 error = None
199 199 data = None
200 200 nmodes = None
201 201
202 202 def __str__(self):
203 203
204 204 return '{} - {}'.format(self.type, self.getDatatime())
205 205
206 206 def getNoise(self):
207 207
208 208 raise NotImplementedError
209 209
210 210 def getNChannels(self):
211 211
212 212 return len(self.channelList)
213 213
214 214 def getChannelIndexList(self):
215 215
216 216 return list(range(self.nChannels))
217 217
218 218 def getNHeights(self):
219 219
220 220 return len(self.heightList)
221 221
222 222 def getHeiRange(self, extrapoints=0):
223 223
224 224 heis = self.heightList
225 225 # deltah = self.heightList[1] - self.heightList[0]
226 226 #
227 227 # heis.append(self.heightList[-1])
228 228
229 229 return heis
230 230
231 231 def getDeltaH(self):
232 232
233 233 delta = self.heightList[1] - self.heightList[0]
234 234
235 235 return delta
236 236
237 237 def getltctime(self):
238 238
239 239 if self.useLocalTime:
240 240 return self.utctime - self.timeZone * 60
241 241
242 242 return self.utctime
243 243
244 244 def getDatatime(self):
245 245
246 246 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
247 247 return datatimeValue
248 248
249 249 def getTimeRange(self):
250 250
251 251 datatime = []
252 252
253 253 datatime.append(self.ltctime)
254 254 datatime.append(self.ltctime + self.timeInterval + 1)
255 255
256 256 datatime = numpy.array(datatime)
257 257
258 258 return datatime
259 259
260 260 def getFmaxTimeResponse(self):
261 261
262 262 period = (10**-6) * self.getDeltaH() / (0.15)
263 263
264 264 PRF = 1. / (period * self.nCohInt)
265 265
266 266 fmax = PRF
267 267
268 268 return fmax
269 269
270 270 def getFmax(self):
271 271 PRF = 1. / (self.ippSeconds * self.nCohInt)
272 272
273 273 fmax = PRF
274 274 return fmax
275 275
276 276 def getVmax(self):
277 277
278 278 _lambda = self.C / self.frequency
279 279
280 280 vmax = self.getFmax() * _lambda / 2
281 281
282 282 return vmax
283 283
284 284 def get_ippSeconds(self):
285 285 '''
286 286 '''
287 287 return self.radarControllerHeaderObj.ippSeconds
288 288
289 289 def set_ippSeconds(self, ippSeconds):
290 290 '''
291 291 '''
292 292
293 293 self.radarControllerHeaderObj.ippSeconds = ippSeconds
294 294
295 295 return
296 296
297 297 def get_dtype(self):
298 298 '''
299 299 '''
300 300 return getNumpyDtype(self.datatype)
301 301
302 302 def set_dtype(self, numpyDtype):
303 303 '''
304 304 '''
305 305
306 306 self.datatype = getDataTypeCode(numpyDtype)
307 307
308 308 def get_code(self):
309 309 '''
310 310 '''
311 311 return self.radarControllerHeaderObj.code
312 312
313 313 def set_code(self, code):
314 314 '''
315 315 '''
316 316 self.radarControllerHeaderObj.code = code
317 317
318 318 return
319 319
320 320 def get_ncode(self):
321 321 '''
322 322 '''
323 323 return self.radarControllerHeaderObj.nCode
324 324
325 325 def set_ncode(self, nCode):
326 326 '''
327 327 '''
328 328 self.radarControllerHeaderObj.nCode = nCode
329 329
330 330 return
331 331
332 332 def get_nbaud(self):
333 333 '''
334 334 '''
335 335 return self.radarControllerHeaderObj.nBaud
336 336
337 337 def set_nbaud(self, nBaud):
338 338 '''
339 339 '''
340 340 self.radarControllerHeaderObj.nBaud = nBaud
341 341
342 342 return
343 343
344 344 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
345 345 channelIndexList = property(
346 346 getChannelIndexList, "I'm the 'channelIndexList' property.")
347 347 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
348 348 #noise = property(getNoise, "I'm the 'nHeights' property.")
349 349 datatime = property(getDatatime, "I'm the 'datatime' property")
350 350 ltctime = property(getltctime, "I'm the 'ltctime' property")
351 351 ippSeconds = property(get_ippSeconds, set_ippSeconds)
352 352 dtype = property(get_dtype, set_dtype)
353 353 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
354 354 code = property(get_code, set_code)
355 355 nCode = property(get_ncode, set_ncode)
356 356 nBaud = property(get_nbaud, set_nbaud)
357 357
358 358
359 359 class Voltage(JROData):
360 360
361 361 # data es un numpy array de 2 dmensiones (canales, alturas)
362 362 data = None
363 363 data_intensity = None
364 364 data_velocity = None
365 365 data_specwidth = None
366 366 def __init__(self):
367 367 '''
368 368 Constructor
369 369 '''
370 370
371 371 self.useLocalTime = True
372 372 self.radarControllerHeaderObj = RadarControllerHeader()
373 373 self.systemHeaderObj = SystemHeader()
374 374 self.type = "Voltage"
375 375 self.data = None
376 376 # self.dtype = None
377 377 # self.nChannels = 0
378 378 # self.nHeights = 0
379 379 self.nProfiles = None
380 380 self.heightList = None
381 381 self.channelList = None
382 382 # self.channelIndexList = None
383 383 self.flagNoData = True
384 384 self.flagDiscontinuousBlock = False
385 385 self.utctime = None
386 386 self.timeZone = None
387 387 self.dstFlag = None
388 388 self.errorCount = None
389 389 self.nCohInt = None
390 390 self.blocksize = None
391 391 self.flagCohInt = False
392 392 self.flagDecodeData = False # asumo q la data no esta decodificada
393 393 self.flagDeflipData = False # asumo q la data no esta sin flip
394 394 self.flagShiftFFT = False
395 395 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
396 396 self.profileIndex = 0
397 397
398 398 def getNoisebyHildebrand(self, channel=None):
399 399 """
400 400 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
401 401
402 402 Return:
403 403 noiselevel
404 404 """
405 405
406 406 if channel != None:
407 407 data = self.data[channel]
408 408 nChannels = 1
409 409 else:
410 410 data = self.data
411 411 nChannels = self.nChannels
412 412
413 413 noise = numpy.zeros(nChannels)
414 414 power = data * numpy.conjugate(data)
415 415
416 416 for thisChannel in range(nChannels):
417 417 if nChannels == 1:
418 418 daux = power[:].real
419 419 else:
420 420 daux = power[thisChannel, :].real
421 421 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
422 422
423 423 return noise
424 424
425 425 def getNoise(self, type=1, channel=None):
426 426
427 427 if type == 1:
428 428 noise = self.getNoisebyHildebrand(channel)
429 429
430 430 return noise
431 431
432 432 def getPower(self, channel=None):
433 433
434 434 if channel != None:
435 435 data = self.data[channel]
436 436 else:
437 437 data = self.data
438 438
439 439 power = data * numpy.conjugate(data)
440 440 powerdB = 10 * numpy.log10(power.real)
441 441 powerdB = numpy.squeeze(powerdB)
442 442
443 443 return powerdB
444 444
445 445 def getTimeInterval(self):
446 446
447 447 timeInterval = self.ippSeconds * self.nCohInt
448 448
449 449 return timeInterval
450 450
451 451 noise = property(getNoise, "I'm the 'nHeights' property.")
452 452 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
453 453
454 454
455 455 class Spectra(JROData):
456 456
457 457 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
458 458 data_spc = None
459 459 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
460 460 data_cspc = None
461 461 # data dc es un numpy array de 2 dmensiones (canales, alturas)
462 462 data_dc = None
463 463 # data power
464 464 data_pwr = None
465 465 nFFTPoints = None
466 466 # nPairs = None
467 467 pairsList = None
468 468 nIncohInt = None
469 469 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
470 470 nCohInt = None # se requiere para determinar el valor de timeInterval
471 471 ippFactor = None
472 472 profileIndex = 0
473 473 plotting = "spectra"
474 474
475 475 def __init__(self):
476 476 '''
477 477 Constructor
478 478 '''
479 479
480 480 self.useLocalTime = True
481 481 self.radarControllerHeaderObj = RadarControllerHeader()
482 482 self.systemHeaderObj = SystemHeader()
483 483 self.type = "Spectra"
484 484 # self.data = None
485 485 # self.dtype = None
486 486 # self.nChannels = 0
487 487 # self.nHeights = 0
488 488 self.nProfiles = None
489 489 self.heightList = None
490 490 self.channelList = None
491 491 # self.channelIndexList = None
492 492 self.pairsList = None
493 493 self.flagNoData = True
494 494 self.flagDiscontinuousBlock = False
495 495 self.utctime = None
496 496 self.nCohInt = None
497 497 self.nIncohInt = None
498 498 self.blocksize = None
499 499 self.nFFTPoints = None
500 500 self.wavelength = None
501 501 self.flagDecodeData = False # asumo q la data no esta decodificada
502 502 self.flagDeflipData = False # asumo q la data no esta sin flip
503 503 self.flagShiftFFT = False
504 504 self.ippFactor = 1
505 505 #self.noise = None
506 506 self.beacon_heiIndexList = []
507 507 self.noise_estimation = None
508 508
509 509 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
510 510 """
511 511 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
512 512
513 513 Return:
514 514 noiselevel
515 515 """
516 516
517 517 noise = numpy.zeros(self.nChannels)
518 518
519 519 for channel in range(self.nChannels):
520 520 daux = self.data_spc[channel,
521 521 xmin_index:xmax_index, ymin_index:ymax_index]
522 522 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
523 523
524 524 return noise
525 525
526 526 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
527 527
528 528 if self.noise_estimation is not None:
529 529 # this was estimated by getNoise Operation defined in jroproc_spectra.py
530 530 return self.noise_estimation
531 531 else:
532 532 noise = self.getNoisebyHildebrand(
533 533 xmin_index, xmax_index, ymin_index, ymax_index)
534 534 return noise
535 535
536 536 def getFreqRangeTimeResponse(self, extrapoints=0):
537 537
538 538 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
539 539 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
540 540
541 541 return freqrange
542 542
543 543 def getAcfRange(self, extrapoints=0):
544 544
545 545 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
546 546 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
547 547
548 548 return freqrange
549 549
550 550 def getFreqRange(self, extrapoints=0):
551 551
552 552 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
553 553 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
554 554
555 555 return freqrange
556 556
557 557 def getVelRange(self, extrapoints=0):
558 558
559 559 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
560 560 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
561 561
562 562 if self.nmodes:
563 563 return velrange/self.nmodes
564 564 else:
565 565 return velrange
566 566
567 567 def getNPairs(self):
568 568
569 569 return len(self.pairsList)
570 570
571 571 def getPairsIndexList(self):
572 572
573 573 return list(range(self.nPairs))
574 574
575 575 def getNormFactor(self):
576 576
577 577 pwcode = 1
578 578
579 579 if self.flagDecodeData:
580 580 pwcode = numpy.sum(self.code[0]**2)
581 581 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
582 582 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
583 583
584 584 return normFactor
585 585
586 586 def getFlagCspc(self):
587 587
588 588 if self.data_cspc is None:
589 589 return True
590 590
591 591 return False
592 592
593 593 def getFlagDc(self):
594 594
595 595 if self.data_dc is None:
596 596 return True
597 597
598 598 return False
599 599
600 600 def getTimeInterval(self):
601 601
602 602 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
603 603 if self.nmodes:
604 604 return self.nmodes*timeInterval
605 605 else:
606 606 return timeInterval
607 607
608 608 def getPower(self):
609 609
610 610 factor = self.normFactor
611 611 z = self.data_spc / factor
612 612 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
613 613 avg = numpy.average(z, axis=1)
614 614
615 615 return 10 * numpy.log10(avg)
616 616
617 617 def getCoherence(self, pairsList=None, phase=False):
618 618
619 619 z = []
620 620 if pairsList is None:
621 621 pairsIndexList = self.pairsIndexList
622 622 else:
623 623 pairsIndexList = []
624 624 for pair in pairsList:
625 625 if pair not in self.pairsList:
626 626 raise ValueError("Pair %s is not in dataOut.pairsList" % (
627 627 pair))
628 628 pairsIndexList.append(self.pairsList.index(pair))
629 629 for i in range(len(pairsIndexList)):
630 630 pair = self.pairsList[pairsIndexList[i]]
631 631 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
632 632 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
633 633 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
634 634 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
635 635 if phase:
636 636 data = numpy.arctan2(avgcoherenceComplex.imag,
637 637 avgcoherenceComplex.real) * 180 / numpy.pi
638 638 else:
639 639 data = numpy.abs(avgcoherenceComplex)
640 640
641 641 z.append(data)
642 642
643 643 return numpy.array(z)
644 644
645 645 def setValue(self, value):
646 646
647 647 print("This property should not be initialized")
648 648
649 649 return
650 650
651 651 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
652 652 pairsIndexList = property(
653 653 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
654 654 normFactor = property(getNormFactor, setValue,
655 655 "I'm the 'getNormFactor' property.")
656 656 flag_cspc = property(getFlagCspc, setValue)
657 657 flag_dc = property(getFlagDc, setValue)
658 658 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
659 659 timeInterval = property(getTimeInterval, setValue,
660 660 "I'm the 'timeInterval' property")
661 661
662 662
663 663 class SpectraHeis(Spectra):
664 664
665 665 data_spc = None
666 666 data_cspc = None
667 667 data_dc = None
668 668 nFFTPoints = None
669 669 # nPairs = None
670 670 pairsList = None
671 671 nCohInt = None
672 672 nIncohInt = None
673 673
674 674 def __init__(self):
675 675
676 676 self.radarControllerHeaderObj = RadarControllerHeader()
677 677
678 678 self.systemHeaderObj = SystemHeader()
679 679
680 680 self.type = "SpectraHeis"
681 681
682 682 # self.dtype = None
683 683
684 684 # self.nChannels = 0
685 685
686 686 # self.nHeights = 0
687 687
688 688 self.nProfiles = None
689 689
690 690 self.heightList = None
691 691
692 692 self.channelList = None
693 693
694 694 # self.channelIndexList = None
695 695
696 696 self.flagNoData = True
697 697
698 698 self.flagDiscontinuousBlock = False
699 699
700 700 # self.nPairs = 0
701 701
702 702 self.utctime = None
703 703
704 704 self.blocksize = None
705 705
706 706 self.profileIndex = 0
707 707
708 708 self.nCohInt = 1
709 709
710 710 self.nIncohInt = 1
711 711
712 712 def getNormFactor(self):
713 713 pwcode = 1
714 714 if self.flagDecodeData:
715 715 pwcode = numpy.sum(self.code[0]**2)
716 716
717 717 normFactor = self.nIncohInt * self.nCohInt * pwcode
718 718
719 719 return normFactor
720 720
721 721 def getTimeInterval(self):
722 722
723 723 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
724 724
725 725 return timeInterval
726 726
727 727 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
728 728 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
729 729
730 730
731 731 class Fits(JROData):
732 732
733 733 heightList = None
734 734 channelList = None
735 735 flagNoData = True
736 736 flagDiscontinuousBlock = False
737 737 useLocalTime = False
738 738 utctime = None
739 739 timeZone = None
740 740 # ippSeconds = None
741 741 # timeInterval = None
742 742 nCohInt = None
743 743 nIncohInt = None
744 744 noise = None
745 745 windowOfFilter = 1
746 746 # Speed of ligth
747 747 C = 3e8
748 748 frequency = 49.92e6
749 749 realtime = False
750 750
751 751 def __init__(self):
752 752
753 753 self.type = "Fits"
754 754
755 755 self.nProfiles = None
756 756
757 757 self.heightList = None
758 758
759 759 self.channelList = None
760 760
761 761 # self.channelIndexList = None
762 762
763 763 self.flagNoData = True
764 764
765 765 self.utctime = None
766 766
767 767 self.nCohInt = 1
768 768
769 769 self.nIncohInt = 1
770 770
771 771 self.useLocalTime = True
772 772
773 773 self.profileIndex = 0
774 774
775 775 # self.utctime = None
776 776 # self.timeZone = None
777 777 # self.ltctime = None
778 778 # self.timeInterval = None
779 779 # self.header = None
780 780 # self.data_header = None
781 781 # self.data = None
782 782 # self.datatime = None
783 783 # self.flagNoData = False
784 784 # self.expName = ''
785 785 # self.nChannels = None
786 786 # self.nSamples = None
787 787 # self.dataBlocksPerFile = None
788 788 # self.comments = ''
789 789 #
790 790
791 791 def getltctime(self):
792 792
793 793 if self.useLocalTime:
794 794 return self.utctime - self.timeZone * 60
795 795
796 796 return self.utctime
797 797
798 798 def getDatatime(self):
799 799
800 800 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
801 801 return datatime
802 802
803 803 def getTimeRange(self):
804 804
805 805 datatime = []
806 806
807 807 datatime.append(self.ltctime)
808 808 datatime.append(self.ltctime + self.timeInterval)
809 809
810 810 datatime = numpy.array(datatime)
811 811
812 812 return datatime
813 813
814 814 def getHeiRange(self):
815 815
816 816 heis = self.heightList
817 817
818 818 return heis
819 819
820 820 def getNHeights(self):
821 821
822 822 return len(self.heightList)
823 823
824 824 def getNChannels(self):
825 825
826 826 return len(self.channelList)
827 827
828 828 def getChannelIndexList(self):
829 829
830 830 return list(range(self.nChannels))
831 831
832 832 def getNoise(self, type=1):
833 833
834 834 #noise = numpy.zeros(self.nChannels)
835 835
836 836 if type == 1:
837 837 noise = self.getNoisebyHildebrand()
838 838
839 839 if type == 2:
840 840 noise = self.getNoisebySort()
841 841
842 842 if type == 3:
843 843 noise = self.getNoisebyWindow()
844 844
845 845 return noise
846 846
847 847 def getTimeInterval(self):
848 848
849 849 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
850 850
851 851 return timeInterval
852 852
853 853 def get_ippSeconds(self):
854 854 '''
855 855 '''
856 856 return self.ipp_sec
857 857
858 858
859 859 datatime = property(getDatatime, "I'm the 'datatime' property")
860 860 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
861 861 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
862 862 channelIndexList = property(
863 863 getChannelIndexList, "I'm the 'channelIndexList' property.")
864 864 noise = property(getNoise, "I'm the 'nHeights' property.")
865 865
866 866 ltctime = property(getltctime, "I'm the 'ltctime' property")
867 867 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
868 868 ippSeconds = property(get_ippSeconds, '')
869 869
870 870 class Correlation(JROData):
871 871
872 872 noise = None
873 873 SNR = None
874 874 #--------------------------------------------------
875 875 mode = None
876 876 split = False
877 877 data_cf = None
878 878 lags = None
879 879 lagRange = None
880 880 pairsList = None
881 881 normFactor = None
882 882 #--------------------------------------------------
883 883 # calculateVelocity = None
884 884 nLags = None
885 885 nPairs = None
886 886 nAvg = None
887 887
888 888 def __init__(self):
889 889 '''
890 890 Constructor
891 891 '''
892 892 self.radarControllerHeaderObj = RadarControllerHeader()
893 893
894 894 self.systemHeaderObj = SystemHeader()
895 895
896 896 self.type = "Correlation"
897 897
898 898 self.data = None
899 899
900 900 self.dtype = None
901 901
902 902 self.nProfiles = None
903 903
904 904 self.heightList = None
905 905
906 906 self.channelList = None
907 907
908 908 self.flagNoData = True
909 909
910 910 self.flagDiscontinuousBlock = False
911 911
912 912 self.utctime = None
913 913
914 914 self.timeZone = None
915 915
916 916 self.dstFlag = None
917 917
918 918 self.errorCount = None
919 919
920 920 self.blocksize = None
921 921
922 922 self.flagDecodeData = False # asumo q la data no esta decodificada
923 923
924 924 self.flagDeflipData = False # asumo q la data no esta sin flip
925 925
926 926 self.pairsList = None
927 927
928 928 self.nPoints = None
929 929
930 930 def getPairsList(self):
931 931
932 932 return self.pairsList
933 933
934 934 def getNoise(self, mode=2):
935 935
936 936 indR = numpy.where(self.lagR == 0)[0][0]
937 937 indT = numpy.where(self.lagT == 0)[0][0]
938 938
939 939 jspectra0 = self.data_corr[:, :, indR, :]
940 940 jspectra = copy.copy(jspectra0)
941 941
942 942 num_chan = jspectra.shape[0]
943 943 num_hei = jspectra.shape[2]
944 944
945 945 freq_dc = jspectra.shape[1] / 2
946 946 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
947 947
948 948 if ind_vel[0] < 0:
949 949 ind_vel[list(range(0, 1))] = ind_vel[list(
950 950 range(0, 1))] + self.num_prof
951 951
952 952 if mode == 1:
953 953 jspectra[:, freq_dc, :] = (
954 954 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
955 955
956 956 if mode == 2:
957 957
958 958 vel = numpy.array([-2, -1, 1, 2])
959 959 xx = numpy.zeros([4, 4])
960 960
961 961 for fil in range(4):
962 962 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
963 963
964 964 xx_inv = numpy.linalg.inv(xx)
965 965 xx_aux = xx_inv[0, :]
966 966
967 967 for ich in range(num_chan):
968 968 yy = jspectra[ich, ind_vel, :]
969 969 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
970 970
971 971 junkid = jspectra[ich, freq_dc, :] <= 0
972 972 cjunkid = sum(junkid)
973 973
974 974 if cjunkid.any():
975 975 jspectra[ich, freq_dc, junkid.nonzero()] = (
976 976 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
977 977
978 978 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
979 979
980 980 return noise
981 981
982 982 def getTimeInterval(self):
983 983
984 984 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
985 985
986 986 return timeInterval
987 987
988 988 def splitFunctions(self):
989 989
990 990 pairsList = self.pairsList
991 991 ccf_pairs = []
992 992 acf_pairs = []
993 993 ccf_ind = []
994 994 acf_ind = []
995 995 for l in range(len(pairsList)):
996 996 chan0 = pairsList[l][0]
997 997 chan1 = pairsList[l][1]
998 998
999 999 # Obteniendo pares de Autocorrelacion
1000 1000 if chan0 == chan1:
1001 1001 acf_pairs.append(chan0)
1002 1002 acf_ind.append(l)
1003 1003 else:
1004 1004 ccf_pairs.append(pairsList[l])
1005 1005 ccf_ind.append(l)
1006 1006
1007 1007 data_acf = self.data_cf[acf_ind]
1008 1008 data_ccf = self.data_cf[ccf_ind]
1009 1009
1010 1010 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1011 1011
1012 1012 def getNormFactor(self):
1013 1013 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1014 1014 acf_pairs = numpy.array(acf_pairs)
1015 1015 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1016 1016
1017 1017 for p in range(self.nPairs):
1018 1018 pair = self.pairsList[p]
1019 1019
1020 1020 ch0 = pair[0]
1021 1021 ch1 = pair[1]
1022 1022
1023 1023 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1024 1024 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1025 1025 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1026 1026
1027 1027 return normFactor
1028 1028
1029 1029 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1030 1030 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1031 1031
1032 1032
1033 1033 class Parameters(Spectra):
1034 1034
1035 1035 experimentInfo = None # Information about the experiment
1036 1036 # Information from previous data
1037 1037 inputUnit = None # Type of data to be processed
1038 1038 operation = None # Type of operation to parametrize
1039 1039 # normFactor = None #Normalization Factor
1040 1040 groupList = None # List of Pairs, Groups, etc
1041 1041 # Parameters
1042 1042 data_param = None # Parameters obtained
1043 1043 data_pre = None # Data Pre Parametrization
1044 1044 data_SNR = None # Signal to Noise Ratio
1045 1045 # heightRange = None #Heights
1046 1046 abscissaList = None # Abscissa, can be velocities, lags or time
1047 1047 # noise = None #Noise Potency
1048 1048 utctimeInit = None # Initial UTC time
1049 1049 paramInterval = None # Time interval to calculate Parameters in seconds
1050 1050 useLocalTime = True
1051 1051 # Fitting
1052 1052 data_error = None # Error of the estimation
1053 1053 constants = None
1054 1054 library = None
1055 1055 # Output signal
1056 1056 outputInterval = None # Time interval to calculate output signal in seconds
1057 1057 data_output = None # Out signal
1058 1058 nAvg = None
1059 1059 noise_estimation = None
1060 1060 GauSPC = None # Fit gaussian SPC
1061 1061
1062 1062 def __init__(self):
1063 1063 '''
1064 1064 Constructor
1065 1065 '''
1066 1066 self.radarControllerHeaderObj = RadarControllerHeader()
1067 1067
1068 1068 self.systemHeaderObj = SystemHeader()
1069 1069
1070 1070 self.type = "Parameters"
1071 1071
1072 1072 def getTimeRange1(self, interval):
1073 1073
1074 1074 datatime = []
1075 1075
1076 1076 if self.useLocalTime:
1077 1077 time1 = self.utctimeInit - self.timeZone * 60
1078 1078 else:
1079 1079 time1 = self.utctimeInit
1080 1080
1081 1081 datatime.append(time1)
1082 1082 datatime.append(time1 + interval)
1083 1083 datatime = numpy.array(datatime)
1084 1084
1085 1085 return datatime
1086 1086
1087 1087 def getTimeInterval(self):
1088 1088
1089 1089 if hasattr(self, 'timeInterval1'):
1090 1090 return self.timeInterval1
1091 1091 else:
1092 1092 return self.paramInterval
1093 1093
1094 1094 def setValue(self, value):
1095 1095
1096 1096 print("This property should not be initialized")
1097 1097
1098 1098 return
1099 1099
1100 1100 def getNoise(self):
1101 1101
1102 1102 return self.spc_noise
1103 1103
1104 1104 timeInterval = property(getTimeInterval)
1105 1105 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1106 1106
1107 1107
1108 1108 class PlotterData(object):
1109 1109 '''
1110 1110 Object to hold data to be plotted
1111 1111 '''
1112 1112
1113 1113 MAXNUMX = 200
1114 1114 MAXNUMY = 200
1115 1115
1116 1116 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1117 1117
1118 1118 self.key = code
1119 1119 self.throttle = throttle_value
1120 1120 self.exp_code = exp_code
1121 1121 self.buffering = buffering
1122 1122 self.ready = False
1123 1123 self.flagNoData = False
1124 1124 self.localtime = False
1125 1125 self.data = {}
1126 1126 self.meta = {}
1127 1127 self.__heights = []
1128 1128
1129 1129 if 'snr' in code:
1130 1130 self.plottypes = ['snr']
1131 1131 elif code == 'spc':
1132 1132 self.plottypes = ['spc', 'noise', 'rti']
1133 1133 elif code == 'cspc':
1134 1134 self.plottypes = ['cspc', 'spc', 'noise', 'rti']
1135 1135 elif code == 'rti':
1136 1136 self.plottypes = ['noise', 'rti']
1137 1137 else:
1138 1138 self.plottypes = [code]
1139 1139
1140 1140 if 'snr' not in self.plottypes and snr:
1141 1141 self.plottypes.append('snr')
1142 1142
1143 1143 for plot in self.plottypes:
1144 1144 self.data[plot] = {}
1145 1145
1146 1146
1147 1147 def __str__(self):
1148 1148 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1149 1149 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1150 1150
1151 1151 def __len__(self):
1152 1152 return len(self.data[self.key])
1153 1153
1154 1154 def __getitem__(self, key):
1155 1155
1156 1156 if key not in self.data:
1157 1157 raise KeyError(log.error('Missing key: {}'.format(key)))
1158 1158 if 'spc' in key or not self.buffering:
1159 1159 ret = self.data[key][self.tm]
1160 1160 elif 'scope' in key:
1161 1161 ret = numpy.array(self.data[key][float(self.tm)])
1162 1162 else:
1163 1163 ret = numpy.array([self.data[key][x] for x in self.times])
1164 1164 if ret.ndim > 1:
1165 1165 ret = numpy.swapaxes(ret, 0, 1)
1166 1166 return ret
1167 1167
1168 1168 def __contains__(self, key):
1169 1169 return key in self.data
1170 1170
1171 1171 def setup(self):
1172 1172 '''
1173 1173 Configure object
1174 1174 '''
1175 1175 self.type = ''
1176 1176 self.ready = False
1177 1177 del self.data
1178 1178 self.data = {}
1179 1179 self.__heights = []
1180 1180 self.__all_heights = set()
1181 1181 for plot in self.plottypes:
1182 1182 if 'snr' in plot:
1183 1183 plot = 'snr'
1184 1184 elif 'spc_moments' == plot:
1185 1185 plot = 'moments'
1186 1186 self.data[plot] = {}
1187 1187
1188 1188 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1189 1189 self.data['noise'] = {}
1190 1190 self.data['rti'] = {}
1191 1191 if 'noise' not in self.plottypes:
1192 1192 self.plottypes.append('noise')
1193 1193 if 'rti' not in self.plottypes:
1194 1194 self.plottypes.append('rti')
1195 1195
1196 1196 def shape(self, key):
1197 1197 '''
1198 1198 Get the shape of the one-element data for the given key
1199 1199 '''
1200 1200
1201 1201 if len(self.data[key]):
1202 1202 if 'spc' in key or not self.buffering:
1203 1203 return self.data[key].shape
1204 1204 return self.data[key][self.times[0]].shape
1205 1205 return (0,)
1206 1206
1207 1207 def update(self, dataOut, tm):
1208 1208 '''
1209 1209 Update data object with new dataOut
1210 1210 '''
1211 1211
1212 1212 self.profileIndex = dataOut.profileIndex
1213 1213 self.tm = tm
1214 1214 self.type = dataOut.type
1215 1215 self.parameters = getattr(dataOut, 'parameters', [])
1216 1216
1217 1217 if hasattr(dataOut, 'meta'):
1218 1218 self.meta.update(dataOut.meta)
1219 1219
1220 1220 if hasattr(dataOut, 'pairsList'):
1221 1221 self.pairs = dataOut.pairsList
1222 1222
1223 1223 self.interval = dataOut.getTimeInterval()
1224 1224 self.localtime = dataOut.useLocalTime
1225 1225 if True in ['spc' in ptype for ptype in self.plottypes]:
1226 1226 self.xrange = (dataOut.getFreqRange(1)/1000.,
1227 1227 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1228 1228 self.__heights.append(dataOut.heightList)
1229 1229 self.__all_heights.update(dataOut.heightList)
1230 1230
1231 1231 for plot in self.plottypes:
1232 1232 if plot in ('spc', 'spc_moments', 'spc_cut'):
1233 1233 z = dataOut.data_spc/dataOut.normFactor
1234 1234 buffer = 10*numpy.log10(z)
1235 1235 if plot == 'cspc':
1236 1236 buffer = (dataOut.data_spc, dataOut.data_cspc)
1237 1237 if plot == 'noise':
1238 1238 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1239 1239 if plot in ('rti', 'spcprofile'):
1240 1240 buffer = dataOut.getPower()
1241 1241 if plot == 'snr_db':
1242 1242 buffer = dataOut.data_SNR
1243 1243 if plot == 'snr':
1244 1244 buffer = 10*numpy.log10(dataOut.data_SNR)
1245 1245 if plot == 'dop':
1246 1246 buffer = dataOut.data_DOP
1247 1247 if plot == 'pow':
1248 1248 buffer = 10*numpy.log10(dataOut.data_POW)
1249 1249 if plot == 'width':
1250 1250 buffer = dataOut.data_WIDTH
1251 1251 if plot == 'coh':
1252 1252 buffer = dataOut.getCoherence()
1253 1253 if plot == 'phase':
1254 1254 buffer = dataOut.getCoherence(phase=True)
1255 1255 if plot == 'output':
1256 1256 buffer = dataOut.data_output
1257 1257 if plot == 'param':
1258 1258 buffer = dataOut.data_param
1259 1259 if plot == 'scope':
1260 1260 buffer = dataOut.data
1261 1261 self.flagDataAsBlock = dataOut.flagDataAsBlock
1262 1262 self.nProfiles = dataOut.nProfiles
1263 1263 if plot == 'pp_power':
1264 1264 buffer = dataOut.data_intensity
1265 1265 self.flagDataAsBlock = dataOut.flagDataAsBlock
1266 1266 self.nProfiles = dataOut.nProfiles
1267 1267 if plot == 'pp_velocity':
1268 1268 buffer = dataOut.data_velocity
1269 1269 self.flagDataAsBlock = dataOut.flagDataAsBlock
1270 1270 self.nProfiles = dataOut.nProfiles
1271 1271 if plot == 'pp_specwidth':
1272 1272 buffer = dataOut.data_specwidth
1273 1273 self.flagDataAsBlock = dataOut.flagDataAsBlock
1274 1274 self.nProfiles = dataOut.nProfiles
1275 1275
1276 1276 if plot == 'spc':
1277 1277 self.data['spc'][tm] = buffer
1278 1278 elif plot == 'cspc':
1279 1279 self.data['cspc'][tm] = buffer
1280 1280 elif plot == 'spc_moments':
1281 1281 self.data['spc'][tm] = buffer
1282 1282 self.data['moments'][tm] = dataOut.moments
1283 1283 else:
1284 1284 if self.buffering:
1285 1285 self.data[plot][tm] = buffer
1286 1286 else:
1287 1287 self.data[plot][tm] = buffer
1288 1288
1289 1289 if dataOut.channelList is None:
1290 1290 self.channels = range(buffer.shape[0])
1291 1291 else:
1292 1292 self.channels = dataOut.channelList
1293 1293
1294 1294 if buffer is None:
1295 1295 self.flagNoData = True
1296 1296 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1297 1297
1298 1298 def normalize_heights(self):
1299 1299 '''
1300 1300 Ensure same-dimension of the data for different heighList
1301 1301 '''
1302 1302
1303 1303 H = numpy.array(list(self.__all_heights))
1304 1304 H.sort()
1305 1305 for key in self.data:
1306 1306 shape = self.shape(key)[:-1] + H.shape
1307 1307 for tm, obj in list(self.data[key].items()):
1308 h = self.__heights[self.times.index(tm)]
1308 h = self.__heights[self.times.tolist().index(tm)]
1309 1309 if H.size == h.size:
1310 1310 continue
1311 1311 index = numpy.where(numpy.in1d(H, h))[0]
1312 1312 dummy = numpy.zeros(shape) + numpy.nan
1313 1313 if len(shape) == 2:
1314 1314 dummy[:, index] = obj
1315 1315 else:
1316 1316 dummy[index] = obj
1317 1317 self.data[key][tm] = dummy
1318 1318
1319 1319 self.__heights = [H for tm in self.times]
1320 1320
1321 1321 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1322 1322 '''
1323 1323 Convert data to json
1324 1324 '''
1325 1325
1326 1326 dy = int(self.heights.size/self.MAXNUMY) + 1
1327 1327 if self.key in ('spc', 'cspc'):
1328 1328 dx = int(self.data[self.key][tm].shape[1]/self.MAXNUMX) + 1
1329 1329 data = self.roundFloats(
1330 1330 self.data[self.key][tm][::, ::dx, ::dy].tolist())
1331 1331 else:
1332 1332 if self.key is 'noise':
1333 1333 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1334 1334 else:
1335 1335 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1336 1336
1337 1337 meta = {}
1338 1338 ret = {
1339 1339 'plot': plot_name,
1340 1340 'code': self.exp_code,
1341 1341 'time': float(tm),
1342 1342 'data': data,
1343 1343 }
1344 1344 meta['type'] = plot_type
1345 1345 meta['interval'] = float(self.interval)
1346 1346 meta['localtime'] = self.localtime
1347 1347 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1348 1348 if 'spc' in self.data or 'cspc' in self.data:
1349 1349 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1350 1350 else:
1351 1351 meta['xrange'] = []
1352 1352
1353 1353 meta.update(self.meta)
1354 1354 ret['metadata'] = meta
1355 1355 return json.dumps(ret)
1356 1356
1357 1357 @property
1358 1358 def times(self):
1359 1359 '''
1360 1360 Return the list of times of the current data
1361 1361 '''
1362 1362
1363 1363 ret = numpy.array([*self.data[self.key]])
1364 1364 if self:
1365 1365 ret.sort()
1366 1366 return ret
1367 1367
1368 1368 @property
1369 1369 def min_time(self):
1370 1370 '''
1371 1371 Return the minimun time value
1372 1372 '''
1373 1373
1374 1374 return self.times[0]
1375 1375
1376 1376 @property
1377 1377 def max_time(self):
1378 1378 '''
1379 1379 Return the maximun time value
1380 1380 '''
1381 1381
1382 1382 return self.times[-1]
1383 1383
1384 1384 @property
1385 1385 def heights(self):
1386 1386 '''
1387 1387 Return the list of heights of the current data
1388 1388 '''
1389 1389
1390 1390 return numpy.array(self.__heights[-1])
1391 1391
1392 1392 @staticmethod
1393 1393 def roundFloats(obj):
1394 1394 if isinstance(obj, list):
1395 1395 return list(map(PlotterData.roundFloats, obj))
1396 1396 elif isinstance(obj, float):
1397 1397 return round(obj, 2)
@@ -1,341 +1,346
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from schainpy.model.graphics.jroplot_base import Plot, plt
6 6 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot
7 7 from schainpy.utils import log
8 8
9 9 EARTH_RADIUS = 6.3710e3
10 10
11 11
12 12 def ll2xy(lat1, lon1, lat2, lon2):
13 13
14 14 p = 0.017453292519943295
15 15 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
16 16 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
17 17 r = 12742 * numpy.arcsin(numpy.sqrt(a))
18 18 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
19 19 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
20 20 theta = -theta + numpy.pi/2
21 21 return r*numpy.cos(theta), r*numpy.sin(theta)
22 22
23 23
24 24 def km2deg(km):
25 25 '''
26 26 Convert distance in km to degrees
27 27 '''
28 28
29 29 return numpy.rad2deg(km/EARTH_RADIUS)
30 30
31 31
32 32
33 33 class SpectralMomentsPlot(SpectraPlot):
34 34 '''
35 35 Plot for Spectral Moments
36 36 '''
37 37 CODE = 'spc_moments'
38 38 colormap = 'jet'
39 39 plot_name = 'SpectralMoments'
40 40 plot_type = 'pcolor'
41 41
42 42
43 43 class SnrPlot(RTIPlot):
44 44 '''
45 45 Plot for SNR Data
46 46 '''
47 47
48 48 CODE = 'snr'
49 49 colormap = 'jet'
50 50 plot_name = 'SNR'
51 51
52 52
53 53 class DopplerPlot(RTIPlot):
54 54 '''
55 55 Plot for DOPPLER Data (1st moment)
56 56 '''
57 57
58 58 CODE = 'dop'
59 59 colormap = 'jet'
60 60 plot_name = 'DopplerShift'
61 61
62 62
63 63 class PowerPlot(RTIPlot):
64 64 '''
65 65 Plot for Power Data (0 moment)
66 66 '''
67 67
68 68 CODE = 'pow'
69 69 colormap = 'jet'
70 70 plot_name = 'TotalPower'
71 71
72 72
73 73 class SpectralWidthPlot(RTIPlot):
74 74 '''
75 75 Plot for Spectral Width Data (2nd moment)
76 76 '''
77 77
78 78 CODE = 'width'
79 79 colormap = 'jet'
80 80 plot_name = 'SpectralWidth'
81 81
82 82
83 83 class SkyMapPlot(Plot):
84 84 '''
85 85 Plot for meteors detection data
86 86 '''
87 87
88 88 CODE = 'param'
89 89
90 90 def setup(self):
91 91
92 92 self.ncols = 1
93 93 self.nrows = 1
94 94 self.width = 7.2
95 95 self.height = 7.2
96 96 self.nplots = 1
97 97 self.xlabel = 'Zonal Zenith Angle (deg)'
98 98 self.ylabel = 'Meridional Zenith Angle (deg)'
99 99 self.polar = True
100 100 self.ymin = -180
101 101 self.ymax = 180
102 102 self.colorbar = False
103 103
104 104 def plot(self):
105 105
106 106 arrayParameters = numpy.concatenate(self.data['param'])
107 107 error = arrayParameters[:, -1]
108 108 indValid = numpy.where(error == 0)[0]
109 109 finalMeteor = arrayParameters[indValid, :]
110 110 finalAzimuth = finalMeteor[:, 3]
111 111 finalZenith = finalMeteor[:, 4]
112 112
113 113 x = finalAzimuth * numpy.pi / 180
114 114 y = finalZenith
115 115
116 116 ax = self.axes[0]
117 117
118 118 if ax.firsttime:
119 119 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
120 120 else:
121 121 ax.plot.set_data(x, y)
122 122
123 123 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
124 124 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
125 125 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
126 126 dt2,
127 127 len(x))
128 128 self.titles[0] = title
129 129
130 130
131 131 class ParametersPlot(RTIPlot):
132 132 '''
133 133 Plot for data_param object
134 134 '''
135 135
136 136 CODE = 'param'
137 137 colormap = 'seismic'
138 138 plot_name = 'Parameters'
139 139
140 140 def setup(self):
141 141 self.xaxis = 'time'
142 142 self.ncols = 1
143 143 self.nrows = self.data.shape(self.CODE)[0]
144 144 self.nplots = self.nrows
145 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
146
147 if not self.xlabel:
148 self.xlabel = 'Time'
149
145 150 if self.showSNR:
146 151 self.nrows += 1
147 152 self.nplots += 1
148 153
149 154 self.ylabel = 'Height [km]'
150 155 if not self.titles:
151 156 self.titles = self.data.parameters \
152 157 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
153 158 if self.showSNR:
154 159 self.titles.append('SNR')
155 160
156 161 def plot(self):
157 162 self.data.normalize_heights()
158 163 self.x = self.data.times
159 164 self.y = self.data.heights
160 165 if self.showSNR:
161 166 self.z = numpy.concatenate(
162 167 (self.data[self.CODE], self.data['snr'])
163 168 )
164 169 else:
165 170 self.z = self.data[self.CODE]
166 171
167 172 self.z = numpy.ma.masked_invalid(self.z)
168 173
169 174 if self.decimation is None:
170 175 x, y, z = self.fill_gaps(self.x, self.y, self.z)
171 176 else:
172 177 x, y, z = self.fill_gaps(*self.decimate())
173 178
174 179 for n, ax in enumerate(self.axes):
175 180
176 181 self.zmax = self.zmax if self.zmax is not None else numpy.max(
177 182 self.z[n])
178 183 self.zmin = self.zmin if self.zmin is not None else numpy.min(
179 184 self.z[n])
180 185
181 186 if ax.firsttime:
182 187 if self.zlimits is not None:
183 188 self.zmin, self.zmax = self.zlimits[n]
184 189
185 190 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
186 191 vmin=self.zmin,
187 192 vmax=self.zmax,
188 193 cmap=self.cmaps[n]
189 194 )
190 195 else:
191 196 if self.zlimits is not None:
192 197 self.zmin, self.zmax = self.zlimits[n]
193 198 ax.collections.remove(ax.collections[0])
194 199 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
195 200 vmin=self.zmin,
196 201 vmax=self.zmax,
197 202 cmap=self.cmaps[n]
198 203 )
199 204
200 205
201 206 class OutputPlot(ParametersPlot):
202 207 '''
203 208 Plot data_output object
204 209 '''
205 210
206 211 CODE = 'output'
207 212 colormap = 'seismic'
208 213 plot_name = 'Output'
209 214
210 215
211 216 class PolarMapPlot(Plot):
212 217 '''
213 218 Plot for weather radar
214 219 '''
215 220
216 221 CODE = 'param'
217 222 colormap = 'seismic'
218 223
219 224 def setup(self):
220 225 self.ncols = 1
221 226 self.nrows = 1
222 227 self.width = 9
223 228 self.height = 8
224 229 self.mode = self.data.meta['mode']
225 230 if self.channels is not None:
226 231 self.nplots = len(self.channels)
227 232 self.nrows = len(self.channels)
228 233 else:
229 234 self.nplots = self.data.shape(self.CODE)[0]
230 235 self.nrows = self.nplots
231 236 self.channels = list(range(self.nplots))
232 237 if self.mode == 'E':
233 238 self.xlabel = 'Longitude'
234 239 self.ylabel = 'Latitude'
235 240 else:
236 241 self.xlabel = 'Range (km)'
237 242 self.ylabel = 'Height (km)'
238 243 self.bgcolor = 'white'
239 244 self.cb_labels = self.data.meta['units']
240 245 self.lat = self.data.meta['latitude']
241 246 self.lon = self.data.meta['longitude']
242 247 self.xmin, self.xmax = float(
243 248 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
244 249 self.ymin, self.ymax = float(
245 250 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
246 251 # self.polar = True
247 252
248 253 def plot(self):
249 254
250 255 for n, ax in enumerate(self.axes):
251 256 data = self.data['param'][self.channels[n]]
252 257
253 258 zeniths = numpy.linspace(
254 259 0, self.data.meta['max_range'], data.shape[1])
255 260 if self.mode == 'E':
256 261 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
257 262 r, theta = numpy.meshgrid(zeniths, azimuths)
258 263 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
259 264 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
260 265 x = km2deg(x) + self.lon
261 266 y = km2deg(y) + self.lat
262 267 else:
263 268 azimuths = numpy.radians(self.data.heights)
264 269 r, theta = numpy.meshgrid(zeniths, azimuths)
265 270 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
266 271 self.y = zeniths
267 272
268 273 if ax.firsttime:
269 274 if self.zlimits is not None:
270 275 self.zmin, self.zmax = self.zlimits[n]
271 276 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
272 277 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
273 278 vmin=self.zmin,
274 279 vmax=self.zmax,
275 280 cmap=self.cmaps[n])
276 281 else:
277 282 if self.zlimits is not None:
278 283 self.zmin, self.zmax = self.zlimits[n]
279 284 ax.collections.remove(ax.collections[0])
280 285 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
281 286 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
282 287 vmin=self.zmin,
283 288 vmax=self.zmax,
284 289 cmap=self.cmaps[n])
285 290
286 291 if self.mode == 'A':
287 292 continue
288 293
289 294 # plot district names
290 295 f = open('/data/workspace/schain_scripts/distrito.csv')
291 296 for line in f:
292 297 label, lon, lat = [s.strip() for s in line.split(',') if s]
293 298 lat = float(lat)
294 299 lon = float(lon)
295 300 # ax.plot(lon, lat, '.b', ms=2)
296 301 ax.text(lon, lat, label.decode('utf8'), ha='center',
297 302 va='bottom', size='8', color='black')
298 303
299 304 # plot limites
300 305 limites = []
301 306 tmp = []
302 307 for line in open('/data/workspace/schain_scripts/lima.csv'):
303 308 if '#' in line:
304 309 if tmp:
305 310 limites.append(tmp)
306 311 tmp = []
307 312 continue
308 313 values = line.strip().split(',')
309 314 tmp.append((float(values[0]), float(values[1])))
310 315 for points in limites:
311 316 ax.add_patch(
312 317 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
313 318
314 319 # plot Cuencas
315 320 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
316 321 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
317 322 values = [line.strip().split(',') for line in f]
318 323 points = [(float(s[0]), float(s[1])) for s in values]
319 324 ax.add_patch(Polygon(points, ec='b', fc='none'))
320 325
321 326 # plot grid
322 327 for r in (15, 30, 45, 60):
323 328 ax.add_artist(plt.Circle((self.lon, self.lat),
324 329 km2deg(r), color='0.6', fill=False, lw=0.2))
325 330 ax.text(
326 331 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
327 332 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
328 333 '{}km'.format(r),
329 334 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
330 335
331 336 if self.mode == 'E':
332 337 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
333 338 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
334 339 else:
335 340 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
336 341 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
337 342
338 343 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
339 344 self.titles = ['{} {}'.format(
340 345 self.data.parameters[x], title) for x in self.channels]
341 346
@@ -1,407 +1,355
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 import schainpy.admin
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator
18 18 from schainpy.model.data.jrodata import Parameters
19 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
19 from schainpy.model.io.jroIO_base import Reader
20 20 from schainpy.utils import log
21 21
22 22 FILE_HEADER_STRUCTURE = numpy.dtype([
23 23 ('FMN', '<u4'),
24 24 ('nrec', '<u4'),
25 25 ('fr_offset', '<u4'),
26 26 ('id', '<u4'),
27 27 ('site', 'u1', (32,))
28 28 ])
29 29
30 30 REC_HEADER_STRUCTURE = numpy.dtype([
31 31 ('rmn', '<u4'),
32 32 ('rcounter', '<u4'),
33 33 ('nr_offset', '<u4'),
34 34 ('tr_offset', '<u4'),
35 35 ('time', '<u4'),
36 36 ('time_msec', '<u4'),
37 37 ('tag', 'u1', (32,)),
38 38 ('comments', 'u1', (32,)),
39 39 ('lat', '<f4'),
40 40 ('lon', '<f4'),
41 41 ('gps_status', '<u4'),
42 42 ('freq', '<u4'),
43 43 ('freq0', '<u4'),
44 44 ('nchan', '<u4'),
45 45 ('delta_r', '<u4'),
46 46 ('nranges', '<u4'),
47 47 ('r0', '<u4'),
48 48 ('prf', '<u4'),
49 49 ('ncoh', '<u4'),
50 50 ('npoints', '<u4'),
51 51 ('polarization', '<i4'),
52 52 ('rx_filter', '<u4'),
53 53 ('nmodes', '<u4'),
54 54 ('dmode_index', '<u4'),
55 55 ('dmode_rngcorr', '<u4'),
56 56 ('nrxs', '<u4'),
57 57 ('acf_length', '<u4'),
58 58 ('acf_lags', '<u4'),
59 59 ('sea_to_atmos', '<f4'),
60 60 ('sea_notch', '<u4'),
61 61 ('lh_sea', '<u4'),
62 62 ('hh_sea', '<u4'),
63 63 ('nbins_sea', '<u4'),
64 64 ('min_snr', '<f4'),
65 65 ('min_cc', '<f4'),
66 66 ('max_time_diff', '<f4')
67 67 ])
68 68
69 69 DATA_STRUCTURE = numpy.dtype([
70 70 ('range', '<u4'),
71 71 ('status', '<u4'),
72 72 ('zonal', '<f4'),
73 73 ('meridional', '<f4'),
74 74 ('vertical', '<f4'),
75 75 ('zonal_a', '<f4'),
76 76 ('meridional_a', '<f4'),
77 77 ('corrected_fading', '<f4'), # seconds
78 78 ('uncorrected_fading', '<f4'), # seconds
79 79 ('time_diff', '<f4'),
80 80 ('major_axis', '<f4'),
81 81 ('axial_ratio', '<f4'),
82 82 ('orientation', '<f4'),
83 83 ('sea_power', '<u4'),
84 84 ('sea_algorithm', '<u4')
85 85 ])
86 86
87 87
88 class BLTRParamReader(JRODataReader, ProcessingUnit):
88 class BLTRParamReader(Reader, ProcessingUnit):
89 89 '''
90 90 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR
91 91 from *.sswma files
92 92 '''
93 93
94 94 ext = '.sswma'
95 95
96 96 def __init__(self):
97 97
98 98 ProcessingUnit.__init__(self)
99 99
100 100 self.dataOut = Parameters()
101 self.dataOut.timezone = 300
101 102 self.counter_records = 0
102 103 self.flagNoMoreFiles = 0
103 104 self.isConfig = False
104 105 self.filename = None
105
106 def setup(self,
107 path=None,
108 startDate=None,
109 endDate=None,
110 ext=None,
111 startTime=datetime.time(0, 0, 0),
112 endTime=datetime.time(23, 59, 59),
113 timezone=0,
114 status_value=0,
115 **kwargs):
116 self.path = path
117 self.startDate = startDate
118 self.endDate = endDate
119 self.startTime = startTime
120 self.endTime = endTime
121 self.status_value = status_value
106 self.status_value = 0
122 107 self.datatime = datetime.datetime(1900,1,1)
123 self.delay = kwargs.get('delay', 10)
124 self.online = kwargs.get('online', False)
125 self.nTries = kwargs.get('nTries', 3)
108 self.filefmt = "*********%Y%m%d******"
109
110 def setup(self, **kwargs):
111
112 self.set_kwargs(**kwargs)
126 113
127 114 if self.path is None:
128 115 raise ValueError("The path is not valid")
129 116
130 if ext is None:
131 ext = self.ext
132
133 self.fileList = self.search_files(self.path, startDate, endDate, ext)
134 self.timezone = timezone
135 self.fileIndex = 0
117 if self.online:
118 log.log("Searching files in online mode...", self.name)
119
120 for nTries in range(self.nTries):
121 fullpath = self.searchFilesOnLine(self.path, self.startDate,
122 self.endDate, self.expLabel, self.ext, self.walk,
123 self.filefmt, self.folderfmt)
124 try:
125 fullpath = next(fullpath)
126 except:
127 fullpath = None
128
129 if fullpath:
130 self.fileSize = os.path.getsize(fullpath)
131 self.filename = fullpath
132 self.flagIsNewFile = 1
133 if self.fp != None:
134 self.fp.close()
135 self.fp = self.open_file(fullpath, self.open_mode)
136 self.flagNoMoreFiles = 0
137 break
136 138
137 if not self.fileList:
138 raise Warning("There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' " % (
139 path))
139 log.warning(
140 'Waiting {} sec for a valid file in {}: try {} ...'.format(
141 self.delay, self.path, nTries + 1),
142 self.name)
143 time.sleep(self.delay)
140 144
141 self.setNextFile()
145 if not(fullpath):
146 raise schainpy.admin.SchainError(
147 'There isn\'t any valid file in {}'.format(self.path))
148 self.readFirstHeader()
149 else:
150 log.log("Searching files in {}".format(self.path), self.name)
151 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
152 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
153 self.setNextFile()
142 154
143 def search_last_file(self):
155 def checkForRealPath(self, nextFile, nextDay):
144 156 '''
145 Get last file and add it to the list
146 157 '''
147 158
148 for n in range(self.nTries+1):
149 if n>0:
150 log.warning(
151 "Waiting %0.2f seconds for the next file, try %03d ..." % (self.delay, n+1),
152 self.name
153 )
154 time.sleep(self.delay)
155 file_list = os.listdir(self.path)
156 file_list.sort()
157 if file_list:
158 if self.filename:
159 if file_list[-1] not in self.filename:
160 return file_list[-1]
161 else:
162 continue
163 return file_list[-1]
164 return 0
165
166 def search_files(self, path, startDate, endDate, ext):
159 dt = self.datatime + datetime.timedelta(1)
160 filename = '{}.{}{}'.format(self.siteFile, dt.strftime('%Y%m%d'), self.ext)
161 fullfilename = os.path.join(self.path, filename)
162 if os.path.exists(fullfilename):
163 return fullfilename, filename
164 return None, filename
165
166
167 def readFirstHeader(self):
167 168 '''
168 Searching for BLTR rawdata file in path
169 Creating a list of file to proces included in [startDate,endDate]
170
171 Input:
172 path - Path to find BLTR rawdata files
173 startDate - Select file from this date
174 enDate - Select file until this date
175 ext - Extension of the file to read
176 169 '''
177
178 log.success('Searching files in {} '.format(path), 'BLTRParamReader')
179 foldercounter = 0
180 fileList0 = glob.glob1(path, "*%s" % ext)
181 fileList0.sort()
182
183 for thisFile in fileList0:
184 year = thisFile[-14:-10]
185 if not isNumber(year):
186 continue
187
188 month = thisFile[-10:-8]
189 if not isNumber(month):
190 continue
191
192 day = thisFile[-8:-6]
193 if not isNumber(day):
194 continue
195 170
196 year, month, day = int(year), int(month), int(day)
197 dateFile = datetime.date(year, month, day)
198
199 if (startDate > dateFile) or (endDate < dateFile):
200 continue
201
202 yield thisFile
203
204 return
205
206 def setNextFile(self):
207
208 if self.online:
209 filename = self.search_last_file()
210 if not filename:
211 self.flagNoMoreFiles = 1
212 return 0
213 else:
214 try:
215 filename = next(self.fileList)
216 except StopIteration:
217 self.flagNoMoreFiles = 1
218 return 0
219
220 log.success('Opening {}'.format(filename), 'BLTRParamReader')
221
222 dirname, name = os.path.split(filename)
223 171 # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
224 self.siteFile = filename.split('.')[0]
225 if self.filename is not None:
226 self.fp.close()
227 self.filename = os.path.join(self.path, filename)
228 self.fp = open(self.filename, 'rb')
172 self.siteFile = self.filename.split('/')[-1].split('.')[0]
229 173 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
230 174 self.nrecords = self.header_file['nrec'][0]
231 self.sizeOfFile = os.path.getsize(self.filename)
232 175 self.counter_records = 0
233 176 self.flagIsNewFile = 0
234 177 self.fileIndex += 1
235 178
236 return 1
237
238 179 def readNextBlock(self):
239 180
240 181 while True:
241 182 if not self.online and self.counter_records == self.nrecords:
242 183 self.flagIsNewFile = 1
243 184 if not self.setNextFile():
244 185 return 0
245
246 186 try:
247 187 pointer = self.fp.tell()
248 188 self.readBlock()
249 189 except:
250 190 if self.online and self.waitDataBlock(pointer, 38512) == 1:
251 191 continue
252 192 else:
253 193 if not self.setNextFile():
254 194 return 0
255 195
256 196 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
257 197 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
258 198 log.warning(
259 199 'Reading Record No. {}/{} -> {} [Skipping]'.format(
260 200 self.counter_records,
261 201 self.nrecords,
262 202 self.datatime.ctime()),
263 203 'BLTRParamReader')
264 204 continue
265 205 break
266 206
267 207 log.log('Reading Record No. {} -> {}'.format(
268 208 self.counter_records,
269 # self.nrecords,
270 209 self.datatime.ctime()), 'BLTRParamReader')
271 210
272 211 return 1
273 212
274 213 def readBlock(self):
275 214
276 215 pointer = self.fp.tell()
277 216 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
278 217 self.nchannels = int(header_rec['nchan'][0] / 2)
279 218 self.kchan = header_rec['nrxs'][0]
280 219 self.nmodes = header_rec['nmodes'][0]
281 220 self.nranges = header_rec['nranges'][0]
282 221 self.fp.seek(pointer)
283 222 self.height = numpy.empty((self.nmodes, self.nranges))
284 223 self.snr = numpy.empty((self.nmodes, int(self.nchannels), self.nranges))
285 224 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
286 225 self.flagDiscontinuousBlock = 0
287 226
288 227 for mode in range(self.nmodes):
289 228 self.readHeader()
290 229 data = self.readData()
291 230 self.height[mode] = (data[0] - self.correction) / 1000.
292 231 self.buffer[mode] = data[1]
293 232 self.snr[mode] = data[2]
294 233
295 234 self.counter_records = self.counter_records + self.nmodes
296 235
297 236 return
298 237
299 238 def readHeader(self):
300 239 '''
301 240 RecordHeader of BLTR rawdata file
302 241 '''
303 242
304 243 header_structure = numpy.dtype(
305 244 REC_HEADER_STRUCTURE.descr + [
306 245 ('antenna_coord', 'f4', (2, int(self.nchannels))),
307 246 ('rx_gains', 'u4', (int(self.nchannels),)),
308 247 ('rx_analysis', 'u4', (int(self.nchannels),))
309 248 ]
310 249 )
311 250
312 251 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
313 252 self.lat = self.header_rec['lat'][0]
314 253 self.lon = self.header_rec['lon'][0]
315 254 self.delta = self.header_rec['delta_r'][0]
316 255 self.correction = self.header_rec['dmode_rngcorr'][0]
317 256 self.imode = self.header_rec['dmode_index'][0]
318 257 self.antenna = self.header_rec['antenna_coord']
319 258 self.rx_gains = self.header_rec['rx_gains']
320 259 self.time = self.header_rec['time'][0]
321 260 dt = datetime.datetime.utcfromtimestamp(self.time)
322 261 if dt.date()>self.datatime.date():
323 262 self.flagDiscontinuousBlock = 1
324 263 self.datatime = dt
325 264
326 265 def readData(self):
327 266 '''
328 267 Reading and filtering data block record of BLTR rawdata file,
329 268 filtering is according to status_value.
330 269
331 270 Input:
332 271 status_value - Array data is set to NAN for values that are not
333 272 equal to status_value
334 273
335 274 '''
336 275 self.nchannels = int(self.nchannels)
337 276
338 277 data_structure = numpy.dtype(
339 278 DATA_STRUCTURE.descr + [
340 279 ('rx_saturation', 'u4', (self.nchannels,)),
341 280 ('chan_offset', 'u4', (2 * self.nchannels,)),
342 281 ('rx_amp', 'u4', (self.nchannels,)),
343 282 ('rx_snr', 'f4', (self.nchannels,)),
344 283 ('cross_snr', 'f4', (self.kchan,)),
345 284 ('sea_power_relative', 'f4', (self.kchan,))]
346 285 )
347 286
348 287 data = numpy.fromfile(self.fp, data_structure, self.nranges)
349 288
350 289 height = data['range']
351 290 winds = numpy.array(
352 291 (data['zonal'], data['meridional'], data['vertical']))
353 292 snr = data['rx_snr'].T
354 293
355 294 winds[numpy.where(winds == -9999.)] = numpy.nan
356 295 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
357 296 snr[numpy.where(snr == -9999.)] = numpy.nan
358 297 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
359 298 snr = numpy.power(10, snr / 10)
360 299
361 300 return height, winds, snr
362 301
363 302 def set_output(self):
364 303 '''
365 304 Storing data from databuffer to dataOut object
366 305 '''
367 306
368 307 self.dataOut.data_SNR = self.snr
369 308 self.dataOut.height = self.height
370 309 self.dataOut.data = self.buffer
371 310 self.dataOut.utctimeInit = self.time
372 311 self.dataOut.utctime = self.dataOut.utctimeInit
373 312 self.dataOut.useLocalTime = False
374 313 self.dataOut.paramInterval = 157
375 self.dataOut.timezone = self.timezone
376 314 self.dataOut.site = self.siteFile
377 315 self.dataOut.nrecords = self.nrecords / self.nmodes
378 self.dataOut.sizeOfFile = self.sizeOfFile
379 316 self.dataOut.lat = self.lat
380 317 self.dataOut.lon = self.lon
381 318 self.dataOut.channelList = list(range(self.nchannels))
382 319 self.dataOut.kchan = self.kchan
383 320 self.dataOut.delta = self.delta
384 321 self.dataOut.correction = self.correction
385 322 self.dataOut.nmodes = self.nmodes
386 323 self.dataOut.imode = self.imode
387 324 self.dataOut.antenna = self.antenna
388 325 self.dataOut.rx_gains = self.rx_gains
389 326 self.dataOut.flagNoData = False
390 327 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
391 328
392 329 def getData(self):
393 330 '''
394 331 Storing data from databuffer to dataOut object
395 332 '''
396 333 if self.flagNoMoreFiles:
397 334 self.dataOut.flagNoData = True
398 raise schainpy.admin.SchainError('No More files to read')
335 return 0
399 336
400 337 if not self.readNextBlock():
401 338 self.dataOut.flagNoData = True
402 raise schainpy.admin.SchainError('Time for wait new file reach!!!')
339 return 0
403 340
404 341 self.set_output()
405 342
406 343 return 1
407 No newline at end of file
344
345 def run(self, **kwargs):
346 '''
347 '''
348
349 if not(self.isConfig):
350 self.setup(**kwargs)
351 self.isConfig = True
352
353 self.getData()
354
355 return No newline at end of file
@@ -1,1575 +1,1577
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 zmq
16 16
17 17 from schainpy.model.proc.jroproc_base import Operation, MPDecorator
18 18 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
19 19 from schainpy.model.data.jroheaderIO import get_dtype_index, get_numpy_dtype, get_procflag_dtype, get_dtype_width
20 20 from schainpy.utils import log
21 21 import schainpy.admin
22 22
23 23 LOCALTIME = True
24 24 DT_DIRECTIVES = {
25 25 '%Y': 4,
26 26 '%y': 2,
27 27 '%m': 2,
28 28 '%d': 2,
29 29 '%j': 3,
30 30 '%H': 2,
31 31 '%M': 2,
32 32 '%S': 2,
33 33 '%f': 6
34 34 }
35 35
36 36
37 37 def isNumber(cad):
38 38 """
39 39 Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero.
40 40
41 41 Excepciones:
42 42 Si un determinado string no puede ser convertido a numero
43 43 Input:
44 44 str, string al cual se le analiza para determinar si convertible a un numero o no
45 45
46 46 Return:
47 47 True : si el string es uno numerico
48 48 False : no es un string numerico
49 49 """
50 50 try:
51 51 float(cad)
52 52 return True
53 53 except:
54 54 return False
55 55
56 56
57 57 def isFileInEpoch(filename, startUTSeconds, endUTSeconds):
58 58 """
59 59 Esta funcion determina si un archivo de datos se encuentra o no dentro del rango de fecha especificado.
60 60
61 61 Inputs:
62 62 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
63 63
64 64 startUTSeconds : fecha inicial del rango seleccionado. La fecha esta dada en
65 65 segundos contados desde 01/01/1970.
66 66 endUTSeconds : fecha final del rango seleccionado. La fecha esta dada en
67 67 segundos contados desde 01/01/1970.
68 68
69 69 Return:
70 70 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
71 71 fecha especificado, de lo contrario retorna False.
72 72
73 73 Excepciones:
74 74 Si el archivo no existe o no puede ser abierto
75 75 Si la cabecera no puede ser leida.
76 76
77 77 """
78 78 basicHeaderObj = BasicHeader(LOCALTIME)
79 79
80 80 try:
81 81 fp = open(filename, 'rb')
82 82 except IOError:
83 83 print("The file %s can't be opened" % (filename))
84 84 return 0
85 85
86 86 sts = basicHeaderObj.read(fp)
87 87 fp.close()
88 88
89 89 if not(sts):
90 90 print("Skipping the file %s because it has not a valid header" % (filename))
91 91 return 0
92 92
93 93 if not ((startUTSeconds <= basicHeaderObj.utc) and (endUTSeconds > basicHeaderObj.utc)):
94 94 return 0
95 95
96 96 return 1
97 97
98 98
99 99 def isTimeInRange(thisTime, startTime, endTime):
100 100 if endTime >= startTime:
101 101 if (thisTime < startTime) or (thisTime > endTime):
102 102 return 0
103 103 return 1
104 104 else:
105 105 if (thisTime < startTime) and (thisTime > endTime):
106 106 return 0
107 107 return 1
108 108
109 109
110 110 def isFileInTimeRange(filename, startDate, endDate, startTime, endTime):
111 111 """
112 112 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
113 113
114 114 Inputs:
115 115 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
116 116
117 117 startDate : fecha inicial del rango seleccionado en formato datetime.date
118 118
119 119 endDate : fecha final del rango seleccionado en formato datetime.date
120 120
121 121 startTime : tiempo inicial del rango seleccionado en formato datetime.time
122 122
123 123 endTime : tiempo final del rango seleccionado en formato datetime.time
124 124
125 125 Return:
126 126 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
127 127 fecha especificado, de lo contrario retorna False.
128 128
129 129 Excepciones:
130 130 Si el archivo no existe o no puede ser abierto
131 131 Si la cabecera no puede ser leida.
132 132
133 133 """
134 134
135 135 try:
136 136 fp = open(filename, 'rb')
137 137 except IOError:
138 138 print("The file %s can't be opened" % (filename))
139 139 return None
140 140
141 141 firstBasicHeaderObj = BasicHeader(LOCALTIME)
142 142 systemHeaderObj = SystemHeader()
143 143 radarControllerHeaderObj = RadarControllerHeader()
144 144 processingHeaderObj = ProcessingHeader()
145 145
146 146 lastBasicHeaderObj = BasicHeader(LOCALTIME)
147 147
148 148 sts = firstBasicHeaderObj.read(fp)
149 149
150 150 if not(sts):
151 151 print("[Reading] Skipping the file %s because it has not a valid header" % (filename))
152 152 return None
153 153
154 154 if not systemHeaderObj.read(fp):
155 155 return None
156 156
157 157 if not radarControllerHeaderObj.read(fp):
158 158 return None
159 159
160 160 if not processingHeaderObj.read(fp):
161 161 return None
162 162
163 163 filesize = os.path.getsize(filename)
164 164
165 165 offset = processingHeaderObj.blockSize + 24 # header size
166 166
167 167 if filesize <= offset:
168 168 print("[Reading] %s: This file has not enough data" % filename)
169 169 return None
170 170
171 171 fp.seek(-offset, 2)
172 172
173 173 sts = lastBasicHeaderObj.read(fp)
174 174
175 175 fp.close()
176 176
177 177 thisDatetime = lastBasicHeaderObj.datatime
178 178 thisTime_last_block = thisDatetime.time()
179 179
180 180 thisDatetime = firstBasicHeaderObj.datatime
181 181 thisDate = thisDatetime.date()
182 182 thisTime_first_block = thisDatetime.time()
183 183
184 184 # General case
185 185 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
186 186 #-----------o----------------------------o-----------
187 187 # startTime endTime
188 188
189 189 if endTime >= startTime:
190 190 if (thisTime_last_block < startTime) or (thisTime_first_block > endTime):
191 191 return None
192 192
193 193 return thisDatetime
194 194
195 195 # If endTime < startTime then endTime belongs to the next day
196 196
197 197 #<<<<<<<<<<<o o>>>>>>>>>>>
198 198 #-----------o----------------------------o-----------
199 199 # endTime startTime
200 200
201 201 if (thisDate == startDate) and (thisTime_last_block < startTime):
202 202 return None
203 203
204 204 if (thisDate == endDate) and (thisTime_first_block > endTime):
205 205 return None
206 206
207 207 if (thisTime_last_block < startTime) and (thisTime_first_block > endTime):
208 208 return None
209 209
210 210 return thisDatetime
211 211
212 212
213 213 def isFolderInDateRange(folder, startDate=None, endDate=None):
214 214 """
215 215 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
216 216
217 217 Inputs:
218 218 folder : nombre completo del directorio.
219 219 Su formato deberia ser "/path_root/?YYYYDDD"
220 220
221 221 siendo:
222 222 YYYY : Anio (ejemplo 2015)
223 223 DDD : Dia del anio (ejemplo 305)
224 224
225 225 startDate : fecha inicial del rango seleccionado en formato datetime.date
226 226
227 227 endDate : fecha final del rango seleccionado en formato datetime.date
228 228
229 229 Return:
230 230 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
231 231 fecha especificado, de lo contrario retorna False.
232 232 Excepciones:
233 233 Si el directorio no tiene el formato adecuado
234 234 """
235 235
236 236 basename = os.path.basename(folder)
237 237
238 238 if not isRadarFolder(basename):
239 239 print("The folder %s has not the rigth format" % folder)
240 240 return 0
241 241
242 242 if startDate and endDate:
243 243 thisDate = getDateFromRadarFolder(basename)
244 244
245 245 if thisDate < startDate:
246 246 return 0
247 247
248 248 if thisDate > endDate:
249 249 return 0
250 250
251 251 return 1
252 252
253 253
254 254 def isFileInDateRange(filename, startDate=None, endDate=None):
255 255 """
256 256 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
257 257
258 258 Inputs:
259 259 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
260 260
261 261 Su formato deberia ser "?YYYYDDDsss"
262 262
263 263 siendo:
264 264 YYYY : Anio (ejemplo 2015)
265 265 DDD : Dia del anio (ejemplo 305)
266 266 sss : set
267 267
268 268 startDate : fecha inicial del rango seleccionado en formato datetime.date
269 269
270 270 endDate : fecha final del rango seleccionado en formato datetime.date
271 271
272 272 Return:
273 273 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
274 274 fecha especificado, de lo contrario retorna False.
275 275 Excepciones:
276 276 Si el archivo no tiene el formato adecuado
277 277 """
278 278
279 279 basename = os.path.basename(filename)
280 280
281 281 if not isRadarFile(basename):
282 282 print("The filename %s has not the rigth format" % filename)
283 283 return 0
284 284
285 285 if startDate and endDate:
286 286 thisDate = getDateFromRadarFile(basename)
287 287
288 288 if thisDate < startDate:
289 289 return 0
290 290
291 291 if thisDate > endDate:
292 292 return 0
293 293
294 294 return 1
295 295
296 296
297 297 def getFileFromSet(path, ext, set):
298 298 validFilelist = []
299 299 fileList = os.listdir(path)
300 300
301 301 # 0 1234 567 89A BCDE
302 302 # H YYYY DDD SSS .ext
303 303
304 304 for thisFile in fileList:
305 305 try:
306 306 year = int(thisFile[1:5])
307 307 doy = int(thisFile[5:8])
308 308 except:
309 309 continue
310 310
311 311 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
312 312 continue
313 313
314 314 validFilelist.append(thisFile)
315 315
316 316 myfile = fnmatch.filter(
317 317 validFilelist, '*%4.4d%3.3d%3.3d*' % (year, doy, set))
318 318
319 319 if len(myfile) != 0:
320 320 return myfile[0]
321 321 else:
322 322 filename = '*%4.4d%3.3d%3.3d%s' % (year, doy, set, ext.lower())
323 323 print('the filename %s does not exist' % filename)
324 324 print('...going to the last file: ')
325 325
326 326 if validFilelist:
327 327 validFilelist = sorted(validFilelist, key=str.lower)
328 328 return validFilelist[-1]
329 329
330 330 return None
331 331
332 332
333 333 def getlastFileFromPath(path, ext):
334 334 """
335 335 Depura el fileList dejando solo los que cumplan el formato de "PYYYYDDDSSS.ext"
336 336 al final de la depuracion devuelve el ultimo file de la lista que quedo.
337 337
338 338 Input:
339 339 fileList : lista conteniendo todos los files (sin path) que componen una determinada carpeta
340 340 ext : extension de los files contenidos en una carpeta
341 341
342 342 Return:
343 343 El ultimo file de una determinada carpeta, no se considera el path.
344 344 """
345 345 validFilelist = []
346 346 fileList = os.listdir(path)
347 347
348 348 # 0 1234 567 89A BCDE
349 349 # H YYYY DDD SSS .ext
350 350
351 351 for thisFile in fileList:
352 352
353 353 year = thisFile[1:5]
354 354 if not isNumber(year):
355 355 continue
356 356
357 357 doy = thisFile[5:8]
358 358 if not isNumber(doy):
359 359 continue
360 360
361 361 year = int(year)
362 362 doy = int(doy)
363 363
364 364 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
365 365 continue
366 366
367 367 validFilelist.append(thisFile)
368 368
369 369 if validFilelist:
370 370 validFilelist = sorted(validFilelist, key=str.lower)
371 371 return validFilelist[-1]
372 372
373 373 return None
374 374
375 375
376 376 def isRadarFolder(folder):
377 377 try:
378 378 year = int(folder[1:5])
379 379 doy = int(folder[5:8])
380 380 except:
381 381 return 0
382 382
383 383 return 1
384 384
385 385
386 386 def isRadarFile(file):
387 387 try:
388 388 year = int(file[1:5])
389 389 doy = int(file[5:8])
390 390 set = int(file[8:11])
391 391 except:
392 392 return 0
393 393
394 394 return 1
395 395
396 396
397 397 def getDateFromRadarFile(file):
398 398 try:
399 399 year = int(file[1:5])
400 400 doy = int(file[5:8])
401 401 set = int(file[8:11])
402 402 except:
403 403 return None
404 404
405 405 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
406 406 return thisDate
407 407
408 408
409 409 def getDateFromRadarFolder(folder):
410 410 try:
411 411 year = int(folder[1:5])
412 412 doy = int(folder[5:8])
413 413 except:
414 414 return None
415 415
416 416 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
417 417 return thisDate
418 418
419 419 def parse_format(s, fmt):
420 420
421 421 for i in range(fmt.count('%')):
422 422 x = fmt.index('%')
423 423 d = DT_DIRECTIVES[fmt[x:x+2]]
424 424 fmt = fmt.replace(fmt[x:x+2], s[x:x+d])
425 425 return fmt
426 426
427 427 class Reader(object):
428 428
429 429 c = 3E8
430 430 isConfig = False
431 431 dtype = None
432 432 pathList = []
433 433 filenameList = []
434 434 datetimeList = []
435 435 filename = None
436 436 ext = None
437 437 flagIsNewFile = 1
438 438 flagDiscontinuousBlock = 0
439 439 flagIsNewBlock = 0
440 440 flagNoMoreFiles = 0
441 441 fp = None
442 442 firstHeaderSize = 0
443 443 basicHeaderSize = 24
444 444 versionFile = 1103
445 445 fileSize = None
446 446 fileSizeByHeader = None
447 447 fileIndex = -1
448 448 profileIndex = None
449 449 blockIndex = 0
450 450 nTotalBlocks = 0
451 451 maxTimeStep = 30
452 452 lastUTTime = None
453 453 datablock = None
454 454 dataOut = None
455 455 getByBlock = False
456 456 path = None
457 457 startDate = None
458 458 endDate = None
459 459 startTime = datetime.time(0, 0, 0)
460 460 endTime = datetime.time(23, 59, 59)
461 461 set = None
462 462 expLabel = ""
463 463 online = False
464 464 delay = 60
465 465 nTries = 3 # quantity tries
466 466 nFiles = 3 # number of files for searching
467 467 walk = True
468 468 getblock = False
469 469 nTxs = 1
470 470 realtime = False
471 471 blocksize = 0
472 472 blocktime = None
473 473 warnings = True
474 474 verbose = True
475 475 server = None
476 476 format = None
477 477 oneDDict = None
478 478 twoDDict = None
479 479 independentParam = None
480 480 filefmt = None
481 481 folderfmt = None
482 482 open_file = open
483 483 open_mode = 'rb'
484 484
485 485 def run(self):
486 486
487 487 raise NotImplementedError
488 488
489 489 def getAllowedArgs(self):
490 490 if hasattr(self, '__attrs__'):
491 491 return self.__attrs__
492 492 else:
493 493 return inspect.getargspec(self.run).args
494 494
495 495 def set_kwargs(self, **kwargs):
496 496
497 497 for key, value in kwargs.items():
498 498 setattr(self, key, value)
499 499
500 500 def find_folders(self, path, startDate, endDate, folderfmt, last=False):
501 501
502 502 folders = [x for f in path.split(',')
503 503 for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))]
504 504 folders.sort()
505 505
506 506 if last:
507 507 folders = [folders[-1]]
508 508
509 509 for folder in folders:
510 510 try:
511 511 dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
512 512 if dt >= startDate and dt <= endDate:
513 513 yield os.path.join(path, folder)
514 514 else:
515 515 log.log('Skiping folder {}'.format(folder), self.name)
516 516 except Exception as e:
517 517 log.log('Skiping folder {}'.format(folder), self.name)
518 518 continue
519 519 return
520 520
521 521 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
522 522 expLabel='', last=False):
523 523
524 524 for path in folders:
525 525 files = glob.glob1(path, '*{}'.format(ext))
526 526 files.sort()
527 527 if last:
528 528 if files:
529 529 fo = files[-1]
530 530 try:
531 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
531 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
532 532 yield os.path.join(path, expLabel, fo)
533 533 except Exception as e:
534 534 pass
535 535 return
536 536 else:
537 537 return
538 538
539 539 for fo in files:
540 540 try:
541 541 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
542 542 if dt >= startDate and dt <= endDate:
543 543 yield os.path.join(path, expLabel, fo)
544 544 else:
545 545 log.log('Skiping file {}'.format(fo), self.name)
546 546 except Exception as e:
547 547 log.log('Skiping file {}'.format(fo), self.name)
548 548 continue
549 549
550 550 def searchFilesOffLine(self, path, startDate, endDate,
551 551 expLabel, ext, walk,
552 552 filefmt, folderfmt):
553 553 """Search files in offline mode for the given arguments
554 554
555 555 Return:
556 556 Generator of files
557 557 """
558 558
559 559 if walk:
560 560 folders = self.find_folders(
561 561 path, startDate, endDate, folderfmt)
562 562 else:
563 563 folders = path.split(',')
564 564
565 565 return self.find_files(
566 566 folders, ext, filefmt, startDate, endDate, expLabel)
567 567
568 568 def searchFilesOnLine(self, path, startDate, endDate,
569 569 expLabel, ext, walk,
570 570 filefmt, folderfmt):
571 571 """Search for the last file of the last folder
572 572
573 573 Arguments:
574 574 path : carpeta donde estan contenidos los files que contiene data
575 575 expLabel : Nombre del subexperimento (subfolder)
576 576 ext : extension de los files
577 577 walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath)
578 578
579 579 Return:
580 580 generator with the full path of last filename
581 581 """
582 582
583 583 if walk:
584 584 folders = self.find_folders(
585 585 path, startDate, endDate, folderfmt, last=True)
586 586 else:
587 587 folders = path.split(',')
588 588
589 589 return self.find_files(
590 590 folders, ext, filefmt, startDate, endDate, expLabel, last=True)
591 591
592 592 def setNextFile(self):
593 593 """Set the next file to be readed open it and parse de file header"""
594 594
595 595 while True:
596 596 if self.fp != None:
597 597 self.fp.close()
598 598
599 599 if self.online:
600 600 newFile = self.setNextFileOnline()
601 601 else:
602 602 newFile = self.setNextFileOffline()
603 603
604 604 if not(newFile):
605 605 if self.online:
606 606 raise schainpy.admin.SchainError('Time to wait for new files reach')
607 607 else:
608 608 if self.fileIndex == -1:
609 609 raise schainpy.admin.SchainWarning('No files found in the given path')
610 610 else:
611 611 raise schainpy.admin.SchainWarning('No more files to read')
612 612
613 613 if self.verifyFile(self.filename):
614 614 break
615 615
616 616 log.log('Opening file: %s' % self.filename, self.name)
617 617
618 618 self.readFirstHeader()
619 619 self.nReadBlocks = 0
620 620
621 621 def setNextFileOnline(self):
622 622 """Check for the next file to be readed in online mode.
623 623
624 624 Set:
625 625 self.filename
626 626 self.fp
627 627 self.filesize
628 628
629 629 Return:
630 630 boolean
631 631
632 632 """
633 633 nextFile = True
634 634 nextDay = False
635 635
636 636 for nFiles in range(self.nFiles+1):
637 637 for nTries in range(self.nTries):
638 638 fullfilename, filename = self.checkForRealPath(nextFile, nextDay)
639 639 if fullfilename is not None:
640 640 break
641 641 log.warning(
642 642 "Waiting %0.2f sec for the next file: \"%s\" , try %02d ..." % (self.delay, filename, nTries + 1),
643 643 self.name)
644 644 time.sleep(self.delay)
645 645 nextFile = False
646 646 continue
647 647
648 648 if fullfilename is not None:
649 649 break
650 650
651 651 self.nTries = 1
652 652 nextFile = True
653 653
654 654 if nFiles == (self.nFiles - 1):
655 655 log.log('Trying with next day...', self.name)
656 656 nextDay = True
657 657 self.nTries = 3
658 658
659 659 if fullfilename:
660 660 self.fileSize = os.path.getsize(fullfilename)
661 661 self.filename = fullfilename
662 662 self.flagIsNewFile = 1
663 663 if self.fp != None:
664 664 self.fp.close()
665 665 self.fp = self.open_file(fullfilename, self.open_mode)
666 666 self.flagNoMoreFiles = 0
667 667 self.fileIndex += 1
668 668 return 1
669 669 else:
670 670 return 0
671 671
672 672 def setNextFileOffline(self):
673 673 """Open the next file to be readed in offline mode"""
674 674
675 675 try:
676 676 filename = next(self.filenameList)
677 677 self.fileIndex +=1
678 678 except StopIteration:
679 679 self.flagNoMoreFiles = 1
680 680 return 0
681 681
682 682 self.filename = filename
683 683 self.fileSize = os.path.getsize(filename)
684 684 self.fp = self.open_file(filename, self.open_mode)
685 685 self.flagIsNewFile = 1
686 686
687 687 return 1
688 688
689 689 @staticmethod
690 690 def isDateTimeInRange(dt, startDate, endDate, startTime, endTime):
691 691 """Check if the given datetime is in range"""
692 692
693 693 if startDate <= dt.date() <= endDate:
694 694 if startTime <= dt.time() <= endTime:
695 695 return True
696 696 return False
697 697
698 698 def verifyFile(self, filename):
699 699 """Check for a valid file
700 700
701 701 Arguments:
702 702 filename -- full path filename
703 703
704 704 Return:
705 705 boolean
706 706 """
707 707
708 708 return True
709 709
710 710 def checkForRealPath(self, nextFile, nextDay):
711 711 """Check if the next file to be readed exists"""
712 712
713 713 raise NotImplementedError
714 714
715 715 def readFirstHeader(self):
716 716 """Parse the file header"""
717 717
718 718 pass
719 719
720 def waitDataBlock(self, pointer_location, blocksize=None):
721 """
722 """
723
724 currentPointer = pointer_location
725 if blocksize is None:
726 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
727 else:
728 neededSize = blocksize
729
730 for nTries in range(self.nTries):
731 self.fp.close()
732 self.fp = open(self.filename, 'rb')
733 self.fp.seek(currentPointer)
734
735 self.fileSize = os.path.getsize(self.filename)
736 currentSize = self.fileSize - currentPointer
737
738 if (currentSize >= neededSize):
739 return 1
740
741 log.warning(
742 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
743 self.name
744 )
745 time.sleep(self.delay)
746
747 return 0
748
720 749 class JRODataReader(Reader):
721 750
722 751 utc = 0
723 752 nReadBlocks = 0
724 753 foldercounter = 0
725 754 firstHeaderSize = 0
726 755 basicHeaderSize = 24
727 756 __isFirstTimeOnline = 1
728 757 filefmt = "*%Y%j***"
729 758 folderfmt = "*%Y%j"
730 759 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'online', 'delay', 'walk']
731 760
732 761 def getDtypeWidth(self):
733 762
734 763 dtype_index = get_dtype_index(self.dtype)
735 764 dtype_width = get_dtype_width(dtype_index)
736 765
737 766 return dtype_width
738 767
739 768 def checkForRealPath(self, nextFile, nextDay):
740 769 """Check if the next file to be readed exists.
741 770
742 771 Example :
743 772 nombre correcto del file es .../.../D2009307/P2009307367.ext
744 773
745 774 Entonces la funcion prueba con las siguientes combinaciones
746 775 .../.../y2009307367.ext
747 776 .../.../Y2009307367.ext
748 777 .../.../x2009307/y2009307367.ext
749 778 .../.../x2009307/Y2009307367.ext
750 779 .../.../X2009307/y2009307367.ext
751 780 .../.../X2009307/Y2009307367.ext
752 781 siendo para este caso, la ultima combinacion de letras, identica al file buscado
753 782
754 783 Return:
755 784 str -- fullpath of the file
756 785 """
757 786
758 787
759 788 if nextFile:
760 789 self.set += 1
761 790 if nextDay:
762 791 self.set = 0
763 792 self.doy += 1
764 793 foldercounter = 0
765 794 prefixDirList = [None, 'd', 'D']
766 795 if self.ext.lower() == ".r": # voltage
767 796 prefixFileList = ['d', 'D']
768 797 elif self.ext.lower() == ".pdata": # spectra
769 798 prefixFileList = ['p', 'P']
770 799
771 800 # barrido por las combinaciones posibles
772 801 for prefixDir in prefixDirList:
773 802 thispath = self.path
774 803 if prefixDir != None:
775 804 # formo el nombre del directorio xYYYYDDD (x=d o x=D)
776 805 if foldercounter == 0:
777 806 thispath = os.path.join(self.path, "%s%04d%03d" %
778 807 (prefixDir, self.year, self.doy))
779 808 else:
780 809 thispath = os.path.join(self.path, "%s%04d%03d_%02d" % (
781 810 prefixDir, self.year, self.doy, foldercounter))
782 811 for prefixFile in prefixFileList: # barrido por las dos combinaciones posibles de "D"
783 812 # formo el nombre del file xYYYYDDDSSS.ext
784 813 filename = "%s%04d%03d%03d%s" % (prefixFile, self.year, self.doy, self.set, self.ext)
785 814 fullfilename = os.path.join(
786 815 thispath, filename)
787 816
788 817 if os.path.exists(fullfilename):
789 818 return fullfilename, filename
790 819
791 820 return None, filename
792 821
793 822 def __waitNewBlock(self):
794 823 """
795 824 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
796 825
797 826 Si el modo de lectura es OffLine siempre retorn 0
798 827 """
799 828 if not self.online:
800 829 return 0
801 830
802 831 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
803 832 return 0
804 833
805 834 currentPointer = self.fp.tell()
806 835
807 836 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
808 837
809 838 for nTries in range(self.nTries):
810 839
811 840 self.fp.close()
812 841 self.fp = open(self.filename, 'rb')
813 842 self.fp.seek(currentPointer)
814 843
815 844 self.fileSize = os.path.getsize(self.filename)
816 845 currentSize = self.fileSize - currentPointer
817 846
818 847 if (currentSize >= neededSize):
819 848 self.basicHeaderObj.read(self.fp)
820 849 return 1
821 850
822 851 if self.fileSize == self.fileSizeByHeader:
823 852 # self.flagEoF = True
824 853 return 0
825 854
826 855 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
827 856 time.sleep(self.delay)
828 857
829 858 return 0
830 859
831 def waitDataBlock(self, pointer_location, blocksize=None):
832
833 currentPointer = pointer_location
834 if blocksize is None:
835 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
836 else:
837 neededSize = blocksize
838
839 for nTries in range(self.nTries):
840 self.fp.close()
841 self.fp = open(self.filename, 'rb')
842 self.fp.seek(currentPointer)
843
844 self.fileSize = os.path.getsize(self.filename)
845 currentSize = self.fileSize - currentPointer
846
847 if (currentSize >= neededSize):
848 return 1
849
850 log.warning(
851 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
852 self.name
853 )
854 time.sleep(self.delay)
855
856 return 0
857
858 860 def __setNewBlock(self):
859 861
860 862 if self.fp == None:
861 863 return 0
862 864
863 865 if self.flagIsNewFile:
864 866 self.lastUTTime = self.basicHeaderObj.utc
865 867 return 1
866 868
867 869 if self.realtime:
868 870 self.flagDiscontinuousBlock = 1
869 871 if not(self.setNextFile()):
870 872 return 0
871 873 else:
872 874 return 1
873 875
874 876 currentSize = self.fileSize - self.fp.tell()
875 877 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
876 878
877 879 if (currentSize >= neededSize):
878 880 self.basicHeaderObj.read(self.fp)
879 881 self.lastUTTime = self.basicHeaderObj.utc
880 882 return 1
881 883
882 884 if self.__waitNewBlock():
883 885 self.lastUTTime = self.basicHeaderObj.utc
884 886 return 1
885 887
886 888 if not(self.setNextFile()):
887 889 return 0
888 890
889 891 deltaTime = self.basicHeaderObj.utc - self.lastUTTime
890 892 self.lastUTTime = self.basicHeaderObj.utc
891 893
892 894 self.flagDiscontinuousBlock = 0
893 895
894 896 if deltaTime > self.maxTimeStep:
895 897 self.flagDiscontinuousBlock = 1
896 898
897 899 return 1
898 900
899 901 def readNextBlock(self):
900 902
901 903 while True:
902 904 if not(self.__setNewBlock()):
903 905 continue
904 906
905 907 if not(self.readBlock()):
906 908 return 0
907 909
908 910 self.getBasicHeader()
909 911
910 912 if not self.isDateTimeInRange(self.dataOut.datatime, self.startDate, self.endDate, self.startTime, self.endTime):
911 913 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.nReadBlocks,
912 914 self.processingHeaderObj.dataBlocksPerFile,
913 915 self.dataOut.datatime.ctime()))
914 916 continue
915 917
916 918 break
917 919
918 920 if self.verbose:
919 921 print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
920 922 self.processingHeaderObj.dataBlocksPerFile,
921 923 self.dataOut.datatime.ctime()))
922 924 return 1
923 925
924 926 def readFirstHeader(self):
925 927
926 928 self.basicHeaderObj.read(self.fp)
927 929 self.systemHeaderObj.read(self.fp)
928 930 self.radarControllerHeaderObj.read(self.fp)
929 931 self.processingHeaderObj.read(self.fp)
930 932 self.firstHeaderSize = self.basicHeaderObj.size
931 933
932 934 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
933 935 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
934 936 if datatype == 0:
935 937 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
936 938 elif datatype == 1:
937 939 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
938 940 elif datatype == 2:
939 941 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
940 942 elif datatype == 3:
941 943 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
942 944 elif datatype == 4:
943 945 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
944 946 elif datatype == 5:
945 947 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
946 948 else:
947 949 raise ValueError('Data type was not defined')
948 950
949 951 self.dtype = datatype_str
950 952 #self.ippSeconds = 2 * 1000 * self.radarControllerHeaderObj.ipp / self.c
951 953 self.fileSizeByHeader = self.processingHeaderObj.dataBlocksPerFile * self.processingHeaderObj.blockSize + \
952 954 self.firstHeaderSize + self.basicHeaderSize * \
953 955 (self.processingHeaderObj.dataBlocksPerFile - 1)
954 956 # self.dataOut.channelList = numpy.arange(self.systemHeaderObj.numChannels)
955 957 # self.dataOut.channelIndexList = numpy.arange(self.systemHeaderObj.numChannels)
956 958 self.getBlockDimension()
957 959
958 960 def verifyFile(self, filename):
959 961
960 962 flag = True
961 963
962 964 try:
963 965 fp = open(filename, 'rb')
964 966 except IOError:
965 967 log.error("File {} can't be opened".format(filename), self.name)
966 968 return False
967 969
968 970 if self.online and self.waitDataBlock(0):
969 971 pass
970 972
971 973 basicHeaderObj = BasicHeader(LOCALTIME)
972 974 systemHeaderObj = SystemHeader()
973 975 radarControllerHeaderObj = RadarControllerHeader()
974 976 processingHeaderObj = ProcessingHeader()
975 977
976 978 if not(basicHeaderObj.read(fp)):
977 979 flag = False
978 980 if not(systemHeaderObj.read(fp)):
979 981 flag = False
980 982 if not(radarControllerHeaderObj.read(fp)):
981 983 flag = False
982 984 if not(processingHeaderObj.read(fp)):
983 985 flag = False
984 986 if not self.online:
985 987 dt1 = basicHeaderObj.datatime
986 988 pos = self.fileSize-processingHeaderObj.blockSize-24
987 989 if pos<0:
988 990 flag = False
989 991 log.error('Invalid size for file: {}'.format(self.filename), self.name)
990 992 else:
991 993 fp.seek(pos)
992 994 if not(basicHeaderObj.read(fp)):
993 995 flag = False
994 996 dt2 = basicHeaderObj.datatime
995 997 if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
996 998 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
997 999 flag = False
998 1000
999 1001 fp.close()
1000 1002 return flag
1001 1003
1002 1004 def findDatafiles(self, path, startDate=None, endDate=None, expLabel='', ext='.r', walk=True, include_path=False):
1003 1005
1004 1006 path_empty = True
1005 1007
1006 1008 dateList = []
1007 1009 pathList = []
1008 1010
1009 1011 multi_path = path.split(',')
1010 1012
1011 1013 if not walk:
1012 1014
1013 1015 for single_path in multi_path:
1014 1016
1015 1017 if not os.path.isdir(single_path):
1016 1018 continue
1017 1019
1018 1020 fileList = glob.glob1(single_path, "*" + ext)
1019 1021
1020 1022 if not fileList:
1021 1023 continue
1022 1024
1023 1025 path_empty = False
1024 1026
1025 1027 fileList.sort()
1026 1028
1027 1029 for thisFile in fileList:
1028 1030
1029 1031 if not os.path.isfile(os.path.join(single_path, thisFile)):
1030 1032 continue
1031 1033
1032 1034 if not isRadarFile(thisFile):
1033 1035 continue
1034 1036
1035 1037 if not isFileInDateRange(thisFile, startDate, endDate):
1036 1038 continue
1037 1039
1038 1040 thisDate = getDateFromRadarFile(thisFile)
1039 1041
1040 1042 if thisDate in dateList or single_path in pathList:
1041 1043 continue
1042 1044
1043 1045 dateList.append(thisDate)
1044 1046 pathList.append(single_path)
1045 1047
1046 1048 else:
1047 1049 for single_path in multi_path:
1048 1050
1049 1051 if not os.path.isdir(single_path):
1050 1052 continue
1051 1053
1052 1054 dirList = []
1053 1055
1054 1056 for thisPath in os.listdir(single_path):
1055 1057
1056 1058 if not os.path.isdir(os.path.join(single_path, thisPath)):
1057 1059 continue
1058 1060
1059 1061 if not isRadarFolder(thisPath):
1060 1062 continue
1061 1063
1062 1064 if not isFolderInDateRange(thisPath, startDate, endDate):
1063 1065 continue
1064 1066
1065 1067 dirList.append(thisPath)
1066 1068
1067 1069 if not dirList:
1068 1070 continue
1069 1071
1070 1072 dirList.sort()
1071 1073
1072 1074 for thisDir in dirList:
1073 1075
1074 1076 datapath = os.path.join(single_path, thisDir, expLabel)
1075 1077 fileList = glob.glob1(datapath, "*" + ext)
1076 1078
1077 1079 if not fileList:
1078 1080 continue
1079 1081
1080 1082 path_empty = False
1081 1083
1082 1084 thisDate = getDateFromRadarFolder(thisDir)
1083 1085
1084 1086 pathList.append(datapath)
1085 1087 dateList.append(thisDate)
1086 1088
1087 1089 dateList.sort()
1088 1090
1089 1091 if walk:
1090 1092 pattern_path = os.path.join(multi_path[0], "[dYYYYDDD]", expLabel)
1091 1093 else:
1092 1094 pattern_path = multi_path[0]
1093 1095
1094 1096 if path_empty:
1095 1097 raise schainpy.admin.SchainError("[Reading] No *%s files in %s for %s to %s" % (ext, pattern_path, startDate, endDate))
1096 1098 else:
1097 1099 if not dateList:
1098 1100 raise schainpy.admin.SchainError("[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" % (startDate, endDate, ext, path))
1099 1101
1100 1102 if include_path:
1101 1103 return dateList, pathList
1102 1104
1103 1105 return dateList
1104 1106
1105 1107 def setup(self, **kwargs):
1106 1108
1107 1109 self.set_kwargs(**kwargs)
1108 1110 if not self.ext.startswith('.'):
1109 1111 self.ext = '.{}'.format(self.ext)
1110 1112
1111 1113 if self.server is not None:
1112 1114 if 'tcp://' in self.server:
1113 1115 address = server
1114 1116 else:
1115 1117 address = 'ipc:///tmp/%s' % self.server
1116 1118 self.server = address
1117 1119 self.context = zmq.Context()
1118 1120 self.receiver = self.context.socket(zmq.PULL)
1119 1121 self.receiver.connect(self.server)
1120 1122 time.sleep(0.5)
1121 1123 print('[Starting] ReceiverData from {}'.format(self.server))
1122 1124 else:
1123 1125 self.server = None
1124 1126 if self.path == None:
1125 1127 raise ValueError("[Reading] The path is not valid")
1126 1128
1127 1129 if self.online:
1128 1130 log.log("[Reading] Searching files in online mode...", self.name)
1129 1131
1130 1132 for nTries in range(self.nTries):
1131 1133 fullpath = self.searchFilesOnLine(self.path, self.startDate,
1132 1134 self.endDate, self.expLabel, self.ext, self.walk,
1133 1135 self.filefmt, self.folderfmt)
1134 1136
1135 1137 try:
1136 1138 fullpath = next(fullpath)
1137 1139 except:
1138 1140 fullpath = None
1139 1141
1140 1142 if fullpath:
1141 1143 break
1142 1144
1143 1145 log.warning(
1144 1146 'Waiting {} sec for a valid file in {}: try {} ...'.format(
1145 1147 self.delay, self.path, nTries + 1),
1146 1148 self.name)
1147 1149 time.sleep(self.delay)
1148 1150
1149 1151 if not(fullpath):
1150 1152 raise schainpy.admin.SchainError(
1151 1153 'There isn\'t any valid file in {}'.format(self.path))
1152 1154
1153 1155 pathname, filename = os.path.split(fullpath)
1154 1156 self.year = int(filename[1:5])
1155 1157 self.doy = int(filename[5:8])
1156 1158 self.set = int(filename[8:11]) - 1
1157 1159 else:
1158 1160 log.log("Searching files in {}".format(self.path), self.name)
1159 1161 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1160 1162 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
1161 1163
1162 1164 self.setNextFile()
1163 1165
1164 1166 return
1165 1167
1166 1168 def getBasicHeader(self):
1167 1169
1168 1170 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
1169 1171 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
1170 1172
1171 1173 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
1172 1174
1173 1175 self.dataOut.timeZone = self.basicHeaderObj.timeZone
1174 1176
1175 1177 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
1176 1178
1177 1179 self.dataOut.errorCount = self.basicHeaderObj.errorCount
1178 1180
1179 1181 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1180 1182
1181 1183 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1182 1184
1183 1185 # self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock*self.nTxs
1184 1186
1185 1187 def getFirstHeader(self):
1186 1188
1187 1189 raise NotImplementedError
1188 1190
1189 1191 def getData(self):
1190 1192
1191 1193 raise NotImplementedError
1192 1194
1193 1195 def hasNotDataInBuffer(self):
1194 1196
1195 1197 raise NotImplementedError
1196 1198
1197 1199 def readBlock(self):
1198 1200
1199 1201 raise NotImplementedError
1200 1202
1201 1203 def isEndProcess(self):
1202 1204
1203 1205 return self.flagNoMoreFiles
1204 1206
1205 1207 def printReadBlocks(self):
1206 1208
1207 1209 print("[Reading] Number of read blocks per file %04d" % self.nReadBlocks)
1208 1210
1209 1211 def printTotalBlocks(self):
1210 1212
1211 1213 print("[Reading] Number of read blocks %04d" % self.nTotalBlocks)
1212 1214
1213 1215 def run(self, **kwargs):
1214 1216 """
1215 1217
1216 1218 Arguments:
1217 1219 path :
1218 1220 startDate :
1219 1221 endDate :
1220 1222 startTime :
1221 1223 endTime :
1222 1224 set :
1223 1225 expLabel :
1224 1226 ext :
1225 1227 online :
1226 1228 delay :
1227 1229 walk :
1228 1230 getblock :
1229 1231 nTxs :
1230 1232 realtime :
1231 1233 blocksize :
1232 1234 blocktime :
1233 1235 skip :
1234 1236 cursor :
1235 1237 warnings :
1236 1238 server :
1237 1239 verbose :
1238 1240 format :
1239 1241 oneDDict :
1240 1242 twoDDict :
1241 1243 independentParam :
1242 1244 """
1243 1245
1244 1246 if not(self.isConfig):
1245 1247 self.setup(**kwargs)
1246 1248 self.isConfig = True
1247 1249 if self.server is None:
1248 1250 self.getData()
1249 1251 else:
1250 1252 self.getFromServer()
1251 1253
1252 1254
1253 1255 class JRODataWriter(Reader):
1254 1256
1255 1257 """
1256 1258 Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura
1257 1259 de los datos siempre se realiza por bloques.
1258 1260 """
1259 1261
1260 1262 setFile = None
1261 1263 profilesPerBlock = None
1262 1264 blocksPerFile = None
1263 1265 nWriteBlocks = 0
1264 1266 fileDate = None
1265 1267
1266 1268 def __init__(self, dataOut=None):
1267 1269 raise NotImplementedError
1268 1270
1269 1271 def hasAllDataInBuffer(self):
1270 1272 raise NotImplementedError
1271 1273
1272 1274 def setBlockDimension(self):
1273 1275 raise NotImplementedError
1274 1276
1275 1277 def writeBlock(self):
1276 1278 raise NotImplementedError
1277 1279
1278 1280 def putData(self):
1279 1281 raise NotImplementedError
1280 1282
1281 1283 def getDtypeWidth(self):
1282 1284
1283 1285 dtype_index = get_dtype_index(self.dtype)
1284 1286 dtype_width = get_dtype_width(dtype_index)
1285 1287
1286 1288 return dtype_width
1287 1289
1288 1290 def getProcessFlags(self):
1289 1291
1290 1292 processFlags = 0
1291 1293
1292 1294 dtype_index = get_dtype_index(self.dtype)
1293 1295 procflag_dtype = get_procflag_dtype(dtype_index)
1294 1296
1295 1297 processFlags += procflag_dtype
1296 1298
1297 1299 if self.dataOut.flagDecodeData:
1298 1300 processFlags += PROCFLAG.DECODE_DATA
1299 1301
1300 1302 if self.dataOut.flagDeflipData:
1301 1303 processFlags += PROCFLAG.DEFLIP_DATA
1302 1304
1303 1305 if self.dataOut.code is not None:
1304 1306 processFlags += PROCFLAG.DEFINE_PROCESS_CODE
1305 1307
1306 1308 if self.dataOut.nCohInt > 1:
1307 1309 processFlags += PROCFLAG.COHERENT_INTEGRATION
1308 1310
1309 1311 if self.dataOut.type == "Spectra":
1310 1312 if self.dataOut.nIncohInt > 1:
1311 1313 processFlags += PROCFLAG.INCOHERENT_INTEGRATION
1312 1314
1313 1315 if self.dataOut.data_dc is not None:
1314 1316 processFlags += PROCFLAG.SAVE_CHANNELS_DC
1315 1317
1316 1318 if self.dataOut.flagShiftFFT:
1317 1319 processFlags += PROCFLAG.SHIFT_FFT_DATA
1318 1320
1319 1321 return processFlags
1320 1322
1321 1323 def setBasicHeader(self):
1322 1324
1323 1325 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1324 1326 self.basicHeaderObj.version = self.versionFile
1325 1327 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1326 1328 utc = numpy.floor(self.dataOut.utctime)
1327 1329 milisecond = (self.dataOut.utctime - utc) * 1000.0
1328 1330 self.basicHeaderObj.utc = utc
1329 1331 self.basicHeaderObj.miliSecond = milisecond
1330 1332 self.basicHeaderObj.timeZone = self.dataOut.timeZone
1331 1333 self.basicHeaderObj.dstFlag = self.dataOut.dstFlag
1332 1334 self.basicHeaderObj.errorCount = self.dataOut.errorCount
1333 1335
1334 1336 def setFirstHeader(self):
1335 1337 """
1336 1338 Obtiene una copia del First Header
1337 1339
1338 1340 Affected:
1339 1341
1340 1342 self.basicHeaderObj
1341 1343 self.systemHeaderObj
1342 1344 self.radarControllerHeaderObj
1343 1345 self.processingHeaderObj self.
1344 1346
1345 1347 Return:
1346 1348 None
1347 1349 """
1348 1350
1349 1351 raise NotImplementedError
1350 1352
1351 1353 def __writeFirstHeader(self):
1352 1354 """
1353 1355 Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader)
1354 1356
1355 1357 Affected:
1356 1358 __dataType
1357 1359
1358 1360 Return:
1359 1361 None
1360 1362 """
1361 1363
1362 1364 # CALCULAR PARAMETROS
1363 1365
1364 1366 sizeLongHeader = self.systemHeaderObj.size + \
1365 1367 self.radarControllerHeaderObj.size + self.processingHeaderObj.size
1366 1368 self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader
1367 1369
1368 1370 self.basicHeaderObj.write(self.fp)
1369 1371 self.systemHeaderObj.write(self.fp)
1370 1372 self.radarControllerHeaderObj.write(self.fp)
1371 1373 self.processingHeaderObj.write(self.fp)
1372 1374
1373 1375 def __setNewBlock(self):
1374 1376 """
1375 1377 Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header
1376 1378
1377 1379 Return:
1378 1380 0 : si no pudo escribir nada
1379 1381 1 : Si escribio el Basic el First Header
1380 1382 """
1381 1383 if self.fp == None:
1382 1384 self.setNextFile()
1383 1385
1384 1386 if self.flagIsNewFile:
1385 1387 return 1
1386 1388
1387 1389 if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile:
1388 1390 self.basicHeaderObj.write(self.fp)
1389 1391 return 1
1390 1392
1391 1393 if not(self.setNextFile()):
1392 1394 return 0
1393 1395
1394 1396 return 1
1395 1397
1396 1398 def writeNextBlock(self):
1397 1399 """
1398 1400 Selecciona el bloque siguiente de datos y los escribe en un file
1399 1401
1400 1402 Return:
1401 1403 0 : Si no hizo pudo escribir el bloque de datos
1402 1404 1 : Si no pudo escribir el bloque de datos
1403 1405 """
1404 1406 if not(self.__setNewBlock()):
1405 1407 return 0
1406 1408
1407 1409 self.writeBlock()
1408 1410
1409 1411 print("[Writing] Block No. %d/%d" % (self.blockIndex,
1410 1412 self.processingHeaderObj.dataBlocksPerFile))
1411 1413
1412 1414 return 1
1413 1415
1414 1416 def setNextFile(self):
1415 1417 """Determina el siguiente file que sera escrito
1416 1418
1417 1419 Affected:
1418 1420 self.filename
1419 1421 self.subfolder
1420 1422 self.fp
1421 1423 self.setFile
1422 1424 self.flagIsNewFile
1423 1425
1424 1426 Return:
1425 1427 0 : Si el archivo no puede ser escrito
1426 1428 1 : Si el archivo esta listo para ser escrito
1427 1429 """
1428 1430 ext = self.ext
1429 1431 path = self.path
1430 1432
1431 1433 if self.fp != None:
1432 1434 self.fp.close()
1433 1435
1434 1436 if not os.path.exists(path):
1435 1437 os.mkdir(path)
1436 1438
1437 1439 timeTuple = time.localtime(self.dataOut.utctime)
1438 1440 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
1439 1441
1440 1442 fullpath = os.path.join(path, subfolder)
1441 1443 setFile = self.setFile
1442 1444
1443 1445 if not(os.path.exists(fullpath)):
1444 1446 os.mkdir(fullpath)
1445 1447 setFile = -1 # inicializo mi contador de seteo
1446 1448 else:
1447 1449 filesList = os.listdir(fullpath)
1448 1450 if len(filesList) > 0:
1449 1451 filesList = sorted(filesList, key=str.lower)
1450 1452 filen = filesList[-1]
1451 1453 # el filename debera tener el siguiente formato
1452 1454 # 0 1234 567 89A BCDE (hex)
1453 1455 # x YYYY DDD SSS .ext
1454 1456 if isNumber(filen[8:11]):
1455 1457 # inicializo mi contador de seteo al seteo del ultimo file
1456 1458 setFile = int(filen[8:11])
1457 1459 else:
1458 1460 setFile = -1
1459 1461 else:
1460 1462 setFile = -1 # inicializo mi contador de seteo
1461 1463
1462 1464 setFile += 1
1463 1465
1464 1466 # If this is a new day it resets some values
1465 1467 if self.dataOut.datatime.date() > self.fileDate:
1466 1468 setFile = 0
1467 1469 self.nTotalBlocks = 0
1468 1470
1469 1471 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1470 1472 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1471 1473
1472 1474 filename = os.path.join(path, subfolder, filen)
1473 1475
1474 1476 fp = open(filename, 'wb')
1475 1477
1476 1478 self.blockIndex = 0
1477 1479 self.filename = filename
1478 1480 self.subfolder = subfolder
1479 1481 self.fp = fp
1480 1482 self.setFile = setFile
1481 1483 self.flagIsNewFile = 1
1482 1484 self.fileDate = self.dataOut.datatime.date()
1483 1485 self.setFirstHeader()
1484 1486
1485 1487 print('[Writing] Opening file: %s' % self.filename)
1486 1488
1487 1489 self.__writeFirstHeader()
1488 1490
1489 1491 return 1
1490 1492
1491 1493 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1492 1494 """
1493 1495 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1494 1496
1495 1497 Inputs:
1496 1498 path : directory where data will be saved
1497 1499 profilesPerBlock : number of profiles per block
1498 1500 set : initial file set
1499 1501 datatype : An integer number that defines data type:
1500 1502 0 : int8 (1 byte)
1501 1503 1 : int16 (2 bytes)
1502 1504 2 : int32 (4 bytes)
1503 1505 3 : int64 (8 bytes)
1504 1506 4 : float32 (4 bytes)
1505 1507 5 : double64 (8 bytes)
1506 1508
1507 1509 Return:
1508 1510 0 : Si no realizo un buen seteo
1509 1511 1 : Si realizo un buen seteo
1510 1512 """
1511 1513
1512 1514 if ext == None:
1513 1515 ext = self.ext
1514 1516
1515 1517 self.ext = ext.lower()
1516 1518
1517 1519 self.path = path
1518 1520
1519 1521 if set is None:
1520 1522 self.setFile = -1
1521 1523 else:
1522 1524 self.setFile = set - 1
1523 1525
1524 1526 self.blocksPerFile = blocksPerFile
1525 1527 self.profilesPerBlock = profilesPerBlock
1526 1528 self.dataOut = dataOut
1527 1529 self.fileDate = self.dataOut.datatime.date()
1528 1530 self.dtype = self.dataOut.dtype
1529 1531
1530 1532 if datatype is not None:
1531 1533 self.dtype = get_numpy_dtype(datatype)
1532 1534
1533 1535 if not(self.setNextFile()):
1534 1536 print("[Writing] There isn't a next file")
1535 1537 return 0
1536 1538
1537 1539 self.setBlockDimension()
1538 1540
1539 1541 return 1
1540 1542
1541 1543 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1542 1544
1543 1545 if not(self.isConfig):
1544 1546
1545 1547 self.setup(dataOut, path, blocksPerFile, profilesPerBlock=profilesPerBlock,
1546 1548 set=set, ext=ext, datatype=datatype, **kwargs)
1547 1549 self.isConfig = True
1548 1550
1549 1551 self.dataOut = dataOut
1550 1552 self.putData()
1551 1553 return self.dataOut
1552 1554
1553 1555 @MPDecorator
1554 1556 class printInfo(Operation):
1555 1557
1556 1558 def __init__(self):
1557 1559
1558 1560 Operation.__init__(self)
1559 1561 self.__printInfo = True
1560 1562
1561 1563 def run(self, dataOut, headers = ['systemHeaderObj', 'radarControllerHeaderObj', 'processingHeaderObj']):
1562 1564 if self.__printInfo == False:
1563 1565 return
1564 1566
1565 1567 for header in headers:
1566 1568 if hasattr(dataOut, header):
1567 1569 obj = getattr(dataOut, header)
1568 1570 if hasattr(obj, 'printInfo'):
1569 1571 obj.printInfo()
1570 1572 else:
1571 1573 print(obj)
1572 1574 else:
1573 1575 log.warning('Header {} Not found in object'.format(header))
1574 1576
1575 1577 self.__printInfo = False
General Comments 0
You need to be logged in to leave comments. Login now