##// END OF EJS Templates
schain
RSS
Add georef for realtime
Juan C. Espinoza -
r1498:675eff1be3dc
parent child
Show More
Side by Side Unified
Context file:
r1498:675eff1be3dc Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1069 +1,1074
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Definition of diferent Data objects for different types of data
6 6
7 7 Here you will find the diferent data objects for the different types
8 8 of data, this data objects must be used as dataIn or dataOut objects in
9 9 processing units and operations. Currently the supported data objects are:
10 10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 11 """
12 12
13 13 import copy
14 14 import numpy
15 15 import datetime
16 16 import json
17 17
18 18 import schainpy.admin
19 19 from schainpy.utils import log
20 20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 21 from schainpy.model.data import _noise
22 22
23 23
24 24 def getNumpyDtype(dataTypeCode):
25 25
26 26 if dataTypeCode == 0:
27 27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 28 elif dataTypeCode == 1:
29 29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 30 elif dataTypeCode == 2:
31 31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 32 elif dataTypeCode == 3:
33 33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 34 elif dataTypeCode == 4:
35 35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 36 elif dataTypeCode == 5:
37 37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 38 else:
39 39 raise ValueError('dataTypeCode was not defined')
40 40
41 41 return numpyDtype
42 42
43 43
44 44 def getDataTypeCode(numpyDtype):
45 45
46 46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 47 datatype = 0
48 48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 49 datatype = 1
50 50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 51 datatype = 2
52 52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 53 datatype = 3
54 54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 55 datatype = 4
56 56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 57 datatype = 5
58 58 else:
59 59 datatype = None
60 60
61 61 return datatype
62 62
63 63
64 64 def hildebrand_sekhon(data, navg):
65 65 """
66 66 This method is for the objective determination of the noise level in Doppler spectra. This
67 67 implementation technique is based on the fact that the standard deviation of the spectral
68 68 densities is equal to the mean spectral density for white Gaussian noise
69 69
70 70 Inputs:
71 71 Data : heights
72 72 navg : numbers of averages
73 73
74 74 Return:
75 75 mean : noise's level
76 76 """
77 77
78 78 sortdata = numpy.sort(data, axis=None)
79 79 '''
80 80 lenOfData = len(sortdata)
81 81 nums_min = lenOfData*0.2
82 82
83 83 if nums_min <= 5:
84 84
85 85 nums_min = 5
86 86
87 87 sump = 0.
88 88 sumq = 0.
89 89
90 90 j = 0
91 91 cont = 1
92 92
93 93 while((cont == 1)and(j < lenOfData)):
94 94
95 95 sump += sortdata[j]
96 96 sumq += sortdata[j]**2
97 97
98 98 if j > nums_min:
99 99 rtest = float(j)/(j-1) + 1.0/navg
100 100 if ((sumq*j) > (rtest*sump**2)):
101 101 j = j - 1
102 102 sump = sump - sortdata[j]
103 103 sumq = sumq - sortdata[j]**2
104 104 cont = 0
105 105
106 106 j += 1
107 107
108 108 lnoise = sump / j
109 109 '''
110 110 return _noise.hildebrand_sekhon(sortdata, navg)
111 111
112 112
113 113 class Beam:
114 114
115 115 def __init__(self):
116 116 self.codeList = []
117 117 self.azimuthList = []
118 118 self.zenithList = []
119 119
120 120
121 121 class GenericData(object):
122 122
123 123 flagNoData = True
124 124
125 125 def copy(self, inputObj=None):
126 126
127 127 if inputObj == None:
128 128 return copy.deepcopy(self)
129 129
130 130 for key in list(inputObj.__dict__.keys()):
131 131
132 132 attribute = inputObj.__dict__[key]
133 133
134 134 # If this attribute is a tuple or list
135 135 if type(inputObj.__dict__[key]) in (tuple, list):
136 136 self.__dict__[key] = attribute[:]
137 137 continue
138 138
139 139 # If this attribute is another object or instance
140 140 if hasattr(attribute, '__dict__'):
141 141 self.__dict__[key] = attribute.copy()
142 142 continue
143 143
144 144 self.__dict__[key] = inputObj.__dict__[key]
145 145
146 146 def deepcopy(self):
147 147
148 148 return copy.deepcopy(self)
149 149
150 150 def isEmpty(self):
151 151
152 152 return self.flagNoData
153 153
154 154 def isReady(self):
155 155
156 156 return not self.flagNoData
157 157
158 158
159 159 class JROData(GenericData):
160 160
161 161 systemHeaderObj = SystemHeader()
162 162 radarControllerHeaderObj = RadarControllerHeader()
163 163 type = None
164 164 datatype = None # dtype but in string
165 165 nProfiles = None
166 166 heightList = None
167 167 channelList = None
168 168 flagDiscontinuousBlock = False
169 169 useLocalTime = False
170 170 utctime = None
171 171 timeZone = None
172 172 dstFlag = None
173 173 errorCount = None
174 174 blocksize = None
175 175 flagDecodeData = False # asumo q la data no esta decodificada
176 176 flagDeflipData = False # asumo q la data no esta sin flip
177 177 flagShiftFFT = False
178 178 nCohInt = None
179 179 windowOfFilter = 1
180 180 C = 3e8
181 181 frequency = 49.92e6
182 182 realtime = False
183 183 beacon_heiIndexList = None
184 184 last_block = None
185 185 blocknow = None
186 186 azimuth = None
187 187 zenith = None
188 188 beam = Beam()
189 189 profileIndex = None
190 190 error = None
191 191 data = None
192 192 nmodes = None
193 193 h0 = 0
194 194 metadata_list = ['heightList', 'timeZone', 'type']
195 195
196 196 def __str__(self):
197 197
198 198 return '{} - {}'.format(self.type, self.datatime())
199 199
200 200 def getNoise(self):
201 201
202 202 raise NotImplementedError
203 203
204 204 @property
205 205 def nChannels(self):
206 206
207 207 return len(self.channelList)
208 208
209 209 @property
210 210 def channelIndexList(self):
211 211
212 212 return list(range(self.nChannels))
213 213
214 214 @property
215 215 def nHeights(self):
216 216
217 217 return len(self.heightList)
218 218
219 219 def getDeltaH(self):
220 220
221 221 return self.heightList[1] - self.heightList[0]
222 222
223 223 @property
224 224 def ltctime(self):
225 225
226 226 if self.useLocalTime:
227 227 return self.utctime - self.timeZone * 60
228 228
229 229 return self.utctime
230 230
231 231 @property
232 232 def datatime(self):
233 233
234 234 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
235 235 return datatimeValue
236 236
237 237 def getTimeRange(self):
238 238
239 239 datatime = []
240 240
241 241 datatime.append(self.ltctime)
242 242 datatime.append(self.ltctime + self.timeInterval + 1)
243 243
244 244 datatime = numpy.array(datatime)
245 245
246 246 return datatime
247 247
248 248 def getFmaxTimeResponse(self):
249 249
250 250 period = (10**-6) * self.getDeltaH() / (0.15)
251 251
252 252 PRF = 1. / (period * self.nCohInt)
253 253
254 254 fmax = PRF
255 255
256 256 return fmax
257 257
258 258 def getFmax(self):
259 259 PRF = 1. / (self.ippSeconds * self.nCohInt)
260 260
261 261 fmax = PRF
262 262 return fmax
263 263
264 264 def getVmax(self):
265 265
266 266 _lambda = self.C / self.frequency
267 267
268 268 vmax = self.getFmax() * _lambda / 2
269 269
270 270 return vmax
271 271
272 272 @property
273 273 def ippSeconds(self):
274 274 '''
275 275 '''
276 276 return self.radarControllerHeaderObj.ippSeconds
277 277
278 278 @ippSeconds.setter
279 279 def ippSeconds(self, ippSeconds):
280 280 '''
281 281 '''
282 282 self.radarControllerHeaderObj.ippSeconds = ippSeconds
283 283
284 284 @property
285 285 def code(self):
286 286 '''
287 287 '''
288 288 return self.radarControllerHeaderObj.code
289 289
290 290 @code.setter
291 291 def code(self, code):
292 292 '''
293 293 '''
294 294 self.radarControllerHeaderObj.code = code
295 295
296 296 @property
297 297 def nCode(self):
298 298 '''
299 299 '''
300 300 return self.radarControllerHeaderObj.nCode
301 301
302 302 @nCode.setter
303 303 def nCode(self, ncode):
304 304 '''
305 305 '''
306 306 self.radarControllerHeaderObj.nCode = ncode
307 307
308 308 @property
309 309 def nBaud(self):
310 310 '''
311 311 '''
312 312 return self.radarControllerHeaderObj.nBaud
313 313
314 314 @nBaud.setter
315 315 def nBaud(self, nbaud):
316 316 '''
317 317 '''
318 318 self.radarControllerHeaderObj.nBaud = nbaud
319 319
320 320 @property
321 321 def ipp(self):
322 322 '''
323 323 '''
324 324 return self.radarControllerHeaderObj.ipp
325 325
326 326 @ipp.setter
327 327 def ipp(self, ipp):
328 328 '''
329 329 '''
330 330 self.radarControllerHeaderObj.ipp = ipp
331 331
332 332 @property
333 333 def metadata(self):
334 334 '''
335 335 '''
336 336
337 337 return {attr: getattr(self, attr) for attr in self.metadata_list}
338 338
339 339
340 340 class Voltage(JROData):
341 341
342 342 dataPP_POW = None
343 343 dataPP_DOP = None
344 344 dataPP_WIDTH = None
345 345 dataPP_SNR = None
346 346
347 347 def __init__(self):
348 348 '''
349 349 Constructor
350 350 '''
351 351
352 352 self.useLocalTime = True
353 353 self.radarControllerHeaderObj = RadarControllerHeader()
354 354 self.systemHeaderObj = SystemHeader()
355 355 self.type = "Voltage"
356 356 self.data = None
357 357 self.nProfiles = None
358 358 self.heightList = None
359 359 self.channelList = None
360 360 self.flagNoData = True
361 361 self.flagDiscontinuousBlock = False
362 362 self.utctime = None
363 363 self.timeZone = 0
364 364 self.dstFlag = None
365 365 self.errorCount = None
366 366 self.nCohInt = None
367 367 self.blocksize = None
368 368 self.flagCohInt = False
369 369 self.flagDecodeData = False # asumo q la data no esta decodificada
370 370 self.flagDeflipData = False # asumo q la data no esta sin flip
371 371 self.flagShiftFFT = False
372 372 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
373 373 self.profileIndex = 0
374 374 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
375 375 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
376 376
377 377 def getNoisebyHildebrand(self, channel=None):
378 378 """
379 379 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
380 380
381 381 Return:
382 382 noiselevel
383 383 """
384 384
385 385 if channel != None:
386 386 data = self.data[channel]
387 387 nChannels = 1
388 388 else:
389 389 data = self.data
390 390 nChannels = self.nChannels
391 391
392 392 noise = numpy.zeros(nChannels)
393 393 power = data * numpy.conjugate(data)
394 394
395 395 for thisChannel in range(nChannels):
396 396 if nChannels == 1:
397 397 daux = power[:].real
398 398 else:
399 399 daux = power[thisChannel, :].real
400 400 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
401 401
402 402 return noise
403 403
404 404 def getNoise(self, type=1, channel=None):
405 405
406 406 if type == 1:
407 407 noise = self.getNoisebyHildebrand(channel)
408 408
409 409 return noise
410 410
411 411 def getPower(self, channel=None):
412 412
413 413 if channel != None:
414 414 data = self.data[channel]
415 415 else:
416 416 data = self.data
417 417
418 418 power = data * numpy.conjugate(data)
419 419 powerdB = 10 * numpy.log10(power.real)
420 420 powerdB = numpy.squeeze(powerdB)
421 421
422 422 return powerdB
423 423
424 424 @property
425 425 def timeInterval(self):
426 426
427 427 return self.ippSeconds * self.nCohInt
428 428
429 429 noise = property(getNoise, "I'm the 'nHeights' property.")
430 430
431 431
432 432 class Spectra(JROData):
433 433
434 434 def __init__(self):
435 435 '''
436 436 Constructor
437 437 '''
438 438
439 439 self.data_dc = None
440 440 self.data_spc = None
441 441 self.data_cspc = None
442 442 self.useLocalTime = True
443 443 self.radarControllerHeaderObj = RadarControllerHeader()
444 444 self.systemHeaderObj = SystemHeader()
445 445 self.type = "Spectra"
446 446 self.timeZone = 0
447 447 self.nProfiles = None
448 448 self.heightList = None
449 449 self.channelList = None
450 450 self.pairsList = None
451 451 self.flagNoData = True
452 452 self.flagDiscontinuousBlock = False
453 453 self.utctime = None
454 454 self.nCohInt = None
455 455 self.nIncohInt = None
456 456 self.blocksize = None
457 457 self.nFFTPoints = None
458 458 self.wavelength = None
459 459 self.flagDecodeData = False # asumo q la data no esta decodificada
460 460 self.flagDeflipData = False # asumo q la data no esta sin flip
461 461 self.flagShiftFFT = False
462 462 self.ippFactor = 1
463 463 self.beacon_heiIndexList = []
464 464 self.noise_estimation = None
465 465 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
466 466 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
467 467
468 468 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
469 469 """
470 470 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
471 471
472 472 Return:
473 473 noiselevel
474 474 """
475 475
476 476 noise = numpy.zeros(self.nChannels)
477 477
478 478 for channel in range(self.nChannels):
479 479 daux = self.data_spc[channel,
480 480 xmin_index:xmax_index, ymin_index:ymax_index]
481 481 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
482 482
483 483 return noise
484 484
485 485 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
486 486
487 487 if self.noise_estimation is not None:
488 488 # this was estimated by getNoise Operation defined in jroproc_spectra.py
489 489 return self.noise_estimation
490 490 else:
491 491 noise = self.getNoisebyHildebrand(
492 492 xmin_index, xmax_index, ymin_index, ymax_index)
493 493 return noise
494 494
495 495 def getFreqRangeTimeResponse(self, extrapoints=0):
496 496
497 497 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
498 498 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
499 499
500 500 return freqrange
501 501
502 502 def getAcfRange(self, extrapoints=0):
503 503
504 504 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
505 505 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
506 506
507 507 return freqrange
508 508
509 509 def getFreqRange(self, extrapoints=0):
510 510
511 511 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
512 512 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
513 513
514 514 return freqrange
515 515
516 516 def getVelRange(self, extrapoints=0):
517 517
518 518 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
519 519 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
520 520
521 521 if self.nmodes:
522 522 return velrange/self.nmodes
523 523 else:
524 524 return velrange
525 525
526 526 @property
527 527 def nPairs(self):
528 528
529 529 return len(self.pairsList)
530 530
531 531 @property
532 532 def pairsIndexList(self):
533 533
534 534 return list(range(self.nPairs))
535 535
536 536 @property
537 537 def normFactor(self):
538 538
539 539 pwcode = 1
540 540
541 541 if self.flagDecodeData:
542 542 pwcode = numpy.sum(self.code[0]**2)
543 543 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
544 544 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
545 545
546 546 return normFactor
547 547
548 548 @property
549 549 def flag_cspc(self):
550 550
551 551 if self.data_cspc is None:
552 552 return True
553 553
554 554 return False
555 555
556 556 @property
557 557 def flag_dc(self):
558 558
559 559 if self.data_dc is None:
560 560 return True
561 561
562 562 return False
563 563
564 564 @property
565 565 def timeInterval(self):
566 566
567 567 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
568 568 if self.nmodes:
569 569 return self.nmodes*timeInterval
570 570 else:
571 571 return timeInterval
572 572
573 573 def getPower(self):
574 574
575 575 factor = self.normFactor
576 576 z = self.data_spc / factor
577 577 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
578 578 avg = numpy.average(z, axis=1)
579 579 return 10 * numpy.log10(avg)
580 580
581 581 def getCoherence(self, pairsList=None, phase=False):
582 582
583 583 z = []
584 584 if pairsList is None:
585 585 pairsIndexList = self.pairsIndexList
586 586 else:
587 587 pairsIndexList = []
588 588 for pair in pairsList:
589 589 if pair not in self.pairsList:
590 590 raise ValueError("Pair %s is not in dataOut.pairsList" % (
591 591 pair))
592 592 pairsIndexList.append(self.pairsList.index(pair))
593 593 for i in range(len(pairsIndexList)):
594 594 pair = self.pairsList[pairsIndexList[i]]
595 595 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
596 596 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
597 597 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
598 598 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
599 599 if phase:
600 600 data = numpy.arctan2(avgcoherenceComplex.imag,
601 601 avgcoherenceComplex.real) * 180 / numpy.pi
602 602 else:
603 603 data = numpy.abs(avgcoherenceComplex)
604 604
605 605 z.append(data)
606 606
607 607 return numpy.array(z)
608 608
609 609 def setValue(self, value):
610 610
611 611 print("This property should not be initialized")
612 612
613 613 return
614 614
615 615 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
616 616
617 617
618 618 class SpectraHeis(Spectra):
619 619
620 620 def __init__(self):
621 621
622 622 self.radarControllerHeaderObj = RadarControllerHeader()
623 623 self.systemHeaderObj = SystemHeader()
624 624 self.type = "SpectraHeis"
625 625 self.nProfiles = None
626 626 self.heightList = None
627 627 self.channelList = None
628 628 self.flagNoData = True
629 629 self.flagDiscontinuousBlock = False
630 630 self.utctime = None
631 631 self.blocksize = None
632 632 self.profileIndex = 0
633 633 self.nCohInt = 1
634 634 self.nIncohInt = 1
635 635
636 636 @property
637 637 def normFactor(self):
638 638 pwcode = 1
639 639 if self.flagDecodeData:
640 640 pwcode = numpy.sum(self.code[0]**2)
641 641
642 642 normFactor = self.nIncohInt * self.nCohInt * pwcode
643 643
644 644 return normFactor
645 645
646 646 @property
647 647 def timeInterval(self):
648 648
649 649 return self.ippSeconds * self.nCohInt * self.nIncohInt
650 650
651 651
652 652 class Fits(JROData):
653 653
654 654 def __init__(self):
655 655
656 656 self.type = "Fits"
657 657 self.nProfiles = None
658 658 self.heightList = None
659 659 self.channelList = None
660 660 self.flagNoData = True
661 661 self.utctime = None
662 662 self.nCohInt = 1
663 663 self.nIncohInt = 1
664 664 self.useLocalTime = True
665 665 self.profileIndex = 0
666 666 self.timeZone = 0
667 667
668 668 def getTimeRange(self):
669 669
670 670 datatime = []
671 671
672 672 datatime.append(self.ltctime)
673 673 datatime.append(self.ltctime + self.timeInterval)
674 674
675 675 datatime = numpy.array(datatime)
676 676
677 677 return datatime
678 678
679 679 def getChannelIndexList(self):
680 680
681 681 return list(range(self.nChannels))
682 682
683 683 def getNoise(self, type=1):
684 684
685 685
686 686 if type == 1:
687 687 noise = self.getNoisebyHildebrand()
688 688
689 689 if type == 2:
690 690 noise = self.getNoisebySort()
691 691
692 692 if type == 3:
693 693 noise = self.getNoisebyWindow()
694 694
695 695 return noise
696 696
697 697 @property
698 698 def timeInterval(self):
699 699
700 700 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
701 701
702 702 return timeInterval
703 703
704 704 @property
705 705 def ippSeconds(self):
706 706 '''
707 707 '''
708 708 return self.ipp_sec
709 709
710 710 noise = property(getNoise, "I'm the 'nHeights' property.")
711 711
712 712
713 713 class Correlation(JROData):
714 714
715 715 def __init__(self):
716 716 '''
717 717 Constructor
718 718 '''
719 719 self.radarControllerHeaderObj = RadarControllerHeader()
720 720 self.systemHeaderObj = SystemHeader()
721 721 self.type = "Correlation"
722 722 self.data = None
723 723 self.dtype = None
724 724 self.nProfiles = None
725 725 self.heightList = None
726 726 self.channelList = None
727 727 self.flagNoData = True
728 728 self.flagDiscontinuousBlock = False
729 729 self.utctime = None
730 730 self.timeZone = 0
731 731 self.dstFlag = None
732 732 self.errorCount = None
733 733 self.blocksize = None
734 734 self.flagDecodeData = False # asumo q la data no esta decodificada
735 735 self.flagDeflipData = False # asumo q la data no esta sin flip
736 736 self.pairsList = None
737 737 self.nPoints = None
738 738
739 739 def getPairsList(self):
740 740
741 741 return self.pairsList
742 742
743 743 def getNoise(self, mode=2):
744 744
745 745 indR = numpy.where(self.lagR == 0)[0][0]
746 746 indT = numpy.where(self.lagT == 0)[0][0]
747 747
748 748 jspectra0 = self.data_corr[:, :, indR, :]
749 749 jspectra = copy.copy(jspectra0)
750 750
751 751 num_chan = jspectra.shape[0]
752 752 num_hei = jspectra.shape[2]
753 753
754 754 freq_dc = jspectra.shape[1] / 2
755 755 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
756 756
757 757 if ind_vel[0] < 0:
758 758 ind_vel[list(range(0, 1))] = ind_vel[list(
759 759 range(0, 1))] + self.num_prof
760 760
761 761 if mode == 1:
762 762 jspectra[:, freq_dc, :] = (
763 763 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
764 764
765 765 if mode == 2:
766 766
767 767 vel = numpy.array([-2, -1, 1, 2])
768 768 xx = numpy.zeros([4, 4])
769 769
770 770 for fil in range(4):
771 771 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
772 772
773 773 xx_inv = numpy.linalg.inv(xx)
774 774 xx_aux = xx_inv[0, :]
775 775
776 776 for ich in range(num_chan):
777 777 yy = jspectra[ich, ind_vel, :]
778 778 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
779 779
780 780 junkid = jspectra[ich, freq_dc, :] <= 0
781 781 cjunkid = sum(junkid)
782 782
783 783 if cjunkid.any():
784 784 jspectra[ich, freq_dc, junkid.nonzero()] = (
785 785 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
786 786
787 787 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
788 788
789 789 return noise
790 790
791 791 @property
792 792 def timeInterval(self):
793 793
794 794 return self.ippSeconds * self.nCohInt * self.nProfiles
795 795
796 796 def splitFunctions(self):
797 797
798 798 pairsList = self.pairsList
799 799 ccf_pairs = []
800 800 acf_pairs = []
801 801 ccf_ind = []
802 802 acf_ind = []
803 803 for l in range(len(pairsList)):
804 804 chan0 = pairsList[l][0]
805 805 chan1 = pairsList[l][1]
806 806
807 807 # Obteniendo pares de Autocorrelacion
808 808 if chan0 == chan1:
809 809 acf_pairs.append(chan0)
810 810 acf_ind.append(l)
811 811 else:
812 812 ccf_pairs.append(pairsList[l])
813 813 ccf_ind.append(l)
814 814
815 815 data_acf = self.data_cf[acf_ind]
816 816 data_ccf = self.data_cf[ccf_ind]
817 817
818 818 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
819 819
820 820 @property
821 821 def normFactor(self):
822 822 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
823 823 acf_pairs = numpy.array(acf_pairs)
824 824 normFactor = numpy.zeros((self.nPairs, self.nHeights))
825 825
826 826 for p in range(self.nPairs):
827 827 pair = self.pairsList[p]
828 828
829 829 ch0 = pair[0]
830 830 ch1 = pair[1]
831 831
832 832 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
833 833 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
834 834 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
835 835
836 836 return normFactor
837 837
838 838
839 839 class Parameters(Spectra):
840 840
841 841 groupList = None # List of Pairs, Groups, etc
842 842 data_param = None # Parameters obtained
843 843 data_pre = None # Data Pre Parametrization
844 844 data_SNR = None # Signal to Noise Ratio
845 845 abscissaList = None # Abscissa, can be velocities, lags or time
846 846 utctimeInit = None # Initial UTC time
847 847 paramInterval = None # Time interval to calculate Parameters in seconds
848 848 useLocalTime = True
849 849 # Fitting
850 850 data_error = None # Error of the estimation
851 851 constants = None
852 852 library = None
853 853 # Output signal
854 854 outputInterval = None # Time interval to calculate output signal in seconds
855 855 data_output = None # Out signal
856 856 nAvg = None
857 857 noise_estimation = None
858 858 GauSPC = None # Fit gaussian SPC
859 859
860 860 def __init__(self):
861 861 '''
862 862 Constructor
863 863 '''
864 864 self.radarControllerHeaderObj = RadarControllerHeader()
865 865 self.systemHeaderObj = SystemHeader()
866 866 self.type = "Parameters"
867 867 self.timeZone = 0
868 868
869 869 def getTimeRange1(self, interval):
870 870
871 871 datatime = []
872 872
873 873 if self.useLocalTime:
874 874 time1 = self.utctimeInit - self.timeZone * 60
875 875 else:
876 876 time1 = self.utctimeInit
877 877
878 878 datatime.append(time1)
879 879 datatime.append(time1 + interval)
880 880 datatime = numpy.array(datatime)
881 881
882 882 return datatime
883 883
884 884 @property
885 885 def timeInterval(self):
886 886
887 887 if hasattr(self, 'timeInterval1'):
888 888 return self.timeInterval1
889 889 else:
890 890 return self.paramInterval
891 891
892 892 def setValue(self, value):
893 893
894 894 print("This property should not be initialized")
895 895
896 896 return
897 897
898 898 def getNoise(self):
899 899
900 900 return self.spc_noise
901 901
902 902 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
903 903
904 904
905 905 class PlotterData(object):
906 906 '''
907 907 Object to hold data to be plotted
908 908 '''
909 909
910 MAXNUMX = 200
911 MAXNUMY = 200
910 MAXNUMX = 1000
911 MAXNUMY = 1000
912 912
913 913 def __init__(self, code, exp_code, localtime=True):
914 914
915 915 self.key = code
916 916 self.exp_code = exp_code
917 917 self.ready = False
918 918 self.flagNoData = False
919 919 self.localtime = localtime
920 920 self.data = {}
921 921 self.meta = {}
922 922 self.__heights = []
923 923
924 924 def __str__(self):
925 925 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
926 926 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
927 927
928 928 def __len__(self):
929 929 return len(self.data)
930 930
931 931 def __getitem__(self, key):
932 932 if isinstance(key, int):
933 933 return self.data[self.times[key]]
934 934 elif isinstance(key, str):
935 935 ret = numpy.array([self.data[x][key] for x in self.times])
936 936 if ret.ndim > 1:
937 937 ret = numpy.swapaxes(ret, 0, 1)
938 938 return ret
939 939
940 940 def __contains__(self, key):
941 941 return key in self.data[self.min_time]
942 942
943 943 def setup(self):
944 944 '''
945 945 Configure object
946 946 '''
947 947 self.type = ''
948 948 self.ready = False
949 949 del self.data
950 950 self.data = {}
951 951 self.__heights = []
952 952 self.__all_heights = set()
953 953
954 954 def shape(self, key):
955 955 '''
956 956 Get the shape of the one-element data for the given key
957 957 '''
958 958
959 959 if len(self.data[self.min_time][key]):
960 960 return self.data[self.min_time][key].shape
961 961 return (0,)
962 962
963 963 def update(self, data, tm, meta={}):
964 964 '''
965 965 Update data object with new dataOut
966 966 '''
967 967
968 968 self.data[tm] = data
969 969
970 970 for key, value in meta.items():
971 971 setattr(self, key, value)
972 972
973 973 def normalize_heights(self):
974 974 '''
975 975 Ensure same-dimension of the data for different heighList
976 976 '''
977 977
978 978 H = numpy.array(list(self.__all_heights))
979 979 H.sort()
980 980 for key in self.data:
981 981 shape = self.shape(key)[:-1] + H.shape
982 982 for tm, obj in list(self.data[key].items()):
983 983 h = self.__heights[self.times.tolist().index(tm)]
984 984 if H.size == h.size:
985 985 continue
986 986 index = numpy.where(numpy.in1d(H, h))[0]
987 987 dummy = numpy.zeros(shape) + numpy.nan
988 988 if len(shape) == 2:
989 989 dummy[:, index] = obj
990 990 else:
991 991 dummy[index] = obj
992 992 self.data[key][tm] = dummy
993 993
994 994 self.__heights = [H for tm in self.times]
995 995
996 def jsonify(self, tm, plot_name, plot_type, decimate=False):
996 def jsonify(self, tm, plot_name, plot_type, key=None, decimate=False):
997 997 '''
998 998 Convert data to json
999 999 '''
1000 1000
1001 if key is None:
1002 key = self.key
1003
1001 1004 meta = {}
1002 1005 meta['xrange'] = []
1003 1006 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1004 tmp = self.data[tm][self.key]
1007 tmp = self.data[tm][key]
1005 1008 shape = tmp.shape
1006 1009 if len(shape) == 2:
1007 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1010 data = self.roundFloats(self.data[tm][key][::, ::dy].tolist())
1008 1011 elif len(shape) == 3:
1009 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1012 dx = int(self.data[tm][key].shape[1]/self.MAXNUMX) + 1
1010 1013 data = self.roundFloats(
1011 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1014 self.data[tm][key][::, ::dx, ::dy].tolist())
1012 1015 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1013 1016 else:
1014 data = self.roundFloats(self.data[tm][self.key].tolist())
1017 data = self.roundFloats(self.data[tm][key].tolist())
1015 1018
1016 1019 ret = {
1017 1020 'plot': plot_name,
1018 1021 'code': self.exp_code,
1019 1022 'time': float(tm),
1020 1023 'data': data,
1021 1024 }
1022 1025 meta['type'] = plot_type
1023 1026 meta['interval'] = float(self.interval)
1024 1027 meta['localtime'] = self.localtime
1025 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1028 #meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1029 meta['yrange'] = self.roundFloats(self.lat[::dy].tolist())
1030 meta['xrange'] = self.roundFloats(self.lon[::dy].tolist())
1026 1031 meta.update(self.meta)
1027 1032 ret['metadata'] = meta
1028 1033 return json.dumps(ret)
1029 1034
1030 1035 @property
1031 1036 def times(self):
1032 1037 '''
1033 1038 Return the list of times of the current data
1034 1039 '''
1035 1040
1036 1041 ret = [t for t in self.data]
1037 1042 ret.sort()
1038 1043 return numpy.array(ret)
1039 1044
1040 1045 @property
1041 1046 def min_time(self):
1042 1047 '''
1043 1048 Return the minimun time value
1044 1049 '''
1045 1050
1046 1051 return self.times[0]
1047 1052
1048 1053 @property
1049 1054 def max_time(self):
1050 1055 '''
1051 1056 Return the maximun time value
1052 1057 '''
1053 1058
1054 1059 return self.times[-1]
1055 1060
1056 1061 # @property
1057 1062 # def heights(self):
1058 1063 # '''
1059 1064 # Return the list of heights of the current data
1060 1065 # '''
1061 1066
1062 1067 # return numpy.array(self.__heights[-1])
1063 1068
1064 1069 @staticmethod
1065 1070 def roundFloats(obj):
1066 1071 if isinstance(obj, list):
1067 1072 return list(map(PlotterData.roundFloats, obj))
1068 1073 elif isinstance(obj, float):
1069 return round(obj, 2)
1074 return round(obj, 4)
Show file before | Show file after | Hide comments
@@ -1,732 +1,732
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Base class to create plot operations
6 6
7 7 """
8 8
9 9 import os
10 10 import sys
11 11 import zmq
12 12 import time
13 13 import numpy
14 14 import datetime
15 15 from collections import deque
16 16 from functools import wraps
17 17 from threading import Thread
18 18 import matplotlib,re
19 19
20 20 if 'BACKEND' in os.environ:
21 21 matplotlib.use(os.environ['BACKEND'])
22 22 elif 'linux' in sys.platform:
23 23 matplotlib.use("TkAgg")
24 24 elif 'darwin' in sys.platform:
25 25 matplotlib.use('MacOSX')
26 26 else:
27 27 from schainpy.utils import log
28 28 log.warning('Using default Backend="Agg"', 'INFO')
29 29 matplotlib.use('Agg')
30 30
31 31 import matplotlib.pyplot as plt
32 32 from matplotlib.patches import Polygon
33 33 from mpl_toolkits.axes_grid1 import make_axes_locatable
34 34 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
35 35
36 36 from .plotting_codes import *
37 37
38 38 from schainpy.model.data.jrodata import PlotterData
39 39 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
40 40 from schainpy.utils import log
41 41
42 42 for name, cb_table in sophy_cb_tables:
43 43 ncmap = matplotlib.colors.ListedColormap(cb_table, name=name)
44 44 matplotlib.pyplot.register_cmap(cmap=ncmap)
45 45
46 46 EARTH_RADIUS = 6.3710e3
47 47
48 48 def ll2xy(lat1, lon1, lat2, lon2):
49 49
50 50 p = 0.017453292519943295
51 51 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
52 52 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
53 53 r = 12742 * numpy.arcsin(numpy.sqrt(a))
54 54 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
55 55 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
56 56 theta = -theta + numpy.pi/2
57 57 return r*numpy.cos(theta), r*numpy.sin(theta)
58 58
59 59
60 60 def km2deg(km):
61 61 '''
62 62 Convert distance in km to degrees
63 63 '''
64 64
65 65 return numpy.rad2deg(km/EARTH_RADIUS)
66 66
67 67
68 68 def figpause(interval):
69 69 backend = plt.rcParams['backend']
70 70 if backend in matplotlib.rcsetup.interactive_bk:
71 71 figManager = matplotlib._pylab_helpers.Gcf.get_active()
72 72 if figManager is not None:
73 73 canvas = figManager.canvas
74 74 if canvas.figure.stale:
75 75 canvas.draw()
76 76 try:
77 77 canvas.start_event_loop(interval)
78 78 except:
79 79 pass
80 80 return
81 81
82 82 def popup(message):
83 83 '''
84 84 '''
85 85
86 86 fig = plt.figure(figsize=(12, 8), facecolor='r')
87 87 text = '\n'.join([s.strip() for s in message.split(':')])
88 88 fig.text(0.01, 0.5, text, ha='left', va='center',
89 89 size='20', weight='heavy', color='w')
90 90 fig.show()
91 91 figpause(1000)
92 92
93 93
94 94 class Throttle(object):
95 95 '''
96 96 Decorator that prevents a function from being called more than once every
97 97 time period.
98 98 To create a function that cannot be called more than once a minute, but
99 99 will sleep until it can be called:
100 100 @Throttle(minutes=1)
101 101 def foo():
102 102 pass
103 103
104 104 for i in range(10):
105 105 foo()
106 106 print "This function has run %s times." % i
107 107 '''
108 108
109 109 def __init__(self, seconds=0, minutes=0, hours=0):
110 110 self.throttle_period = datetime.timedelta(
111 111 seconds=seconds, minutes=minutes, hours=hours
112 112 )
113 113
114 114 self.time_of_last_call = datetime.datetime.min
115 115
116 116 def __call__(self, fn):
117 117 @wraps(fn)
118 118 def wrapper(*args, **kwargs):
119 119 coerce = kwargs.pop('coerce', None)
120 120 if coerce:
121 121 self.time_of_last_call = datetime.datetime.now()
122 122 return fn(*args, **kwargs)
123 123 else:
124 124 now = datetime.datetime.now()
125 125 time_since_last_call = now - self.time_of_last_call
126 126 time_left = self.throttle_period - time_since_last_call
127 127
128 128 if time_left > datetime.timedelta(seconds=0):
129 129 return
130 130
131 131 self.time_of_last_call = datetime.datetime.now()
132 132 return fn(*args, **kwargs)
133 133
134 134 return wrapper
135 135
136 136 def apply_throttle(value):
137 137
138 138 @Throttle(seconds=value)
139 139 def fnThrottled(fn):
140 140 fn()
141 141
142 142 return fnThrottled
143 143
144 144
145 145 @MPDecorator
146 146 class Plot(Operation):
147 147 """Base class for Schain plotting operations
148 148
149 149 This class should never be use directtly you must subclass a new operation,
150 150 children classes must be defined as follow:
151 151
152 152 ExamplePlot(Plot):
153 153
154 154 CODE = 'code'
155 155 colormap = 'jet'
156 156 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
157 157
158 158 def setup(self):
159 159 pass
160 160
161 161 def plot(self):
162 162 pass
163 163
164 164 """
165 165
166 166 CODE = 'Figure'
167 167 colormap = 'jet'
168 168 bgcolor = 'white'
169 169 buffering = True
170 170 __missing = 1E30
171 171
172 172 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
173 173 'showprofile']
174 174
175 175 def __init__(self):
176 176
177 177 Operation.__init__(self)
178 178 self.isConfig = False
179 179 self.isPlotConfig = False
180 180 self.save_time = 0
181 181 self.sender_time = 0
182 182 self.data = None
183 183 self.firsttime = True
184 184 self.sender_queue = deque(maxlen=10)
185 185 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
186 186
187 187 def __fmtTime(self, x, pos):
188 188 '''
189 189 '''
190 190
191 191 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
192 192
193 193 def __setup(self, **kwargs):
194 194 '''
195 195 Initialize variables
196 196 '''
197 197
198 198 self.figures = []
199 199 self.axes = []
200 200 self.cb_axes = []
201 201 self.localtime = kwargs.pop('localtime', True)
202 202 self.show = kwargs.get('show', True)
203 203 self.save = kwargs.get('save', False)
204 204 self.save_period = kwargs.get('save_period', 0)
205 205 self.colormap = kwargs.get('colormap', self.colormap)
206 206 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
207 207 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
208 208 self.colormaps = kwargs.get('colormaps', None)
209 209 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
210 210 self.showprofile = kwargs.get('showprofile', False)
211 211 self.title = kwargs.get('wintitle', self.CODE.upper())
212 212 self.cb_label = kwargs.get('cb_label', None)
213 213 self.cb_labels = kwargs.get('cb_labels', None)
214 214 self.labels = kwargs.get('labels', None)
215 215 self.xaxis = kwargs.get('xaxis', 'frequency')
216 216 self.zmin = kwargs.get('zmin', None)
217 217 self.zmax = kwargs.get('zmax', None)
218 218 self.zlimits = kwargs.get('zlimits', None)
219 219 self.xmin = kwargs.get('xmin', None)
220 220 self.xmax = kwargs.get('xmax', None)
221 221 self.xrange = kwargs.get('xrange', 12)
222 222 self.xscale = kwargs.get('xscale', None)
223 223 self.ymin = kwargs.get('ymin', None)
224 224 self.ymax = kwargs.get('ymax', None)
225 225 self.yscale = kwargs.get('yscale', None)
226 226 self.xlabel = kwargs.get('xlabel', None)
227 227 self.attr_time = kwargs.get('attr_time', 'utctime')
228 228 self.attr_data = kwargs.get('attr_data', 'data_param')
229 229 self.decimation = kwargs.get('decimation', None)
230 230 self.oneFigure = kwargs.get('oneFigure', True)
231 231 self.width = kwargs.get('width', None)
232 232 self.height = kwargs.get('height', None)
233 233 self.colorbar = kwargs.get('colorbar', True)
234 234 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
235 235 self.channels = kwargs.get('channels', None)
236 236 self.titles = kwargs.get('titles', [])
237 237 self.polar = False
238 238 self.type = kwargs.get('type', 'iq')
239 239 self.grid = kwargs.get('grid', False)
240 240 self.pause = kwargs.get('pause', False)
241 241 self.save_code = kwargs.get('save_code', self.CODE)
242 242 self.throttle = kwargs.get('throttle', 0)
243 243 self.exp_code = kwargs.get('exp_code', None)
244 244 self.server = kwargs.get('server', False)
245 245 self.sender_period = kwargs.get('sender_period', 60)
246 246 self.tag = kwargs.get('tag', '')
247 247 self.height_index = kwargs.get('height_index', None)
248 248 self.__throttle_plot = apply_throttle(self.throttle)
249 249 code = self.attr_data if self.attr_data else self.CODE
250 250 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
251 251 self.ang_min = kwargs.get('ang_min', None)
252 252 self.ang_max = kwargs.get('ang_max', None)
253 253 self.mode = kwargs.get('mode', None)
254 254 self.snr_threshold = kwargs.get('snr_threshold', 0)
255 255
256 256
257 257 if self.server:
258 258 if not self.server.startswith('tcp://'):
259 259 self.server = 'tcp://{}'.format(self.server)
260 260 log.success(
261 261 'Sending to server: {}'.format(self.server),
262 262 self.name
263 263 )
264 264
265 265 if isinstance(self.attr_data, str):
266 266 self.attr_data = [self.attr_data]
267 267
268 268 def __setup_plot(self):
269 269 '''
270 270 Common setup for all figures, here figures and axes are created
271 271 '''
272 272
273 273 self.setup()
274 274
275 275 self.time_label = 'LT' if self.localtime else 'UTC'
276 276
277 277 if self.width is None:
278 278 self.width = 8
279 279
280 280 self.figures = []
281 281 self.axes = []
282 282 self.cb_axes = []
283 283 self.pf_axes = []
284 284 self.cmaps = []
285 285
286 286 size = '15%' if self.ncols == 1 else '30%'
287 287 pad = '4%' if self.ncols == 1 else '8%'
288 288
289 289 if self.oneFigure:
290 290 if self.height is None:
291 291 self.height = 1.4 * self.nrows + 1
292 292 fig = plt.figure(figsize=(self.width, self.height),
293 293 edgecolor='k',
294 294 facecolor='w')
295 295 self.figures.append(fig)
296 296 for n in range(self.nplots):
297 297 ax = fig.add_subplot(self.nrows, self.ncols,
298 298 n + 1, polar=self.polar)
299 299 ax.tick_params(labelsize=8)
300 300 ax.firsttime = True
301 301 ax.index = 0
302 302 ax.press = None
303 303 self.axes.append(ax)
304 304 if self.showprofile:
305 305 cax = self.__add_axes(ax, size=size, pad=pad)
306 306 cax.tick_params(labelsize=8)
307 307 self.pf_axes.append(cax)
308 308 else:
309 309 if self.height is None:
310 310 self.height = 3
311 311 for n in range(self.nplots):
312 312 fig = plt.figure(figsize=(self.width, self.height),
313 313 edgecolor='k',
314 314 facecolor='w')
315 315 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
316 316 ax.tick_params(labelsize=8)
317 317 ax.firsttime = True
318 318 ax.index = 0
319 319 ax.press = None
320 320 self.figures.append(fig)
321 321 self.axes.append(ax)
322 322 if self.showprofile:
323 323 cax = self.__add_axes(ax, size=size, pad=pad)
324 324 cax.tick_params(labelsize=8)
325 325 self.pf_axes.append(cax)
326 326
327 327 for n in range(self.nrows):
328 328 if self.colormaps is not None:
329 329 cmap = plt.get_cmap(self.colormaps[n])
330 330 else:
331 331 cmap = plt.get_cmap(self.colormap)
332 332 cmap.set_bad(self.bgcolor, 1.)
333 333 self.cmaps.append(cmap)
334 334
335 335 def __add_axes(self, ax, size='30%', pad='8%'):
336 336 '''
337 337 Add new axes to the given figure
338 338 '''
339 339 divider = make_axes_locatable(ax)
340 340 nax = divider.new_horizontal(size=size, pad=pad)
341 341 ax.figure.add_axes(nax)
342 342 return nax
343 343
344 344 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
345 345 '''
346 346 Create a masked array for missing data
347 347 '''
348 348 if x_buffer.shape[0] < 2:
349 349 return x_buffer, y_buffer, z_buffer
350 350
351 351 deltas = x_buffer[1:] - x_buffer[0:-1]
352 352 x_median = numpy.median(deltas)
353 353
354 354 index = numpy.where(deltas > 5 * x_median)
355 355
356 356 if len(index[0]) != 0:
357 357 z_buffer[::, index[0], ::] = self.__missing
358 358 z_buffer = numpy.ma.masked_inside(z_buffer,
359 359 0.99 * self.__missing,
360 360 1.01 * self.__missing)
361 361
362 362 return x_buffer, y_buffer, z_buffer
363 363
364 364 def decimate(self):
365 365
366 366 # dx = int(len(self.x)/self.__MAXNUMX) + 1
367 367 dy = int(len(self.y) / self.decimation) + 1
368 368
369 369 # x = self.x[::dx]
370 370 x = self.x
371 371 y = self.y[::dy]
372 372 z = self.z[::, ::, ::dy]
373 373
374 374 return x, y, z
375 375
376 376 def format(self):
377 377 '''
378 378 Set min and max values, labels, ticks and titles
379 379 '''
380 380
381 381 for n, ax in enumerate(self.axes):
382 382 if ax.firsttime:
383 383 if self.xaxis != 'time':
384 384 xmin = self.xmin
385 385 xmax = self.xmax
386 386 else:
387 387 xmin = self.tmin
388 388 xmax = self.tmin + self.xrange*60*60
389 389 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
390 390 ax.xaxis.set_major_locator(LinearLocator(9))
391 391 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
392 392 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
393 393 ax.set_facecolor(self.bgcolor)
394 394 if self.xscale:
395 395 ax.xaxis.set_major_formatter(FuncFormatter(
396 396 lambda x, pos: '{0:g}'.format(x*self.xscale)))
397 397 if self.yscale:
398 398 ax.yaxis.set_major_formatter(FuncFormatter(
399 399 lambda x, pos: '{0:g}'.format(x*self.yscale)))
400 400 if self.xlabel is not None:
401 401 ax.set_xlabel(self.xlabel)
402 402 if self.ylabel is not None:
403 403 ax.set_ylabel(self.ylabel)
404 404 if self.showprofile:
405 405 self.pf_axes[n].set_ylim(ymin, ymax)
406 406 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
407 407 self.pf_axes[n].set_xlabel('dB')
408 408 self.pf_axes[n].grid(b=True, axis='x')
409 409 [tick.set_visible(False)
410 410 for tick in self.pf_axes[n].get_yticklabels()]
411 411 if self.colorbar:
412 412 ax.cbar = plt.colorbar(
413 413 ax.plt, ax=ax, fraction=0.05, pad=0.06, aspect=10)
414 414 ax.cbar.ax.tick_params(labelsize=8)
415 415 ax.cbar.ax.press = None
416 416 if self.cb_label:
417 417 ax.cbar.set_label(self.cb_label, size=8)
418 418 elif self.cb_labels:
419 419 ax.cbar.set_label(self.cb_labels[n], size=8)
420 420 else:
421 421 ax.cbar = None
422 422 ax.set_xlim(xmin, xmax)
423 423 ax.set_ylim(ymin, ymax)
424 424 ax.firsttime = False
425 425 if self.grid:
426 426 ax.grid(True)
427 427 if not self.polar:
428 428 ax.set_title('{} {} {}'.format(
429 429 self.titles[n],
430 430 self.getDateTime(self.data.max_time).strftime(
431 431 '%Y-%m-%d %H:%M:%S'),
432 432 self.time_label),
433 433 size=8)
434 434 else:
435 435 #ax.set_title('{}'.format(self.titles[n]), size=8)
436 436 ax.set_title('{} {} {}'.format(
437 437 self.titles[n],
438 438 self.getDateTime(self.data.max_time).strftime(
439 439 '%Y-%m-%d %H:%M:%S'),
440 440 self.time_label),
441 441 size=8)
442 442 ax.set_ylim(0, self.ymax)
443 443 ax.set_yticks(ax.get_yticks(), labels=ax.get_yticks(), color='white')
444 444 ax.yaxis.labelpad = 28
445 445
446 446 if self.firsttime:
447 447 for n, fig in enumerate(self.figures):
448 448 fig.subplots_adjust(**self.plots_adjust)
449 449 self.firsttime = False
450 450
451 451 def clear_figures(self):
452 452 '''
453 453 Reset axes for redraw plots
454 454 '''
455 455
456 456 for ax in self.axes+self.pf_axes+self.cb_axes:
457 457 ax.clear()
458 458 ax.firsttime = True
459 459 if hasattr(ax, 'cbar') and ax.cbar:
460 460 ax.cbar.remove()
461 461
462 462 def __plot(self):
463 463 '''
464 464 Main function to plot, format and save figures
465 465 '''
466 466
467 467 self.plot()
468 468 self.format()
469 469
470 470 for n, fig in enumerate(self.figures):
471 471 if self.nrows == 0 or self.nplots == 0:
472 472 log.warning('No data', self.name)
473 473 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
474 474 fig.canvas.manager.set_window_title(self.CODE)
475 475 continue
476 476
477 477 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
478 478 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
479 479 fig.canvas.draw()
480 480 if self.show:
481 481 fig.show()
482 482 figpause(0.01)
483 483
484 484 if self.save:
485 485 if self.CODE=="PPI" or self.CODE=="RHI":
486 486 self.save_figure(n,stitle =self.titles)
487 487 else:
488 488 self.save_figure(n)
489 489
490 490 if self.server:
491 491 self.send_to_server()
492 492
493 493 def __update(self, dataOut, timestamp):
494 494 '''
495 495 '''
496 496
497 497 metadata = {
498 498 'yrange': dataOut.heightList,
499 499 'interval': dataOut.timeInterval,
500 500 'channels': dataOut.channelList
501 501 }
502 502
503 503 data, meta = self.update(dataOut)
504 504 metadata.update(meta)
505 505 self.data.update(data, timestamp, metadata)
506 506
507 507 def save_figure(self, n,stitle=None):
508 508 '''
509 509 '''
510 510 if stitle is not None:
511 511 s_string = re.sub(r"[^A-Z0-9.]","",str(stitle))
512 512 new_string=s_string[:3]+"_"+s_string[4:6]+"_"+s_string[6:]
513 513
514 514 if self.oneFigure:
515 515 if (self.data.max_time - self.save_time) <= self.save_period:
516 516 return
517 517
518 518 self.save_time = self.data.max_time
519 519
520 520 fig = self.figures[n]
521 521
522 522 if self.throttle == 0:
523 523 if self.oneFigure:
524 524 if stitle is not None:
525 525 figname = os.path.join(
526 526 self.save,
527 527 self.save_code + '_' + new_string,
528 528 '{}_{}_{}.png'.format(
529 529 self.save_code,
530 530 new_string,
531 531 self.getDateTime(self.data.max_time).strftime(
532 532 '%Y%m%d_%H%M%S',
533 533 ),
534 534 )
535 535 )
536 536 else:
537 537 figname = os.path.join(
538 538 self.save,
539 539 self.save_code,
540 540 '{}_{}.png'.format(
541 541 self.save_code,
542 542 self.getDateTime(self.data.max_time).strftime(
543 543 '%Y%m%d_%H%M%S'
544 544 ),
545 545 )
546 546 )
547 547 else:
548 548 figname = os.path.join(
549 549 self.save,
550 550 self.save_code,
551 551 '{}_ch{}_{}.png'.format(
552 552 self.save_code,n,
553 553 self.getDateTime(self.data.max_time).strftime(
554 554 '%Y%m%d_%H%M%S'
555 555 ),
556 556 )
557 557 )
558 558 log.log('Saving figure: {}'.format(figname), self.name)
559 559 if not os.path.isdir(os.path.dirname(figname)):
560 560 os.makedirs(os.path.dirname(figname))
561 561 fig.savefig(figname)
562 562
563 563 figname = os.path.join(
564 564 self.save,
565 565 '{}_{}.png'.format(
566 566 self.save_code,
567 567 self.getDateTime(self.data.min_time).strftime(
568 568 '%Y%m%d'
569 569 ),
570 570 )
571 571 )
572 572
573 573 log.log('Saving figure: {}'.format(figname), self.name)
574 574 if not os.path.isdir(os.path.dirname(figname)):
575 575 os.makedirs(os.path.dirname(figname))
576 576 fig.savefig(figname)
577 577
578 578 def send_to_server(self):
579 579 '''
580 580 '''
581 581
582 582 if self.exp_code == None:
583 583 log.warning('Missing `exp_code` skipping sending to server...')
584 584
585 585 last_time = self.data.max_time
586 586 interval = last_time - self.sender_time
587 587 if interval < self.sender_period:
588 588 return
589 589
590 590 self.sender_time = last_time
591 591
592 592 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
593 593 for attr in attrs:
594 594 value = getattr(self, attr)
595 595 if value:
596 596 if isinstance(value, (numpy.float32, numpy.float64)):
597 597 value = round(float(value), 2)
598 598 self.data.meta[attr] = value
599 if self.colormap == 'jet':
599 if self.colormap == 'jet' or self.colormap == 'sophy_w':
600 600 self.data.meta['colormap'] = 'Jet'
601 elif 'RdBu' in self.colormap:
601 elif 'sophy_v' in self.colormap:
602 602 self.data.meta['colormap'] = 'RdBu'
603 603 else:
604 604 self.data.meta['colormap'] = 'Viridis'
605 605 self.data.meta['interval'] = int(interval)
606 606
607 607 self.sender_queue.append(last_time)
608 608
609 609 while True:
610 610 try:
611 611 tm = self.sender_queue.popleft()
612 612 except IndexError:
613 613 break
614 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
614 msg = self.data.jsonify(tm, self.save_code, self.plot_type, key='var')
615 615 self.socket.send_string(msg)
616 616 socks = dict(self.poll.poll(2000))
617 617 if socks.get(self.socket) == zmq.POLLIN:
618 618 reply = self.socket.recv_string()
619 619 if reply == 'ok':
620 620 log.log("Response from server ok", self.name)
621 621 time.sleep(0.1)
622 622 continue
623 623 else:
624 624 log.warning(
625 625 "Malformed reply from server: {}".format(reply), self.name)
626 626 else:
627 627 log.warning(
628 628 "No response from server, retrying...", self.name)
629 629 self.sender_queue.appendleft(tm)
630 630 self.socket.setsockopt(zmq.LINGER, 0)
631 631 self.socket.close()
632 632 self.poll.unregister(self.socket)
633 633 self.socket = self.context.socket(zmq.REQ)
634 634 self.socket.connect(self.server)
635 635 self.poll.register(self.socket, zmq.POLLIN)
636 636 break
637 637
638 638 def setup(self):
639 639 '''
640 640 This method should be implemented in the child class, the following
641 641 attributes should be set:
642 642
643 643 self.nrows: number of rows
644 644 self.ncols: number of cols
645 645 self.nplots: number of plots (channels or pairs)
646 646 self.ylabel: label for Y axes
647 647 self.titles: list of axes title
648 648
649 649 '''
650 650 raise NotImplementedError
651 651
652 652 def plot(self):
653 653 '''
654 654 Must be defined in the child class, the actual plotting method
655 655 '''
656 656 raise NotImplementedError
657 657
658 658 def update(self, dataOut):
659 659 '''
660 660 Must be defined in the child class, update self.data with new data
661 661 '''
662 662
663 663 data = {
664 664 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
665 665 }
666 666 meta = {}
667 667
668 668 return data, meta
669 669
670 670 def run(self, dataOut, **kwargs):
671 671 '''
672 672 Main plotting routine
673 673 '''
674 674
675 675 if self.isConfig is False:
676 676 self.__setup(**kwargs)
677 677
678 678 if self.localtime:
679 679 self.getDateTime = datetime.datetime.fromtimestamp
680 680 else:
681 681 self.getDateTime = datetime.datetime.utcfromtimestamp
682 682
683 683 self.data.setup()
684 684 self.isConfig = True
685 685 if self.server:
686 686 self.context = zmq.Context()
687 687 self.socket = self.context.socket(zmq.REQ)
688 688 self.socket.connect(self.server)
689 689 self.poll = zmq.Poller()
690 690 self.poll.register(self.socket, zmq.POLLIN)
691 691
692 692 tm = getattr(dataOut, self.attr_time)
693 693
694 694 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
695 695 self.save_time = tm
696 696 self.__plot()
697 697 self.tmin += self.xrange*60*60
698 698 self.data.setup()
699 699 self.clear_figures()
700 700
701 701 self.__update(dataOut, tm)
702 702
703 703 if self.isPlotConfig is False:
704 704 self.__setup_plot()
705 705 self.isPlotConfig = True
706 706 if self.xaxis == 'time':
707 707 dt = self.getDateTime(tm)
708 708 if self.xmin is None:
709 709 self.tmin = tm
710 710 self.xmin = dt.hour
711 711 minutes = (self.xmin-int(self.xmin)) * 60
712 712 seconds = (minutes - int(minutes)) * 60
713 713 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
714 714 datetime.datetime(1970, 1, 1)).total_seconds()
715 715 if self.localtime:
716 716 self.tmin += time.timezone
717 717
718 718 if self.xmin is not None and self.xmax is not None:
719 719 self.xrange = self.xmax - self.xmin
720 720
721 721 if self.throttle == 0:
722 722 self.__plot()
723 723 else:
724 724 self.__throttle_plot(self.__plot)#, coerce=coerce)
725 725
726 726 def close(self):
727 727
728 728 if self.data and not self.data.flagNoData:
729 729 self.save_time = 0
730 730 self.__plot()
731 731 if self.data and not self.data.flagNoData and self.pause:
732 732 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,512 +1,519
1 1 import os
2 2 import datetime
3 3 import numpy
4 4 from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
5 5
6 6 from schainpy.model.graphics.jroplot_base import Plot, plt
7 7 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
8 8 from schainpy.utils import log
9 # libreria wradlib
10 #import wradlib as wrl
9
10 import wradlib.georef as georef
11 11
12 12 EARTH_RADIUS = 6.3710e3
13 13
14 14
15 15 def ll2xy(lat1, lon1, lat2, lon2):
16 16
17 17 p = 0.017453292519943295
18 18 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
19 19 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
20 20 r = 12742 * numpy.arcsin(numpy.sqrt(a))
21 21 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
22 22 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
23 23 theta = -theta + numpy.pi/2
24 24 return r*numpy.cos(theta), r*numpy.sin(theta)
25 25
26 26
27 27 def km2deg(km):
28 28 '''
29 29 Convert distance in km to degrees
30 30 '''
31 31
32 32 return numpy.rad2deg(km/EARTH_RADIUS)
33 33
34 34
35 35
36 36 class SpectralMomentsPlot(SpectraPlot):
37 37 '''
38 38 Plot for Spectral Moments
39 39 '''
40 40 CODE = 'spc_moments'
41 41 # colormap = 'jet'
42 42 # plot_type = 'pcolor'
43 43
44 44 class DobleGaussianPlot(SpectraPlot):
45 45 '''
46 46 Plot for Double Gaussian Plot
47 47 '''
48 48 CODE = 'gaussian_fit'
49 49 # colormap = 'jet'
50 50 # plot_type = 'pcolor'
51 51
52 52 class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
53 53 '''
54 54 Plot SpectraCut with Double Gaussian Fit
55 55 '''
56 56 CODE = 'cut_gaussian_fit'
57 57
58 58 class SnrPlot(RTIPlot):
59 59 '''
60 60 Plot for SNR Data
61 61 '''
62 62
63 63 CODE = 'snr'
64 64 colormap = 'jet'
65 65
66 66 def update(self, dataOut):
67 67
68 68 data = {
69 69 'snr': 10*numpy.log10(dataOut.data_snr)
70 70 }
71 71
72 72 return data, {}
73 73
74 74 class DopplerPlot(RTIPlot):
75 75 '''
76 76 Plot for DOPPLER Data (1st moment)
77 77 '''
78 78
79 79 CODE = 'dop'
80 80 colormap = 'jet'
81 81
82 82 def update(self, dataOut):
83 83
84 84 data = {
85 85 'dop': 10*numpy.log10(dataOut.data_dop)
86 86 }
87 87
88 88 return data, {}
89 89
90 90 class PowerPlot(RTIPlot):
91 91 '''
92 92 Plot for Power Data (0 moment)
93 93 '''
94 94
95 95 CODE = 'pow'
96 96 colormap = 'jet'
97 97
98 98 def update(self, dataOut):
99 99 data = {
100 100 'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
101 101 }
102 102 return data, {}
103 103
104 104 class SpectralWidthPlot(RTIPlot):
105 105 '''
106 106 Plot for Spectral Width Data (2nd moment)
107 107 '''
108 108
109 109 CODE = 'width'
110 110 colormap = 'jet'
111 111
112 112 def update(self, dataOut):
113 113
114 114 data = {
115 115 'width': dataOut.data_width
116 116 }
117 117
118 118 return data, {}
119 119
120 120 class SkyMapPlot(Plot):
121 121 '''
122 122 Plot for meteors detection data
123 123 '''
124 124
125 125 CODE = 'param'
126 126
127 127 def setup(self):
128 128
129 129 self.ncols = 1
130 130 self.nrows = 1
131 131 self.width = 7.2
132 132 self.height = 7.2
133 133 self.nplots = 1
134 134 self.xlabel = 'Zonal Zenith Angle (deg)'
135 135 self.ylabel = 'Meridional Zenith Angle (deg)'
136 136 self.polar = True
137 137 self.ymin = -180
138 138 self.ymax = 180
139 139 self.colorbar = False
140 140
141 141 def plot(self):
142 142
143 143 arrayParameters = numpy.concatenate(self.data['param'])
144 144 error = arrayParameters[:, -1]
145 145 indValid = numpy.where(error == 0)[0]
146 146 finalMeteor = arrayParameters[indValid, :]
147 147 finalAzimuth = finalMeteor[:, 3]
148 148 finalZenith = finalMeteor[:, 4]
149 149
150 150 x = finalAzimuth * numpy.pi / 180
151 151 y = finalZenith
152 152
153 153 ax = self.axes[0]
154 154
155 155 if ax.firsttime:
156 156 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
157 157 else:
158 158 ax.plot.set_data(x, y)
159 159
160 160 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
161 161 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
162 162 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
163 163 dt2,
164 164 len(x))
165 165 self.titles[0] = title
166 166
167 167
168 168 class GenericRTIPlot(Plot):
169 169 '''
170 170 Plot for data_xxxx object
171 171 '''
172 172
173 173 CODE = 'param'
174 174 colormap = 'viridis'
175 175 plot_type = 'pcolorbuffer'
176 176
177 177 def setup(self):
178 178 self.xaxis = 'time'
179 179 self.ncols = 1
180 180 self.nrows = self.data.shape('param')[0]
181 181 self.nplots = self.nrows
182 182 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
183 183
184 184 if not self.xlabel:
185 185 self.xlabel = 'Time'
186 186
187 187 self.ylabel = 'Range [km]'
188 188 if not self.titles:
189 189 self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
190 190
191 191 def update(self, dataOut):
192 192
193 193 data = {
194 194 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
195 195 }
196 196
197 197 meta = {}
198 198
199 199 return data, meta
200 200
201 201 def plot(self):
202 202 # self.data.normalize_heights()
203 203 self.x = self.data.times
204 204 self.y = self.data.yrange
205 205 self.z = self.data['param']
206 206 self.z = 10*numpy.log10(self.z)
207 207 self.z = numpy.ma.masked_invalid(self.z)
208 208
209 209 if self.decimation is None:
210 210 x, y, z = self.fill_gaps(self.x, self.y, self.z)
211 211 else:
212 212 x, y, z = self.fill_gaps(*self.decimate())
213 213
214 214 for n, ax in enumerate(self.axes):
215 215
216 216 self.zmax = self.zmax if self.zmax is not None else numpy.max(
217 217 self.z[n])
218 218 self.zmin = self.zmin if self.zmin is not None else numpy.min(
219 219 self.z[n])
220 220
221 221 if ax.firsttime:
222 222 if self.zlimits is not None:
223 223 self.zmin, self.zmax = self.zlimits[n]
224 224
225 225 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
226 226 vmin=self.zmin,
227 227 vmax=self.zmax,
228 228 cmap=self.cmaps[n]
229 229 )
230 230 else:
231 231 if self.zlimits is not None:
232 232 self.zmin, self.zmax = self.zlimits[n]
233 233 ax.collections.remove(ax.collections[0])
234 234 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
235 235 vmin=self.zmin,
236 236 vmax=self.zmax,
237 237 cmap=self.cmaps[n]
238 238 )
239 239
240 240
241 241 class PolarMapPlot(Plot):
242 242 '''
243 243 Plot for weather radar
244 244 '''
245 245
246 246 CODE = 'param'
247 247 colormap = 'seismic'
248 248
249 249 def setup(self):
250 250 self.ncols = 1
251 251 self.nrows = 1
252 252 self.width = 9
253 253 self.height = 8
254 254 self.mode = self.data.meta['mode']
255 255 if self.channels is not None:
256 256 self.nplots = len(self.channels)
257 257 self.nrows = len(self.channels)
258 258 else:
259 259 self.nplots = self.data.shape(self.CODE)[0]
260 260 self.nrows = self.nplots
261 261 self.channels = list(range(self.nplots))
262 262 if self.mode == 'E':
263 263 self.xlabel = 'Longitude'
264 264 self.ylabel = 'Latitude'
265 265 else:
266 266 self.xlabel = 'Range (km)'
267 267 self.ylabel = 'Height (km)'
268 268 self.bgcolor = 'white'
269 269 self.cb_labels = self.data.meta['units']
270 270 self.lat = self.data.meta['latitude']
271 271 self.lon = self.data.meta['longitude']
272 272 self.xmin, self.xmax = float(
273 273 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
274 274 self.ymin, self.ymax = float(
275 275 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
276 276 # self.polar = True
277 277
278 278 def plot(self):
279 279
280 280 for n, ax in enumerate(self.axes):
281 281 data = self.data['param'][self.channels[n]]
282 282
283 283 zeniths = numpy.linspace(
284 284 0, self.data.meta['max_range'], data.shape[1])
285 285 if self.mode == 'E':
286 286 azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
287 287 r, theta = numpy.meshgrid(zeniths, azimuths)
288 288 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
289 289 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
290 290 x = km2deg(x) + self.lon
291 291 y = km2deg(y) + self.lat
292 292 else:
293 293 azimuths = numpy.radians(self.data.yrange)
294 294 r, theta = numpy.meshgrid(zeniths, azimuths)
295 295 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
296 296 self.y = zeniths
297 297
298 298 if ax.firsttime:
299 299 if self.zlimits is not None:
300 300 self.zmin, self.zmax = self.zlimits[n]
301 301 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
302 302 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
303 303 vmin=self.zmin,
304 304 vmax=self.zmax,
305 305 cmap=self.cmaps[n])
306 306 else:
307 307 if self.zlimits is not None:
308 308 self.zmin, self.zmax = self.zlimits[n]
309 309 ax.collections.remove(ax.collections[0])
310 310 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
311 311 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
312 312 vmin=self.zmin,
313 313 vmax=self.zmax,
314 314 cmap=self.cmaps[n])
315 315
316 316 if self.mode == 'A':
317 317 continue
318 318
319 319 # plot district names
320 320 f = open('/data/workspace/schain_scripts/distrito.csv')
321 321 for line in f:
322 322 label, lon, lat = [s.strip() for s in line.split(',') if s]
323 323 lat = float(lat)
324 324 lon = float(lon)
325 325 # ax.plot(lon, lat, '.b', ms=2)
326 326 ax.text(lon, lat, label.decode('utf8'), ha='center',
327 327 va='bottom', size='8', color='black')
328 328
329 329 # plot limites
330 330 limites = []
331 331 tmp = []
332 332 for line in open('/data/workspace/schain_scripts/lima.csv'):
333 333 if '#' in line:
334 334 if tmp:
335 335 limites.append(tmp)
336 336 tmp = []
337 337 continue
338 338 values = line.strip().split(',')
339 339 tmp.append((float(values[0]), float(values[1])))
340 340 for points in limites:
341 341 ax.add_patch(
342 342 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
343 343
344 344 # plot Cuencas
345 345 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
346 346 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
347 347 values = [line.strip().split(',') for line in f]
348 348 points = [(float(s[0]), float(s[1])) for s in values]
349 349 ax.add_patch(Polygon(points, ec='b', fc='none'))
350 350
351 351 # plot grid
352 352 for r in (15, 30, 45, 60):
353 353 ax.add_artist(plt.Circle((self.lon, self.lat),
354 354 km2deg(r), color='0.6', fill=False, lw=0.2))
355 355 ax.text(
356 356 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
357 357 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
358 358 '{}km'.format(r),
359 359 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
360 360
361 361 if self.mode == 'E':
362 362 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
363 363 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
364 364 else:
365 365 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
366 366 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
367 367
368 368 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
369 369 self.titles = ['{} {}'.format(
370 370 self.data.parameters[x], title) for x in self.channels]
371 371
372 372 class WeatherParamsPlot(Plot):
373 373 #CODE = 'RHI'
374 374 #plot_name = 'RHI'
375 #plot_type = 'rhistyle'
375 plot_type = 'scattermap'
376 376 buffering = False
377 377
378 378 def setup(self):
379 379
380 380 self.ncols = 1
381 381 self.nrows = 1
382 382 self.nplots= 1
383 383 self.ylabel= 'Range [km]'
384 384 self.xlabel= 'Range [km]'
385 385 self.polar = True
386 386 self.grid = True
387 387 if self.channels is not None:
388 388 self.nplots = len(self.channels)
389 389 self.nrows = len(self.channels)
390 390 else:
391 391 self.nplots = self.data.shape(self.CODE)[0]
392 392 self.nrows = self.nplots
393 393 self.channels = list(range(self.nplots))
394 394
395 395 self.colorbar=True
396 396 self.width =8
397 397 self.height =8
398 398 self.ini =0
399 399 self.len_azi =0
400 400 self.buffer_ini = None
401 401 self.buffer_ele = None
402 402 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
403 403 self.flag =0
404 404 self.indicador= 0
405 405 self.last_data_ele = None
406 406 self.val_mean = None
407 407
408 408 def update(self, dataOut):
409 409
410 410 data = {}
411 411 meta = {}
412 412 if hasattr(dataOut, 'dataPP_POWER'):
413 413 factor = 1
414 414 if hasattr(dataOut, 'nFFTPoints'):
415 415 factor = dataOut.normFactor
416 416
417 417 mask = dataOut.data_snr<self.snr_threshold
418 418
419 419 if 'pow' in self.attr_data[0].lower():
420 420 # data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
421 data['data'] = numpy.ma.masked_array(10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor)), mask=mask)
421 tmp = numpy.ma.masked_array(10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor)), mask=mask)
422 422 else:
423 data['data'] = numpy.ma.masked_array(getattr(dataOut, self.attr_data[0]), mask=mask)
424 # data['data'] = getattr(dataOut, self.attr_data[0])
423 tmp = numpy.ma.masked_array(getattr(dataOut, self.attr_data[0]), mask=mask)
424 # tmp = getattr(dataOut, self.attr_data[0])
425
426 r = dataOut.heightList
427 delta_height = r[1]-r[0]
428 valid = numpy.where(r>=0)[0]
429 data['r'] = numpy.arange(len(valid))*delta_height
430
431 try:
432 data['data'] = tmp[self.channels[0]][:,valid]
433 except:
434 data['data'] = tmp[0][:,valid]
425 435
426 436 if dataOut.mode_op == 'PPI':
427 437 self.CODE = 'PPI'
428 438 self.title = self.CODE
429 439 elif dataOut.mode_op == 'RHI':
430 440 self.CODE = 'RHI'
431 441 self.title = self.CODE
432 442
433 443 data['azi'] = dataOut.data_azi
434 444 data['ele'] = dataOut.data_ele
435 445 data['mode_op'] = dataOut.mode_op
446 var = data['data'].flatten()
447 r = numpy.tile(data['r'], data['data'].shape[0]).reshape(data['data'].shape)*1000
448 lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))
449 meta['lat'] = lla[:,:,1].flatten()[var.mask==False]
450 meta['lon'] = lla[:,:,0].flatten()[var.mask==False]
451 data['var'] = numpy.array([var[var.mask==False]])
436 452
437 453 return data, meta
438 454
439 455 def plot(self):
440 456 data = self.data[-1]
441 r = self.data.yrange
442 delta_height = r[1]-r[0]
443 r_mask = numpy.where(r>=0)[0]
444 r = numpy.arange(len(r_mask))*delta_height
445 self.y = 2*r
446
447 try:
448 z = data['data'][self.channels[0]][:,r_mask]
449 except:
450 z = data['data'][0][:,r_mask]
451
457 z = data['data']
458 r = data['r']
452 459 self.titles = []
453 460
454 461 self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
455 462 self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
456 463 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
457 464 self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)
458 465
459 466 if data['mode_op'] == 'RHI':
460 467 try:
461 468 if self.data['mode_op'][-2] == 'PPI':
462 469 self.ang_min = None
463 470 self.ang_max = None
464 471 except:
465 472 pass
466 473 self.ang_min = self.ang_min if self.ang_min else 0
467 474 self.ang_max = self.ang_max if self.ang_max else 90
468 475 r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
469 476 elif data['mode_op'] == 'PPI':
470 477 try:
471 478 if self.data['mode_op'][-2] == 'RHI':
472 479 self.ang_min = None
473 480 self.ang_max = None
474 481 except:
475 482 pass
476 483 self.ang_min = self.ang_min if self.ang_min else 0
477 484 self.ang_max = self.ang_max if self.ang_max else 360
478 485 r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
479 486
480 487 self.clear_figures()
481 488
482 489 for i,ax in enumerate(self.axes):
483 490
484 491 if ax.firsttime:
485 492 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
486 493 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
487 494 if data['mode_op'] == 'PPI':
488 495 ax.set_theta_direction(-1)
489 496 ax.set_theta_offset(numpy.pi/2)
490 497
491 498 else:
492 499 ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
493 500 ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
494 501 if data['mode_op'] == 'PPI':
495 502 ax.set_theta_direction(-1)
496 503 ax.set_theta_offset(numpy.pi/2)
497 504
498 505 ax.grid(True)
499 506 if data['mode_op'] == 'RHI':
500 507 len_aux = int(data['azi'].shape[0]/4)
501 508 mean = numpy.mean(data['azi'][len_aux:-len_aux])
502 509 if len(self.channels) !=1:
503 510 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
504 511 else:
505 512 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
506 513 elif data['mode_op'] == 'PPI':
507 514 len_aux = int(data['ele'].shape[0]/4)
508 515 mean = numpy.mean(data['ele'][len_aux:-len_aux])
509 516 if len(self.channels) !=1:
510 517 self.titles = ['PPI {} at EL: {} CH {}'.format(self.self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
511 518 else:
512 519 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
General Comments 0
You need to be logged in to leave comments. Login now