##// END OF EJS Templates
schain
RSS
changes from v4
joabAM -
r1664:a9d92c6e95c5
parent child
Show More
Side by Side Unified
Context file:
r1664:a9d92c6e95c5 Loading diff...:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1420 lines changed) Show them Hide them
@@ -0,0 +1,1420
1 """
2 The module ASTRO_COORDS.py gathers classes and functions for coordinates transformation. Additiona-
3 lly a class EquatorialCorrections and celestial bodies are defined. The first of these is to correct
4 any error in the location of the body and the second to know the location of certain celestial bo-
5 dies in the sky.
6
7 MODULES CALLED:
8 OS, NUMPY, NUMERIC, SCIPY, TIME_CONVERSIONS
9
10 MODIFICATION HISTORY:
11 Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Sep 20, 2009.
12 """
13
14 import numpy
15 #import Numeric
16 import scipy.interpolate
17 import os
18 import sys
19 from schainpy.model.utils import TimeTools
20 from schainpy.model.utils import Misc_Routines
21
22 class EquatorialCorrections():
23 def __init__(self):
24 """
25 EquatorialCorrections class creates an object to call methods to correct the loca-
26 tion of the celestial bodies.
27
28 Modification History
29 --------------------
30 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009.
31 """
32
33 pass
34
35 def co_nutate(self,jd,ra,dec):
36 """
37 co_nutate calculates changes in RA and Dec due to nutation of the Earth's rotation
38 Additionally it returns the obliquity of the ecliptic (eps), nutation in the longi-
39 tude of the ecliptic (d_psi) and nutation in the pbliquity of the ecliptic (d_eps).
40
41 Parameters
42 ----------
43 jd = Julian Date (Scalar or array).
44 RA = A scalar o array giving the Right Ascention of interest.
45 Dec = A scalar o array giving the Right Ascention of interest.
46
47 Return
48 ------
49 d_ra = Correction to ra due to nutation.
50 d_dec = Correction to dec due to nutation.
51
52 Examples
53 --------
54 >> Julian = 2462088.7
55 >> Ra = 41.547213
56 >> Dec = 49.348483
57 >> [d_ra,d_dec,eps,d_psi,d_eps] = co_nutate(julian,Ra,Dec)
58 >> print d_ra, d_dec, eps, d_psi, d_eps
59 [ 15.84276651] [ 6.21641029] [ 0.4090404] [ 14.85990198] [ 2.70408658]
60
61 Modification history
62 --------------------
63 Written by Chris O'Dell, 2002.
64 Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
65 """
66
67 jd = numpy.atleast_1d(jd)
68 ra = numpy.atleast_1d(ra)
69 dec = numpy.atleast_1d(dec)
70
71 # Useful transformation constants
72 d2as = numpy.pi/(180.*3600.)
73
74 # Julian centuries from J2000 of jd
75 T = (jd - 2451545.0)/36525.0
76
77 # Must calculate obliquity of ecliptic
78 [d_psi, d_eps] = self.nutate(jd)
79 d_psi = numpy.atleast_1d(d_psi)
80 d_eps = numpy.atleast_1d(d_eps)
81
82 eps0 = (23.4392911*3600.) - (46.8150*T) - (0.00059*T**2) + (0.001813*T**3)
83 # True obliquity of the ecliptic in radians
84 eps = (eps0 + d_eps)/3600.*Misc_Routines.CoFactors.d2r
85
86 # Useful numbers
87 ce = numpy.cos(eps)
88 se = numpy.sin(eps)
89
90 # Convert Ra-Dec to equatorial rectangular coordinates
91 x = numpy.cos(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r)
92 y = numpy.sin(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r)
93 z = numpy.sin(dec*Misc_Routines.CoFactors.d2r)
94
95 # Apply corrections to each rectangular coordinate
96 x2 = x - (y*ce + z*se)*d_psi*Misc_Routines.CoFactors.s2r
97 y2 = y + (x*ce*d_psi - z*d_eps)*Misc_Routines.CoFactors.s2r
98 z2 = z + (x*se*d_psi + y*d_eps)*Misc_Routines.CoFactors.s2r
99
100 # Convert bask to equatorial spherical coordinates
101 r = numpy.sqrt(x2**2. + y2**2. + z2**2.)
102 xyproj =numpy.sqrt(x2**2. + y2**2.)
103
104 ra2 = x2*0.0
105 dec2 = x2*0.0
106
107 xyproj = numpy.atleast_1d(xyproj)
108 z = numpy.atleast_1d(z)
109 r = numpy.atleast_1d(r)
110 x2 = numpy.atleast_1d(x2)
111 y2 = numpy.atleast_1d(y2)
112 z2 = numpy.atleast_1d(z2)
113 ra2 = numpy.atleast_1d(ra2)
114 dec2 = numpy.atleast_1d(dec2)
115
116 w1 = numpy.where((xyproj==0) & (z!=0))
117 w2 = numpy.where(xyproj!=0)
118
119 # Calculate Ra and Dec in radians (later convert to degrees)
120 if w1[0].size>0:
121 # Places where xyproj=0 (point at NCP or SCP)
122 dec2[w1] = numpy.arcsin(z2[w1]/r[w1])
123 ra2[w1] = 0
124
125 if w2[0].size>0:
126 # Places other than NCP or SCP
127 ra2[w2] = numpy.arctan2(y2[w2],x2[w2])
128 dec2[w2] = numpy.arcsin(z2[w2]/r[w2])
129
130 # Converting to degree
131 ra2 = ra2/Misc_Routines.CoFactors.d2r
132 dec2 = dec2/Misc_Routines.CoFactors.d2r
133
134 w = numpy.where(ra2<0.)
135 if w[0].size>0:
136 ra2[w] = ra2[w] + 360.
137
138 # Return changes in Ra and Dec in arcseconds
139 d_ra = (ra2 -ra)*3600.
140 d_dec = (dec2 - dec)*3600.
141
142 return d_ra, d_dec, eps, d_psi, d_eps
143
144 def nutate(self,jd):
145 """
146 nutate returns the nutation in longitude and obliquity for a given Julian date.
147
148 Parameters
149 ----------
150 jd = Julian ephemeris date, scalar or vector.
151
152 Return
153 ------
154 nut_long = The nutation in longitude.
155 nut_obliq = The nutation in latitude.
156
157 Example
158 -------
159 >> julian = 2446895.5
160 >> [nut_long,nut_obliq] = nutate(julian)
161 >> print nut_long, nut_obliq
162 -3.78793107711 9.44252069864
163
164 >> julians = 2415020.5 + numpy.arange(50)
165 >> [nut_long,nut_obliq] = nutate(julians)
166
167 Modification History
168 --------------------
169 Written by W.Landsman (Goddard/HSTX), June 1996.
170 Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
171 """
172
173 jd = numpy.atleast_1d(jd)
174
175 # Form time in Julian centuries from 1900
176 t = (jd - 2451545.0)/36525.0
177
178 # Mean elongation of the moon
179 coeff1 = numpy.array([1/189474.0,-0.0019142,445267.111480,297.85036])
180 d = numpy.poly1d(coeff1)
181 d = d(t)*Misc_Routines.CoFactors.d2r
182 d = self.cirrange(d,rad=1)
183
184 # Sun's mean elongation
185 coeff2 = numpy.array([-1./3e5,-0.0001603,35999.050340,357.52772])
186 m = numpy.poly1d(coeff2)
187 m = m(t)*Misc_Routines.CoFactors.d2r
188 m = self.cirrange(m,rad=1)
189
190 # Moon's mean elongation
191 coeff3 = numpy.array([1.0/5.625e4,0.0086972,477198.867398,134.96298])
192 mprime = numpy.poly1d(coeff3)
193 mprime = mprime(t)*Misc_Routines.CoFactors.d2r
194 mprime = self.cirrange(mprime,rad=1)
195
196 # Moon's argument of latitude
197 coeff4 = numpy.array([-1.0/3.27270e5,-0.0036825,483202.017538,93.27191])
198 f = numpy.poly1d(coeff4)
199 f = f(t)*Misc_Routines.CoFactors.d2r
200 f = self.cirrange(f,rad=1)
201
202 # Longitude fo the ascending node of the Moon's mean orbit on the ecliptic, measu-
203 # red from the mean equinox of the date.
204 coeff5 = numpy.array([1.0/4.5e5,0.0020708,-1934.136261,125.04452])
205 omega = numpy.poly1d(coeff5)
206 omega = omega(t)*Misc_Routines.CoFactors.d2r
207 omega = self.cirrange(omega,rad=1)
208
209 d_lng = numpy.array([0,-2,0,0,0,0,-2,0,0,-2,-2,-2,0,2,0,2,0,0,-2,0,2,0,0,-2,0,-2,0,0,\
210 2,-2,0,-2,0,0,2,2,0,-2,0,2,2,-2,-2,2,2,0,-2,-2,0,-2,-2,0,-1,-2,1,0,0,-1,0,\
211 0,2,0,2])
212
213 m_lng = numpy.array([0,0,0,0,1,0,1,0,0,-1])
214 m_lng = numpy.append(m_lng,numpy.zeros(17))
215 m_lng = numpy.append(m_lng,numpy.array([2,0,2,1,0,-1,0,0,0,1,1,-1,0,0,0,0,0,0,-1,-1,0,0,\
216 0,1,0,0,1,0,0,0,-1,1,-1,-1,0,-1]))
217
218 mp_lng = numpy.array([0,0,0,0,0,1,0,0,1,0,1,0,-1,0,1,-1,-1,1,2,-2,0,2,2,1,0,0, -1, 0,\
219 -1,0,0,1,0,2,-1,1,0,1,0,0,1,2,1,-2,0,1,0,0,2,2,0,1,1,0,0,1,-2,1,1,1,-1,3,0])
220
221 f_lng = numpy.array([0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,2,0,2,0,2,2,2,0,2,2,2,2,0,0,2,0,0,\
222 0,-2,2,2,2,0,2,2,0,2,2,0,0,0,2,0,2,0,2,-2,0,0,0,2,2,0,0,2,2,2,2])
223
224 om_lng = numpy.array([1,2,2,2,0,0,2,1,2,2,0,1,2,0,1,2,1,1,0,1,2,2,0,2,0,0,1,0,1,2,1, \
225 1,1,0,1,2,2,0,2,1,0,2,1,1,1,0,1,1,1,1,1,0,0,0,0,0,2,0,0,2,2,2,2])
226
227 sin_lng = numpy.array([-171996,-13187,-2274,2062,1426,712,-517,-386,-301, 217, -158, \
228 129,123,63,63,-59,-58,-51,48,46,-38,-31,29,29,26,-22,21,17,16,-16,-15,-13,\
229 -12,11,-10,-8,7,-7,-7,-7,6,6,6,-6,-6,5,-5,-5,-5,4,4,4,-4,-4,-4,3,-3,-3,-3,\
230 -3,-3,-3,-3])
231
232 sdelt = numpy.array([-174.2,-1.6,-0.2,0.2,-3.4,0.1,1.2,-0.4,0,-0.5,0, 0.1, 0, 0, 0.1,\
233 0,-0.1])
234 sdelt = numpy.append(sdelt,numpy.zeros(10))
235 sdelt = numpy.append(sdelt,numpy.array([-0.1, 0, 0.1]))
236 sdelt = numpy.append(sdelt,numpy.zeros(33))
237
238 cos_lng = numpy.array([92025,5736,977,-895,54,-7,224,200,129,-95,0,-70,-53,0,-33,26, \
239 32,27,0,-24,16,13,0,-12,0,0,-10,0,-8,7,9,7,6,0,5,3,-3,0,3,3,0,-3,-3,3,3,0,\
240 3,3,3])
241 cos_lng = numpy.append(cos_lng,numpy.zeros(14))
242
243 cdelt = numpy.array([8.9,-3.1,-0.5,0.5,-0.1,0.0,-0.6,0.0,-0.1,0.3])
244 cdelt = numpy.append(cdelt,numpy.zeros(53))
245
246 # Sum the periodic terms.
247 n = numpy.size(jd)
248 nut_long = numpy.zeros(n)
249 nut_obliq = numpy.zeros(n)
250
251 d_lng = d_lng.reshape(numpy.size(d_lng),1)
252 d = d.reshape(numpy.size(d),1)
253 matrix_d_lng = numpy.dot(d_lng,d.transpose())
254
255 m_lng = m_lng.reshape(numpy.size(m_lng),1)
256 m = m.reshape(numpy.size(m),1)
257 matrix_m_lng = numpy.dot(m_lng,m.transpose())
258
259 mp_lng = mp_lng.reshape(numpy.size(mp_lng),1)
260 mprime = mprime.reshape(numpy.size(mprime),1)
261 matrix_mp_lng = numpy.dot(mp_lng,mprime.transpose())
262
263 f_lng = f_lng.reshape(numpy.size(f_lng),1)
264 f = f.reshape(numpy.size(f),1)
265 matrix_f_lng = numpy.dot(f_lng,f.transpose())
266
267 om_lng = om_lng.reshape(numpy.size(om_lng),1)
268 omega = omega.reshape(numpy.size(omega),1)
269 matrix_om_lng = numpy.dot(om_lng,omega.transpose())
270
271 arg = matrix_d_lng + matrix_m_lng + matrix_mp_lng + matrix_f_lng + matrix_om_lng
272
273 sarg = numpy.sin(arg)
274 carg = numpy.cos(arg)
275
276 for ii in numpy.arange(n):
277 nut_long[ii] = 0.0001*numpy.sum((sdelt*t[ii] + sin_lng)*sarg[:,ii])
278 nut_obliq[ii] = 0.0001*numpy.sum((cdelt*t[ii] + cos_lng)*carg[:,ii])
279
280 if numpy.size(jd)==1:
281 nut_long = nut_long[0]
282 nut_obliq = nut_obliq[0]
283
284 return nut_long, nut_obliq
285
286 def co_aberration(self,jd,ra,dec):
287 """
288 co_aberration calculates changes to Ra and Dec due to "the effect of aberration".
289
290 Parameters
291 ----------
292 jd = Julian Date (Scalar or vector).
293 ra = A scalar o vector giving the Right Ascention of interest.
294 dec = A scalar o vector giving the Declination of interest.
295
296 Return
297 ------
298 d_ra = The correction to right ascension due to aberration (must be added to ra to
299 get the correct value).
300 d_dec = The correction to declination due to aberration (must be added to the dec
301 to get the correct value).
302 eps = True obliquity of the ecliptic (in radians).
303
304 Examples
305 --------
306 >> Julian = 2462088.7
307 >> Ra = 41.547213
308 >> Dec = 49.348483
309 >> [d_ra,d_dec,eps] = co_aberration(julian,Ra,Dec)
310 >> print d_ra, d_dec, eps
311 [ 30.04441796] [ 6.69837858] [ 0.40904059]
312
313 Modification history
314 --------------------
315 Written by Chris O'Dell , Univ. of Wisconsin, June 2002.
316 Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
317 """
318
319 # Julian centuries from J2000 of jd.
320 T = (jd - 2451545.0)/36525.0
321
322 # Getting obliquity of ecliptic
323 njd = numpy.size(jd)
324 jd = numpy.atleast_1d(jd)
325 ra = numpy.atleast_1d(ra)
326 dec = numpy.atleast_1d(dec)
327
328 d_psi = numpy.zeros(njd)
329 d_epsilon = d_psi
330 for ii in numpy.arange(njd):
331 [dp,de] = self.nutate(jd[ii])
332 d_psi[ii] = dp
333 d_epsilon[ii] = de
334
335 coeff = 23 + 26/60. + 21.488/3600.
336 eps0 = coeff*3600. - 46.8150*T - 0.00059*T**2. + 0.001813*T**3.
337 # True obliquity of the ecliptic in radians
338 eps = (eps0 + d_epsilon)/3600*Misc_Routines.CoFactors.d2r
339
340 celestialbodies = CelestialBodies()
341 [sunra,sundec,sunlon,sunobliq] = celestialbodies.sunpos(jd)
342
343 # Earth's orbital eccentricity
344 e = 0.016708634 - 0.000042037*T - 0.0000001267*T**2.
345
346 # longitude of perihelion, in degrees
347 pi = 102.93735 + 1.71946*T + 0.00046*T**2
348
349 # Constant of aberration, in arcseconds
350 k = 20.49552
351
352 cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r) ; sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r)
353 ce = numpy.cos(eps) ; te = numpy.tan(eps)
354 cp = numpy.cos(pi*Misc_Routines.CoFactors.d2r) ; sp = numpy.sin(pi*Misc_Routines.CoFactors.d2r)
355 cs = numpy.cos(sunlon*Misc_Routines.CoFactors.d2r) ; ss = numpy.sin(sunlon*Misc_Routines.CoFactors.d2r)
356 ca = numpy.cos(ra*Misc_Routines.CoFactors.d2r) ; sa = numpy.sin(ra*Misc_Routines.CoFactors.d2r)
357
358 term1 = (ca*cs*ce + sa*ss)/cd
359 term2 = (ca*cp*ce + sa*sp)/cd
360 term3 = (cs*ce*(te*cd - sa*sd) + ca*sd*ss)
361 term4 = (cp*ce*(te*cd - sa*sd) + ca*sd*sp)
362
363 d_ra = -k*term1 + e*k*term2
364 d_dec = -k*term3 + e*k*term4
365
366 return d_ra, d_dec, eps
367
368 def precess(self,ra,dec,equinox1=None,equinox2=None,FK4=0,rad=0):
369 """
370 precess coordinates from EQUINOX1 to EQUINOX2
371
372 Parameters
373 -----------
374 ra = A scalar o vector giving the Right Ascention of interest.
375 dec = A scalar o vector giving the Declination of interest.
376 equinox1 = Original equinox of coordinates, numeric scalar. If omitted, the __Pre-
377 cess will query for equinox1 and equinox2.
378 equinox2 = Original equinox of coordinates.
379 FK4 = If this keyword is set and non-zero, the FK4 (B1950) system will be used
380 otherwise FK5 (J2000) will be used instead.
381 rad = If this keyword is set and non-zero, then the input and output RAD and DEC
382 vectors are in radian rather than degree.
383
384 Return
385 ------
386 ra = Right ascension after precession (scalar or vector) in degrees, unless the rad
387 keyword is set.
388 dec = Declination after precession (scalar or vector) in degrees, unless the rad
389 keyword is set.
390
391 Examples
392 --------
393 >> Ra = 329.88772
394 >> Dec = -56.992515
395 >> [p_ra,p_dec] = precess(Ra,Dec,1950,1975,FK4=1)
396 >> print p_ra, p_dec
397 [ 330.31442971] [-56.87186154]
398
399 Modification history
400 --------------------
401 Written by Wayne Landsman, STI Corporation, August 1986.
402 Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
403 """
404
405 npts = numpy.size(ra)
406 ra = numpy.atleast_1d(ra)
407 dec = numpy.atleast_1d(dec)
408
409 if rad==0:
410 ra_rad = ra*Misc_Routines.CoFactors.d2r
411 dec_rad = dec*Misc_Routines.CoFactors.d2r
412 else:
413 ra_rad = ra
414 dec_rad = dec
415
416 x = numpy.zeros((npts,3))
417 x[:,0] = numpy.cos(dec_rad)*numpy.cos(ra_rad)
418 x[:,1] = numpy.cos(dec_rad)*numpy.sin(ra_rad)
419 x[:,2] = numpy.sin(dec_rad)
420
421 # Use premat function to get precession matrix from equinox1 to equinox2
422 r = self.premat(equinox1,equinox2,FK4)
423
424 x2 = numpy.dot(r,x.transpose())
425
426 ra_rad = numpy.arctan2(x2[1,:],x2[0,:])
427 dec_rad = numpy.arcsin(x2[2,:])
428
429 if rad==0:
430 ra = ra_rad/Misc_Routines.CoFactors.d2r
431 ra = ra + (ra<0)*360.
432 dec = dec_rad/Misc_Routines.CoFactors.d2r
433 else:
434 ra = ra_rad
435 ra = ra + (ra<0)*numpy.pi*2.
436 dec = dec_rad
437
438 return ra, dec
439
440 def premat(self,equinox1,equinox2,FK4=0):
441 """
442 premat returns the precession matrix needed to go from EQUINOX1 to EQUINOX2.
443
444 Parameters
445 ----------
446 equinox1 = Original equinox of coordinates, numeric scalar.
447 equinox2 = Equinox of precessed coordinates.
448 FK4 = If this keyword is set and non-zero, the FK4 (B1950) system precession angles
449 are used to compute the precession matrix. The default is to use FK5 (J2000) pre-
450 cession angles.
451
452 Return
453 ------
454 r = Precession matrix, used to precess equatorial rectangular coordinates.
455
456 Examples
457 --------
458 >> matrix = premat(1950.0,1975.0,FK4=1)
459 >> print matrix
460 [[ 9.99981438e-01 -5.58774959e-03 -2.42908517e-03]
461 [ 5.58774959e-03 9.99984388e-01 -6.78691471e-06]
462 [ 2.42908517e-03 -6.78633095e-06 9.99997050e-01]]
463
464 Modification history
465 --------------------
466 Written by Wayne Landsman, HSTX Corporation, June 1994.
467 Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
468 """
469
470 t = 0.001*(equinox2 - equinox1)
471
472 if FK4==0:
473 st=0.001*(equinox1 - 2000.)
474 # Computing 3 rotation angles.
475 A=Misc_Routines.CoFactors.s2r*t*(23062.181+st*(139.656+0.0139*st)+t*(30.188-0.344*st+17.998*t))
476 B=Misc_Routines.CoFactors.s2r*t*t*(79.280+0.410*st+0.205*t)+A
477 C=Misc_Routines.CoFactors.s2r*t*(20043.109-st*(85.33+0.217*st)+ t*(-42.665-0.217*st-41.833*t))
478 else:
479 st=0.001*(equinox1 - 1900)
480 # Computing 3 rotation angles
481 A=Misc_Routines.CoFactors.s2r*t*(23042.53+st*(139.75+0.06*st)+t*(30.23-0.27*st+18.0*t))
482 B=Misc_Routines.CoFactors.s2r*t*t*(79.27+0.66*st+0.32*t)+A
483 C=Misc_Routines.CoFactors.s2r*t*(20046.85-st*(85.33+0.37*st)+t*(-42.67-0.37*st-41.8*t))
484
485 sina = numpy.sin(A); sinb = numpy.sin(B); sinc = numpy.sin(C)
486 cosa = numpy.cos(A); cosb = numpy.cos(B); cosc = numpy.cos(C)
487
488 r = numpy.zeros((3,3))
489 r[:,0] = numpy.array([cosa*cosb*cosc-sina*sinb,sina*cosb+cosa*sinb*cosc,cosa*sinc])
490 r[:,1] = numpy.array([-cosa*sinb-sina*cosb*cosc,cosa*cosb-sina*sinb*cosc,-sina*sinc])
491 r[:,2] = numpy.array([-cosb*sinc,-sinb*sinc,cosc])
492
493 return r
494
495 def cirrange(self,angle,rad=0):
496 """
497 cirrange forces an angle into the range 0<= angle < 360.
498
499 Parameters
500 ----------
501 angle = The angle to modify, in degrees. Can be scalar or vector.
502 rad = Set to 1 if the angle is specified in radians rather than degrees. It is for-
503 ced into the range 0 <= angle < 2 PI
504
505 Return
506 ------
507 angle = The angle after the modification.
508
509 Example
510 -------
511 >> angle = cirrange(numpy.array([420,400,361]))
512 >> print angle
513 >> [60, 40, 1]
514
515 Modification History
516 --------------------
517 Written by Michael R. Greason, Hughes STX, 10 February 1994.
518 Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
519 """
520
521 angle = numpy.atleast_1d(angle)
522
523 if rad==1:
524 cnst = numpy.pi*2.
525 elif rad==0:
526 cnst = 360.
527
528 # Deal with the lower limit.
529 angle = angle % cnst
530
531 # Deal with negative values, if way
532 neg = numpy.where(angle<0.0)
533 if neg[0].size>0: angle[neg] = angle[neg] + cnst
534
535 return angle
536
537
538 class CelestialBodies(EquatorialCorrections):
539 def __init__(self):
540 """
541 CelestialBodies class creates a object to call methods of celestial bodies location.
542
543 Modification History
544 --------------------
545 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009.
546 """
547
548 EquatorialCorrections.__init__(self)
549
550 def sunpos(self,jd,rad=0):
551 """
552 sunpos method computes the RA and Dec of the Sun at a given date.
553
554 Parameters
555 ----------
556 jd = The julian date of the day (and time), scalar or vector.
557 rad = If this keyword is set and non-zero, then the input and output RAD and DEC
558 vectors are in radian rather than degree.
559
560 Return
561 ------
562 ra = The right ascension of the sun at that date in degrees.
563 dec = The declination of the sun at that date in degrees.
564 elong = Ecliptic longitude of the sun at that date in degrees.
565 obliquity = The declination of the sun at that date in degrees.
566
567 Examples
568 --------
569 >> jd = 2466880
570 >> [ra,dec,elong,obliquity] = sunpos(jd)
571 >> print ra, dec, elong, obliquity
572 [ 275.53499556] [-23.33840558] [ 275.08917968] [ 23.43596165]
573
574 >> [ra,dec,elong,obliquity] = sunpos(jd,rad=1)
575 >> print ra, dec, elong, obliquity
576 [ 4.80899288] [-0.40733202] [ 4.80121192] [ 0.40903469]
577
578 >> jd = 2450449.5 + numpy.arange(365)
579 >> [ra,dec,elong,obliquity] = sunpos(jd)
580
581 Modification history
582 --------------------
583 Written by Micheal R. Greason, STX Corporation, 28 October 1988.
584 Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
585 """
586
587 jd = numpy.atleast_1d(jd)
588
589 # Form time in Julian centuries from 1900.
590 t = (jd -2415020.0)/36525.0
591
592 # Form sun's mean longitude
593 l = (279.696678+((36000.768925*t) % 360.0))*3600.0
594
595 # Allow for ellipticity of the orbit (equation of centre) using the Earth's mean
596 # anomoly ME
597 me = 358.475844 + ((35999.049750*t) % 360.0)
598 ellcor = (6910.1 - 17.2*t)*numpy.sin(me*Misc_Routines.CoFactors.d2r) + 72.3*numpy.sin(2.0*me*Misc_Routines.CoFactors.d2r)
599 l = l + ellcor
600
601 # Allow for the Venus perturbations using the mean anomaly of Venus MV
602 mv = 212.603219 + ((58517.803875*t) % 360.0)
603 vencorr = 4.8*numpy.cos((299.1017 + mv - me)*Misc_Routines.CoFactors.d2r) + \
604 5.5*numpy.cos((148.3133 + 2.0*mv - 2.0*me )*Misc_Routines.CoFactors.d2r) + \
605 2.5*numpy.cos((315.9433 + 2.0*mv - 3.0*me )*Misc_Routines.CoFactors.d2r) + \
606 1.6*numpy.cos((345.2533 + 3.0*mv - 4.0*me )*Misc_Routines.CoFactors.d2r) + \
607 1.0*numpy.cos((318.15 + 3.0*mv - 5.0*me )*Misc_Routines.CoFactors.d2r)
608 l = l + vencorr
609
610 # Allow for the Mars perturbations using the mean anomaly of Mars MM
611 mm = 319.529425 + ((19139.858500*t) % 360.0)
612 marscorr = 2.0*numpy.cos((343.8883 - 2.0*mm + 2.0*me)*Misc_Routines.CoFactors.d2r ) + \
613 1.8*numpy.cos((200.4017 - 2.0*mm + me)*Misc_Routines.CoFactors.d2r)
614 l = l + marscorr
615
616 # Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ
617 mj = 225.328328 + ((3034.6920239*t) % 360.0)
618 jupcorr = 7.2*numpy.cos((179.5317 - mj + me )*Misc_Routines.CoFactors.d2r) + \
619 2.6*numpy.cos((263.2167 - mj)*Misc_Routines.CoFactors.d2r) + \
620 2.7*numpy.cos((87.1450 - 2.0*mj + 2.0*me)*Misc_Routines.CoFactors.d2r) + \
621 1.6*numpy.cos((109.4933 - 2.0*mj + me)*Misc_Routines.CoFactors.d2r)
622 l = l + jupcorr
623
624 # Allow for Moons perturbations using mean elongation of the Moon from the Sun D
625 d = 350.7376814 + ((445267.11422*t) % 360.0)
626 mooncorr = 6.5*numpy.sin(d*Misc_Routines.CoFactors.d2r)
627 l = l + mooncorr
628
629 # Allow for long period terms
630 longterm = + 6.4*numpy.sin((231.19 + 20.20*t)*Misc_Routines.CoFactors.d2r)
631 l = l + longterm
632 l = (l + 2592000.0) % 1296000.0
633 longmed = l/3600.0
634
635 # Allow for Aberration
636 l = l - 20.5
637
638 # Allow for Nutation using the longitude of the Moons mean node OMEGA
639 omega = 259.183275 - ((1934.142008*t) % 360.0)
640 l = l - 17.2*numpy.sin(omega*Misc_Routines.CoFactors.d2r)
641
642 # Form the True Obliquity
643 oblt = 23.452294 - 0.0130125*t + (9.2*numpy.cos(omega*Misc_Routines.CoFactors.d2r))/3600.0
644
645 # Form Right Ascension and Declination
646 l = l/3600.0
647 ra = numpy.arctan2((numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.cos(oblt*Misc_Routines.CoFactors.d2r)),numpy.cos(l*Misc_Routines.CoFactors.d2r))
648
649 neg = numpy.where(ra < 0.0)
650 if neg[0].size > 0: ra[neg] = ra[neg] + 2.0*numpy.pi
651
652 dec = numpy.arcsin(numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.sin(oblt*Misc_Routines.CoFactors.d2r))
653
654 if rad==1:
655 oblt = oblt*Misc_Routines.CoFactors.d2r
656 longmed = longmed*Misc_Routines.CoFactors.d2r
657 else:
658 ra = ra/Misc_Routines.CoFactors.d2r
659 dec = dec/Misc_Routines.CoFactors.d2r
660
661 return ra, dec, longmed, oblt
662
663 def moonpos(self,jd,rad=0):
664 """
665 moonpos method computes the RA and Dec of the Moon at specified Julian date(s).
666
667 Parameters
668 ----------
669 jd = The julian date of the day (and time), scalar or vector.
670 rad = If this keyword is set and non-zero, then the input and output RAD and DEC
671 vectors are in radian rather than degree.
672
673 Return
674 ------
675 ra = The right ascension of the sun at that date in degrees.
676 dec = The declination of the sun at that date in degrees.
677 dist = The Earth-moon distance in kilometers (between the center of the Earth and
678 the center of the moon).
679 geolon = Apparent longitude of the moon in degrees, referred to the ecliptic of the
680 specified date(s).
681 geolat = Apparent latitude the moon in degrees, referred to the ecliptic of the
682 specified date(s).
683
684 Examples
685 --------
686 >> jd = 2448724.5
687 >> [ra,dec,dist,geolon,geolat] = sunpos(jd)
688 >> print ra, dec, dist, geolon, geolat
689 [ 134.68846855] [ 13.76836663] [ 368409.68481613] [ 133.16726428] [-3.22912642]
690
691 >> [ra,dec,dist,geolon, geolat] = sunpos(jd,rad=1)
692 >> print ra, dec, dist, geolon, geolat
693 [ 2.35075724] [ 0.24030333] [ 368409.68481613] [ 2.32420722] [-0.05635889]
694
695 >> jd = 2450449.5 + numpy.arange(365)
696 >> [ra,dec,dist,geolon, geolat] = sunpos(jd)
697
698 Modification history
699 --------------------
700 Written by Micheal R. Greason, STX Corporation, 31 October 1988.
701 Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
702 """
703
704 jd = numpy.atleast_1d(jd)
705
706 # Form time in Julian centuries from 1900.
707 t = (jd - 2451545.0)/36525.0
708
709 d_lng = numpy.array([0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,\
710 0,0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2])
711
712 m_lng = numpy.array([0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,0,1,-1,0,0,0,1,0,-1,\
713 0,-2,1,2,-2,0,0,-1,0,0,1,-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0])
714
715 mp_lng = numpy.array([1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,-1,0,-1,0,1,2,0,-3,\
716 -2,-1,-2,1,0,2,0,-1,1,0,-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,2,1,-1,3,-1])
717
718 f_lng = numpy.array([0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,\
719 0,0,0,0,-2,2,0,2,0,0,0,0,0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2])
720
721 sin_lng = numpy.array([6288774,1274027,658314,213618,-185116,-114332,58793,57066,\
722 53322,45758,-40923,-34720,-30383,15327,-12528,10980,10675,10034,8548,-7888,\
723 -6766,-5163,4987,4036,3994,3861,3665,-2689,-2602,2390,-2348,2236,-2120,-2069,\
724 2048,-1773,-1595,1215,-1110,-892,-810,759,-713,-700,691,596,549,537,520,-487,\
725 -399,-381,351,-340,330,327,-323,299,294,0.0])
726
727 cos_lng = numpy.array([-20905355,-3699111,-2955968,-569925,48888,-3149,246158,-152138,\
728 -170733,-204586,-129620,108743,104755,10321,0,79661,-34782,-23210,-21636,24208,\
729 30824,-8379,-16675,-12831,-10445,-11650,14403,-7003,0,10056,6322, -9884,5751,0,\
730 -4950,4130,0,-3958,0,3258,2616,-1897,-2117,2354,0,0,-1423,-1117,-1571,-1739,0, \
731 -4421,0,0,0,0,1165,0,0,8752.0])
732
733 d_lat = numpy.array([0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,1,0,0,0,1,0,4,4,0,4,\
734 2,2,2,2,0,2,2,2,2,4,2,2,0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2])
735
736 m_lat = numpy.array([0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,0,1,0,1,1,1,0,0,0,0,\
737 0,0,0,0,-1,0,0,0,0,1,1,0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2])
738
739 mp_lat = numpy.array([0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,0,0,-1,0,1,1,0,0,\
740 3,0,-1,1,-2,0,2,1,-2,3,2,-3,-1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-1,0,0])
741
742 f_lat = numpy.array([1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,1,3,1,1,1,-1,\
743 -1,-1,1,-1,1,-3,1,-3,-1,-1,1,-1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,-1, \
744 -1,-1,-1,-1,1])
745
746 sin_lat = numpy.array([5128122,280602,277693,173237,55413,46271, 32573, 17198, 9266, \
747 8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794, -1749, -1565, -1491, \
748 -1475,-1410,-1344,-1335,1107,1021,833,777,671,607,596,491,-451,439,422,421,-366,\
749 -351,331,315,302,-283,-229,223,223,-220,-220,-185,181,-177,176, 166, -164, 132, \
750 -119,115,107.0])
751
752 # Mean longitude of the moon refered to mean equinox of the date.
753 coeff0 = numpy.array([-1./6.5194e7,1./538841.,-0.0015786,481267.88123421,218.3164477])
754 lprimed = numpy.poly1d(coeff0)
755 lprimed = lprimed(t)
756 lprimed = self.cirrange(lprimed,rad=0)
757 lprime = lprimed*Misc_Routines.CoFactors.d2r
758
759 # Mean elongation of the moon
760 coeff1 = numpy.array([-1./1.13065e8,1./545868.,-0.0018819,445267.1114034,297.8501921])
761 d = numpy.poly1d(coeff1)
762 d = d(t)*Misc_Routines.CoFactors.d2r
763 d = self.cirrange(d,rad=1)
764
765 # Sun's mean anomaly
766 coeff2 = numpy.array([1.0/2.449e7,-0.0001536,35999.0502909,357.5291092])
767 M = numpy.poly1d(coeff2)
768 M = M(t)*Misc_Routines.CoFactors.d2r
769 M = self.cirrange(M,rad=1)
770
771 # Moon's mean anomaly
772 coeff3 = numpy.array([-1.0/1.4712e7,1.0/6.9699e4,0.0087414,477198.8675055,134.9633964])
773 Mprime = numpy.poly1d(coeff3)
774 Mprime = Mprime(t)*Misc_Routines.CoFactors.d2r
775 Mprime = self.cirrange(Mprime,rad=1)
776
777 # Moon's argument of latitude
778 coeff4 = numpy.array([1.0/8.6331e8,-1.0/3.526e7,-0.0036539,483202.0175233,93.2720950])
779 F = numpy.poly1d(coeff4)
780 F = F(t)*Misc_Routines.CoFactors.d2r
781 F = self.cirrange(F,rad=1)
782
783 # Eccentricity of Earth's orbit around the sun
784 e = 1 - 0.002516*t - 7.4e-6*(t**2.)
785 e2 = e**2.
786
787 ecorr1 = numpy.where((numpy.abs(m_lng))==1)
788 ecorr2 = numpy.where((numpy.abs(m_lat))==1)
789 ecorr3 = numpy.where((numpy.abs(m_lng))==2)
790 ecorr4 = numpy.where((numpy.abs(m_lat))==2)
791
792 # Additional arguments.
793 A1 = (119.75 + 131.849*t)*Misc_Routines.CoFactors.d2r
794 A2 = (53.09 + 479264.290*t)*Misc_Routines.CoFactors.d2r
795 A3 = (313.45 + 481266.484*t)*Misc_Routines.CoFactors.d2r
796 suml_add = 3958.*numpy.sin(A1) + 1962.*numpy.sin(lprime - F) + 318*numpy.sin(A2)
797 sumb_add = -2235.*numpy.sin(lprime) + 382.*numpy.sin(A3) + 175.*numpy.sin(A1-F) + \
798 175.*numpy.sin(A1 + F) + 127.*numpy.sin(lprime - Mprime) - 115.*numpy.sin(lprime + Mprime)
799
800 # Sum the periodic terms
801 geolon = numpy.zeros(jd.size)
802 geolat = numpy.zeros(jd.size)
803 dist = numpy.zeros(jd.size)
804
805 for i in numpy.arange(jd.size):
806 sinlng = sin_lng
807 coslng = cos_lng
808 sinlat = sin_lat
809
810 sinlng[ecorr1] = e[i]*sinlng[ecorr1]
811 coslng[ecorr1] = e[i]*coslng[ecorr1]
812 sinlat[ecorr2] = e[i]*sinlat[ecorr2]
813 sinlng[ecorr3] = e2[i]*sinlng[ecorr3]
814 coslng[ecorr3] = e2[i]*coslng[ecorr3]
815 sinlat[ecorr4] = e2[i]*sinlat[ecorr4]
816
817 arg = d_lng*d[i] + m_lng*M[i] + mp_lng*Mprime[i] + f_lng*F[i]
818 geolon[i] = lprimed[i] + (numpy.sum(sinlng*numpy.sin(arg)) + suml_add[i] )/1.e6
819 dist[i] = 385000.56 + numpy.sum(coslng*numpy.cos(arg))/1.e3
820 arg = d_lat*d[i] + m_lat*M[i] + mp_lat*Mprime[i] + f_lat*F[i]
821 geolat[i] = (numpy.sum(sinlat*numpy.sin(arg)) + sumb_add[i])/1.e6
822
823 [nlon, elon] = self.nutate(jd)
824 geolon = geolon + nlon/3.6e3
825 geolon = self.cirrange(geolon,rad=0)
826 lamb = geolon*Misc_Routines.CoFactors.d2r
827 beta = geolat*Misc_Routines.CoFactors.d2r
828
829 # Find mean obliquity and convert lamb, beta to RA, Dec
830 c = numpy.array([2.45,5.79,27.87,7.12,-39.05,-249.67,-51.38,1999.25,-1.55,-4680.93, \
831 21.448])
832 junk = numpy.poly1d(c);
833 epsilon = 23. + (26./60.) + (junk(t/1.e2)/3600.)
834 # True obliquity in radians
835 eps = (epsilon + elon/3600. )*Misc_Routines.CoFactors.d2r
836
837 ra = numpy.arctan2(numpy.sin(lamb)*numpy.cos(eps)-numpy.tan(beta)*numpy.sin(eps),numpy.cos(lamb))
838 ra = self.cirrange(ra,rad=1)
839
840 dec = numpy.arcsin(numpy.sin(beta)*numpy.cos(eps) + numpy.cos(beta)*numpy.sin(eps)*numpy.sin(lamb))
841
842 if rad==1:
843 geolon = lamb
844 geolat = beta
845 else:
846 ra = ra/Misc_Routines.CoFactors.d2r
847 dec = dec/Misc_Routines.CoFactors.d2r
848
849 return ra, dec, dist, geolon, geolat
850
851 def hydrapos(self):
852 """
853 hydrapos method returns RA and Dec provided by Bill Coles (Oct 2003).
854
855 Parameters
856 ----------
857 None
858
859 Return
860 ------
861 ra = The right ascension of the sun at that date in degrees.
862 dec = The declination of the sun at that date in degrees.
863 Examples
864 --------
865 >> [ra,dec] = hydrapos()
866 >> print ra, dec
867 139.45 -12.0833333333
868
869 Modification history
870 --------------------
871 Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
872 """
873
874 ra = (9. + 17.8/60.)*15.
875 dec = -(12. + 5./60.)
876
877 return ra, dec
878
879
880 def skynoise_jro(self,dec_cut=-11.95,filename='skynoise_jro.dat',filepath=None):
881 """
882 hydrapos returns RA and Dec provided by Bill Coles (Oct 2003).
883
884 Parameters
885 ----------
886 dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica-
887 marca. The default value is -11.95.
888 filename = A string to specify name the skynoise file. The default value is skynoi-
889 se_jro.dat
890
891 Return
892 ------
893 maxra = The maximum right ascension to the declination used to get a cut.
894 ra = The right ascension.
895 Examples
896 --------
897 >> [maxra,ra] = skynoise_jro()
898 >> print maxra, ra
899 139.45 -12.0833333333
900
901 Modification history
902 --------------------
903 Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
904 """
905
906 if filepath==None:
907 filepath = '/app/utils/'
908
909 f = open(os.path.join(filepath,filename),'rb')
910
911 # Reading SkyNoise Power (lineal scale)
912 ha_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20)
913 ha_sky = ha_sky['var'].reshape(20,480).transpose()
914
915 dec_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20)
916 dec_sky = dec_sky['var'].reshape((20,480)).transpose()
917
918 tmp_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20)
919 tmp_sky = tmp_sky['var'].reshape((20,480)).transpose()
920
921 f.close()
922
923 nha = 480
924 tmp_cut = numpy.zeros(nha)
925 for iha in numpy.arange(nha):
926 tck = scipy.interpolate.splrep(dec_sky[iha,:],tmp_sky[iha,:],s=0)
927 tmp_cut[iha] = scipy.interpolate.splev(dec_cut,tck,der=0)
928
929 ptr = numpy.nanargmax(tmp_cut)
930
931 maxra = ha_sky[ptr,0]
932 ra = ha_sky[:,0]
933
934 return maxra, ra
935
936 def skyNoise(self,jd,ut=-5.0,longitude=-76.87,filename='galaxy.txt',filepath=None):
937 """
938 hydrapos returns RA and Dec provided by Bill Coles (Oct 2003).
939
940 Parameters
941 ----------
942 jd = The julian date of the day (and time), scalar or vector.
943
944 dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica-
945 marca. The default value is -11.95.
946 filename = A string to specify name the skynoise file. The default value is skynoi-
947 se_jro.dat
948
949 Return
950 ------
951 maxra = The maximum right ascension to the declination used to get a cut.
952 ra = The right ascension.
953
954 Examples
955 --------
956 >> [maxra,ra] = skynoise_jro()
957 >> print maxra, ra
958 139.45 -12.0833333333
959
960 Modification history
961 --------------------
962 Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
963 """
964
965 # Defining date to compute SkyNoise.
966 [year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time()
967 is_dom = (month==9) & (dom==21)
968 if is_dom:
969 tmp = jd
970 jd = TimeTools.Time(year,9,22).change2julian()
971 dom = 22
972
973 # Reading SkyNoise
974 if filepath==None:filepath='./resource'
975 f = open(os.path.join(filepath,filename))
976
977 lines = f.read()
978 f.close()
979
980 nlines = 99
981 lines = lines.split('\n')
982 data = numpy.zeros((2,nlines))*numpy.float32(0.)
983 for ii in numpy.arange(nlines):
984 line = numpy.array([lines[ii][0:6],lines[ii][6:]])
985 data[:,ii] = numpy.float32(line)
986
987 # Getting SkyNoise to the date desired.
988 otime = data[0,:]*60.0
989 opowr = data[1,:]
990
991 hour = numpy.array([0,23]);
992 mins = numpy.array([0,59]);
993 secs = numpy.array([0,59]);
994 LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday()
995 LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1))
996 lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst()
997
998 ipowr = lst*0.0
999 # Interpolating using scipy (inside max and min "x")
1000 otime = otime/3600.
1001 val = numpy.where((lst>numpy.min(otime)) & (lst<numpy.max(otime))); val = val[0]
1002 tck = scipy.interpolate.interp1d(otime,opowr)
1003 ipowr[val] = tck(lst[val])
1004
1005 # Extrapolating above maximum time data (23.75).
1006 uval = numpy.where(lst>numpy.max(otime))
1007 if uval[0].size>0:
1008 ii = numpy.min(uval[0])
1009 m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2])
1010 b = ipowr[ii-1] - m*lst[ii-1]
1011 ipowr[uval] = m*lst[uval] + b
1012
1013 if is_dom:
1014 lst = numpy.roll(lst,4)
1015 ipowr = numpy.roll(ipowr,4)
1016
1017 new_lst = numpy.int32(lst*3600.)
1018 new_pow = ipowr
1019
1020 return ipowr, LTtime, lst
1021
1022
1023 class AltAz(EquatorialCorrections):
1024 def __init__(self,alt,az,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\
1025 aberration_=0,B1950=0):
1026 """
1027 The AltAz class creates an object which represents the target position in horizontal
1028 coordinates (alt-az) and allows to convert (using the methods) from this coordinate
1029 system to others (e.g. Equatorial).
1030
1031 Parameters
1032 ----------
1033 alt = Altitude in degrees. Scalar or vector.
1034 az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca-
1035 lar or vector.
1036 jd = Julian date. Scalar or vector.
1037 lat = North geodetic latitude of location in degrees. The default value is -11.95.
1038 lon = East longitude of location in degrees. The default value is -76.8667.
1039 WS = Set this to 1 to get the azimuth measured westward from south.
1040 altitude = The altitude of the observing location, in meters. The default 500.
1041 nutate_ = Set this to 1 to force nutation, 0 for no nutation.
1042 precess_ = Set this to 1 to force precession, 0 for no precession.
1043 aberration_ = Set this to 1 to force aberration correction, 0 for no correction.
1044 B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins-
1045 tead of J2000, FK5)
1046
1047 Modification History
1048 --------------------
1049 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 26 September 2009.
1050 """
1051
1052 EquatorialCorrections.__init__(self)
1053
1054 self.alt = numpy.atleast_1d(alt)
1055 self.az = numpy.atleast_1d(az)
1056 self.jd = numpy.atleast_1d(jd)
1057 self.lat = lat
1058 self.lon = lon
1059 self.WS = WS
1060 self.altitude = altitude
1061
1062 self.nutate_ = nutate_
1063 self.aberration_ = aberration_
1064 self.precess_ = precess_
1065 self.B1950 = B1950
1066
1067 def change2equatorial(self):
1068 """
1069 change2equatorial method converts horizon (Alt-Az) coordinates to equatorial coordi-
1070 nates (ra-dec).
1071
1072 Return
1073 ------
1074 ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector.
1075 dec = Declination of object (J2000), in degrees (FK5). Scalar or vector.
1076 ha = Hour angle in degrees.
1077
1078 Example
1079 -------
1080 >> alt = 88.5401
1081 >> az = -128.990
1082 >> jd = 2452640.5
1083 >> ObjAltAz = AltAz(alt,az,jd)
1084 >> [ra, dec, ha] = ObjAltAz.change2equatorial()
1085 >> print ra, dec, ha
1086 [ 22.20280632] [-12.86610025] [ 1.1638927]
1087
1088 Modification History
1089 --------------------
1090 Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002.
1091 Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
1092 """
1093
1094 az = self.az
1095 alt = self.alt
1096 if self.WS>0:az = az -180.
1097 ra_tmp = numpy.zeros(numpy.size(self.jd)) + 45.
1098 dec_tmp = numpy.zeros(numpy.size(self.jd)) + 45.
1099 [dra1,ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra_tmp, dec_tmp)
1100
1101 # Getting local mean sidereal time (lmst)
1102 lmst = TimeTools.Julian(self.jd[0]).change2lst()
1103 lmst = lmst*Misc_Routines.CoFactors.h2d
1104 # Getting local apparent sidereal time (last)
1105 last = lmst + d_psi*numpy.cos(eps)/3600.
1106
1107 # Now do the spherical trig to get APPARENT hour angle and declination (Degrees).
1108 [ha, dec] = self.change2HaDec()
1109
1110 # Finding Right Ascension (in degrees, from 0 to 360.)
1111 ra = (last - ha + 360.) % 360.
1112
1113 # Calculate NUTATION and ABERRATION Correction to Ra-Dec
1114 [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec)
1115 [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec)
1116
1117 # Make Nutation and Aberration correction (if wanted)
1118 ra = ra - (dra1*self.nutate_ + dra2*self.aberration_)/3600.
1119 dec = dec - (ddec1*self.nutate_ + ddec2*self.aberration_)/3600.
1120
1121 # Computing current equinox
1122 j_now = (self.jd - 2451545.)/365.25 + 2000
1123
1124 # Precess coordinates to current date
1125 if self.precess_==1:
1126 njd = numpy.size(self.jd)
1127 for ii in numpy.arange(njd):
1128 ra_i = ra[ii]
1129 dec_i = dec[ii]
1130 now = j_now[ii]
1131
1132 if self.B1950==1:
1133 [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1)
1134 elif self.B1950==0:
1135 [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0)
1136
1137 ra[ii] = ra_i
1138 dec[ii] = dec_i
1139
1140 return ra, dec, ha
1141
1142 def change2HaDec(self):
1143 """
1144 change2HaDec method converts from horizon (Alt-Az) coordinates to hour angle and de-
1145 clination.
1146
1147 Return
1148 ------
1149 ha = The local apparent hour angle, in degrees. The hour angle is the time that ri-
1150 ght ascension of 0 hours crosses the local meridian. It is unambiguisoly defined.
1151 dec = The local apparent declination, in degrees.
1152
1153 Example
1154 -------
1155 >> alt = 88.5401
1156 >> az = -128.990
1157 >> jd = 2452640.5
1158 >> ObjAltAz = AltAz(alt,az,jd)
1159 >> [ha, dec] = ObjAltAz.change2HaDec()
1160 >> print ha, dec
1161 [ 1.1638927] [-12.86610025]
1162
1163 Modification History
1164 --------------------
1165 Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002.
1166 Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
1167 """
1168
1169 alt_r = numpy.atleast_1d(self.alt*Misc_Routines.CoFactors.d2r)
1170 az_r = numpy.atleast_1d(self.az*Misc_Routines.CoFactors.d2r)
1171 lat_r = numpy.atleast_1d(self.lat*Misc_Routines.CoFactors.d2r)
1172
1173 # Find local hour angle (in degrees, from 0 to 360.)
1174 y_ha = -1*numpy.sin(az_r)*numpy.cos(alt_r)
1175 x_ha = -1*numpy.cos(az_r)*numpy.sin(lat_r)*numpy.cos(alt_r) + numpy.sin(alt_r)*numpy.cos(lat_r)
1176
1177 ha = numpy.arctan2(y_ha,x_ha)
1178 ha = ha/Misc_Routines.CoFactors.d2r
1179
1180 w = numpy.where(ha<0.)
1181 if w[0].size>0:ha[w] = ha[w] + 360.
1182 ha = ha % 360.
1183
1184 # Find declination (positive if north of celestial equatorial, negative if south)
1185 sindec = numpy.sin(lat_r)*numpy.sin(alt_r) + numpy.cos(lat_r)*numpy.cos(alt_r)*numpy.cos(az_r)
1186 dec = numpy.arcsin(sindec)/Misc_Routines.CoFactors.d2r
1187
1188 return ha, dec
1189
1190
1191 class Equatorial(EquatorialCorrections):
1192 def __init__(self,ra,dec,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\
1193 aberration_=0,B1950=0):
1194 """
1195 The Equatorial class creates an object which represents the target position in equa-
1196 torial coordinates (ha-dec) and allows to convert (using the class methods) from
1197 this coordinate system to others (e.g. AltAz).
1198
1199 Parameters
1200 ----------
1201 ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector.
1202 dec = Declination of object (J2000), in degrees (FK5). Scalar or vector.
1203 jd = Julian date. Scalar or vector.
1204 lat = North geodetic latitude of location in degrees. The default value is -11.95.
1205 lon = East longitude of location in degrees. The default value is -76.8667.
1206 WS = Set this to 1 to get the azimuth measured westward from south.
1207 altitude = The altitude of the observing location, in meters. The default 500.
1208 nutate = Set this to 1 to force nutation, 0 for no nutation.
1209 precess = Set this to 1 to force precession, 0 for no precession.
1210 aberration = Set this to 1 to force aberration correction, 0 for no correction.
1211 B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins-
1212 tead of J2000, FK5)
1213
1214 Modification History
1215 --------------------
1216 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 29 September 2009.
1217 """
1218
1219 EquatorialCorrections.__init__(self)
1220
1221 self.ra = numpy.atleast_1d(ra)
1222 self.dec = numpy.atleast_1d(dec)
1223 self.jd = numpy.atleast_1d(jd)
1224 self.lat = lat
1225 self.lon = lon
1226 self.WS = WS
1227 self.altitude = altitude
1228
1229 self.nutate_ = nutate_
1230 self.aberration_ = aberration_
1231 self.precess_ = precess_
1232 self.B1950 = B1950
1233
1234 def change2AltAz(self):
1235 """
1236 change2AltAz method converts from equatorial coordinates (ha-dec) to horizon coordi-
1237 nates (alt-az).
1238
1239 Return
1240 ------
1241 alt = Altitude in degrees. Scalar or vector.
1242 az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca-
1243 lar or vector.
1244 ha = Hour angle in degrees.
1245
1246 Example
1247 -------
1248 >> ra = 43.370609
1249 >> dec = -28.0000
1250 >> jd = 2452640.5
1251 >> ObjEq = Equatorial(ra,dec,jd)
1252 >> [alt, az, ha] = ObjEq.change2AltAz()
1253 >> print alt, az, ha
1254 [ 65.3546497] [ 133.58753124] [ 339.99609002]
1255
1256 Modification History
1257 --------------------
1258 Written Chris O'Dell Univ. of Wisconsin-Madison. May 2002
1259 Converted to Python by Freddy R. Galindo, ROJ, 29 September 2009.
1260 """
1261
1262 ra = self.ra
1263 dec = self.dec
1264
1265 # Computing current equinox
1266 j_now = (self.jd - 2451545.)/365.25 + 2000
1267
1268 # Precess coordinates to current date
1269 if self.precess_==1:
1270 njd = numpy.size(self.jd)
1271 for ii in numpy.arange(njd):
1272 ra_i = ra[ii]
1273 dec_i = dec[ii]
1274 now = j_now[ii]
1275
1276 if self.B1950==1:
1277 [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1)
1278 elif self.B1950==0:
1279 [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0)
1280
1281 ra[ii] = ra_i
1282 dec[ii] = dec_i
1283
1284 # Calculate NUTATION and ABERRATION Correction to Ra-Dec
1285 [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec)
1286 [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec)
1287
1288 # Make Nutation and Aberration correction (if wanted)
1289 ra = ra + (dra1*self.nutate_ + dra2*self.aberration_)/3600.
1290 dec = dec + (ddec1*self.nutate_ + ddec2*self.aberration_)/3600.
1291
1292 # Getting local mean sidereal time (lmst)
1293 lmst = TimeTools.Julian(self.jd).change2lst()
1294
1295 lmst = lmst*Misc_Routines.CoFactors.h2d
1296 # Getting local apparent sidereal time (last)
1297 last = lmst + d_psi*numpy.cos(eps)/3600.
1298
1299 # Finding Hour Angle (in degrees, from 0 to 360.)
1300 ha = last - ra
1301 w = numpy.where(ha<0.)
1302 if w[0].size>0:ha[w] = ha[w] + 360.
1303 ha = ha % 360.
1304
1305 # Now do the spherical trig to get APPARENT hour angle and declination (Degrees).
1306 [alt, az] = self.HaDec2AltAz(ha,dec)
1307
1308 return alt, az, ha
1309
1310 def HaDec2AltAz(self,ha,dec):
1311 """
1312 HaDec2AltAz convert hour angle and declination (ha-dec) to horizon coords (alt-az).
1313
1314 Parameters
1315 ----------
1316 ha = The local apparent hour angle, in DEGREES, scalar or vector.
1317 dec = The local apparent declination, in DEGREES, scalar or vector.
1318
1319 Return
1320 ------
1321 alt = Altitude in degrees. Scalar or vector.
1322 az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca-
1323 lar or vector.
1324
1325 Modification History
1326 --------------------
1327 Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002.
1328 Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
1329 """
1330
1331 sh = numpy.sin(ha*Misc_Routines.CoFactors.d2r) ; ch = numpy.cos(ha*Misc_Routines.CoFactors.d2r)
1332 sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r) ; cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r)
1333 sl = numpy.sin(self.lat*Misc_Routines.CoFactors.d2r) ; cl = numpy.cos(self.lat*Misc_Routines.CoFactors.d2r)
1334
1335 x = -1*ch*cd*sl + sd*cl
1336 y = -1*sh*cd
1337 z = ch*cd*cl + sd*sl
1338 r = numpy.sqrt(x**2. + y**2.)
1339
1340 az = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r
1341 alt = numpy.arctan2(z,r)/Misc_Routines.CoFactors.d2r
1342
1343 # correct for negative az.
1344 w = numpy.where(az<0.)
1345 if w[0].size>0:az[w] = az[w] + 360.
1346
1347 # Convert az to West from South, if desired
1348 if self.WS==1: az = (az + 180.) % 360.
1349
1350 return alt, az
1351
1352
1353 class Geodetic():
1354 def __init__(self,lat=-11.95,alt=0):
1355 """
1356 The Geodetic class creates an object which represents the real position on earth of
1357 a target (Geodetic Coordinates: lat-alt) and allows to convert (using the class me-
1358 thods) from this coordinate system to others (e.g. geocentric).
1359
1360 Parameters
1361 ----------
1362 lat = Geodetic latitude of location in degrees. The default value is -11.95.
1363
1364 alt = Geodetic altitude (km). The default value is 0.
1365
1366 Modification History
1367 --------------------
1368 Converted to Object-oriented Programming by Freddy R. Galindo, ROJ, 02 October 2009.
1369 """
1370
1371 self.lat = numpy.atleast_1d(lat)
1372 self.alt = numpy.atleast_1d(alt)
1373
1374 self.a = 6378.16
1375 self.ab2 = 1.0067397
1376 self.ep2 = 0.0067397
1377
1378 def change2geocentric(self):
1379 """
1380 change2geocentric method converts from Geodetic to Geocentric coordinates. The re-
1381 ference geoid is that adopted by the IAU in 1964.
1382
1383 Return
1384 ------
1385 gclat = Geocentric latitude (in degrees), scalar or vector.
1386 gcalt = Geocentric radial distance (km), scalar or vector.
1387
1388 Example
1389 -------
1390 >> ObjGeoid = Geodetic(lat=-11.95,alt=0)
1391 >> [gclat, gcalt] = ObjGeoid.change2geocentric()
1392 >> print gclat, gcalt
1393 [-11.87227742] [ 6377.25048195]
1394
1395 Modification History
1396 --------------------
1397 Converted to Python by Freddy R. Galindo, ROJ, 02 October 2009.
1398 """
1399
1400 gdl = self.lat*Misc_Routines.CoFactors.d2r
1401 slat = numpy.sin(gdl)
1402 clat = numpy.cos(gdl)
1403 slat2 = slat**2.
1404 clat2 = (self.ab2*clat)**2.
1405
1406 sbet = slat/numpy.sqrt(slat2 + clat2)
1407 sbet2 = (sbet**2.) # < 1
1408 noval = numpy.where(sbet2>1)
1409 if noval[0].size>0:sbet2[noval] = 1
1410 cbet = numpy.sqrt(1. - sbet2)
1411
1412 rgeoid = self.a/numpy.sqrt(1. + self.ep2*sbet2)
1413
1414 x = rgeoid*cbet + self.alt*clat
1415 y = rgeoid*sbet + self.alt*slat
1416
1417 gcalt = numpy.sqrt(x**2. + y**2.)
1418 gclat = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r
1419
1420 return gclat, gcalt
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (507 lines changed) Show them Hide them
@@ -0,0 +1,507
1 """
2
3 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009.
4 Added to signal Chain by Joab Apaza, ROJ, Jun 2023.
5 """
6
7 import numpy
8 import time
9 import os
10 from scipy.special import lpmn
11 from schainpy.model.utils import Astro_Coords
12
13 class BField():
14 def __init__(self,year=None,doy=None,site=1,heights=None,alpha_i=90):
15 """
16 BField class creates an object to get the Magnetic field for a specific date and
17 height(s).
18
19 Parameters
20 ----------
21 year = A scalar giving the desired year. If the value is None (default value) then
22 the current year will be used.
23 doy = A scalar giving the desired day of the year. If the value is None (default va-
24 lue) then the current doy will be used.
25 site = An integer to choose the geographic coordinates of the place where the magne-
26 tic field will be computed. The default value is over Jicamarca (site=1)
27 heights = An array giving the heights (km) where the magnetic field will be modeled By default the magnetic field will be computed at 100, 500 and 1000km.
28 alpha_i = Angle to interpolate the magnetic field.
29
30 Modification History
31 --------------------
32 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009.
33 Added to signal Chain by Joab Apaza, ROJ, Jun 2023.
34 """
35
36 tmp = time.localtime()
37 if year==None: year = tmp[0]
38 if doy==None: doy = tmp[7]
39 self.year = year
40 self.doy = doy
41 self.site = site
42 if heights is None:
43 heights = numpy.array([100,500,1000])
44 else:
45 heights = numpy.array(heights)
46 self.heights = heights
47 self.alpha_i = alpha_i
48
49 def getBField(self,maglimits=numpy.array([-37,-37,37,37])):
50 """
51 getBField models the magnetic field for a different heights in a specific date.
52
53 Parameters
54 ----------
55 maglimits = An 4-elements array giving ..... The default value is [-7,-7,7,7].
56
57 Return
58 ------
59 dcos = An 4-dimensional array giving the directional cosines of the magnetic field
60 over the desired place.
61 alpha = An 3-dimensional array giving the angle of the magnetic field over the desi-
62 red place.
63
64 Modification History
65 --------------------
66 Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009.
67 """
68
69 x_ant = numpy.array([1,0,0])
70 y_ant = numpy.array([0,1,0])
71 z_ant = numpy.array([0,0,1])
72
73 if self.site==0:
74 title_site = "Magnetic equator"
75 coord_site = numpy.array([-76+52./60.,-11+57/60.,0.5])
76 elif self.site==1:
77 title_site = 'Jicamarca'
78 coord_site = [-76-52./60.,-11-57/60.,0.5]
79 heta = (45+5.35)*numpy.pi/180. # (50.35 and 1.46 from Fleish Thesis)
80 delta = -1.46*numpy.pi/180
81
82
83 x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1)
84 y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1)
85 z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1)
86
87 ang0 = -1*coord_site[0]*numpy.pi/180.
88 ang1 = coord_site[1]*numpy.pi/180.
89 x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0)
90 y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0)
91 z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0)
92
93 elif self.site==2: #AMISR
94 title_site = 'AMISR 14'
95
96 coord_site = [-76.874913, -11.953371, 0.52984]
97
98 theta = (0.0977)*numpy.pi/180. # 0.0977
99 delta = 0.110*numpy.pi/180 # 0.11
100
101 x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1)
102 y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1)
103 z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1)
104
105 ang0 = -1*coord_site[0]*numpy.pi/180.
106 ang1 = coord_site[1]*numpy.pi/180.
107 x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0)
108 y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0)
109 z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0)
110 else:
111 # print "No defined Site. Skip..."
112 return None
113
114 nhei = self.heights.size
115 pt_intercep = numpy.zeros((nhei,2))
116 nfields = 1
117
118 grid_res = 2.5
119 nlon = int(int(maglimits[2] - maglimits[0])/grid_res + 1)
120 nlat = int(int(maglimits[3] - maglimits[1])/grid_res + 1)
121
122 location = numpy.zeros((nlon,nlat,2))
123 mlon = numpy.atleast_2d(numpy.arange(nlon)*grid_res + maglimits[0])
124 mrep = numpy.atleast_2d(numpy.zeros(nlat) + 1)
125 location0 = numpy.dot(mlon.transpose(),mrep)
126
127 mlat = numpy.atleast_2d(numpy.arange(nlat)*grid_res + maglimits[1])
128 mrep = numpy.atleast_2d(numpy.zeros(nlon) + 1)
129 location1 = numpy.dot(mrep.transpose(),mlat)
130
131 location[:,:,0] = location0
132 location[:,:,1] = location1
133
134 alpha = numpy.zeros((nlon,nlat,nhei))
135 rr = numpy.zeros((nlon,nlat,nhei,3))
136 dcos = numpy.zeros((nlon,nlat,nhei,2))
137
138 global first_time
139
140 first_time = None
141 for ilon in numpy.arange(nlon):
142 for ilat in numpy.arange(nlat):
143 outs = self.__bdotk(self.heights,
144 self.year + self.doy/366.,
145 coord_site[1],
146 coord_site[0],
147 coord_site[2],
148 coord_site[1]+location[ilon,ilat,1],
149 location[ilon,ilat,0]*720./180.)
150
151 alpha[ilon, ilat,:] = outs[1]
152 rr[ilon, ilat,:,:] = outs[3]
153
154 mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose()
155 tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(x_ant))
156 tmp = tmp.sum(axis=1)
157 dcos[ilon,ilat,:,0] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1))
158
159 mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose()
160 tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(y_ant))
161 tmp = tmp.sum(axis=1)
162 dcos[ilon,ilat,:,1] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1))
163
164 return dcos, alpha, nlon, nlat
165
166
167 def __bdotk(self,heights,tm,gdlat=-11.95,gdlon=-76.8667,gdalt=0.0,decd=-12.88, ham=-4.61666667):
168
169 global first_time
170 # Mean Earth radius in Km WGS 84
171 a_igrf = 6371.2
172
173 bk = numpy.zeros(heights.size)
174 alpha = numpy.zeros(heights.size)
175 bfm = numpy.zeros(heights.size)
176 rr = numpy.zeros((heights.size,3))
177 rgc = numpy.zeros((heights.size,3))
178
179 ObjGeodetic = Astro_Coords.Geodetic(gdlat,gdalt)
180 [gclat,gcalt] = ObjGeodetic.change2geocentric()
181
182 gclat = gclat*numpy.pi/180.
183 gclon = gdlon*numpy.pi/180.
184
185 # Antenna position from center of Earth
186 ca_vector = [numpy.cos(gclat)*numpy.cos(gclon),numpy.cos(gclat)*numpy.sin(gclon),numpy.sin(gclat)]
187 ca_vector = gcalt*numpy.array(ca_vector)
188
189 dec = decd*numpy.pi/180.
190
191 # K vector respect to the center of earth.
192 klon = gclon + ham*numpy.pi/720.
193 k_vector = [numpy.cos(dec)*numpy.cos(klon),numpy.cos(dec)*numpy.sin(klon),numpy.sin(dec)]
194 k_vector = numpy.array(k_vector)
195
196 for ih in numpy.arange(heights.size):
197 # Vector from Earth's center to volume of interest
198 rr[ih,:] = k_vector*heights[ih]
199 cv_vector = numpy.squeeze(ca_vector) + rr[ih,:]
200
201 cv_gcalt = numpy.sqrt(numpy.sum(cv_vector**2.))
202 cvxy = numpy.sqrt(numpy.sum(cv_vector[0:2]**2.))
203
204 radial = cv_vector/cv_gcalt
205 east = numpy.array([-1*cv_vector[1],cv_vector[0],0])/cvxy
206 comp1 = east[1]*radial[2] - radial[1]*east[2]
207 comp2 = east[2]*radial[0] - radial[2]*east[0]
208 comp3 = east[0]*radial[1] - radial[0]*east[1]
209 north = -1*numpy.array([comp1, comp2, comp3])
210
211 rr_k = cv_vector - numpy.squeeze(ca_vector)
212 u_rr = rr_k/numpy.sqrt(numpy.sum(rr_k**2.))
213
214 cv_gclat = numpy.arctan2(cv_vector[2],cvxy)
215 cv_gclon = numpy.arctan2(cv_vector[1],cv_vector[0])
216
217 bhei = cv_gcalt-a_igrf
218 blat = cv_gclat*180./numpy.pi
219 blon = cv_gclon*180./numpy.pi
220 bfield = self.__igrfkudeki(bhei,tm,blat,blon)
221
222 B = (bfield[0]*north + bfield[1]*east - bfield[2]*radial)*1.0e-5
223
224 bfm[ih] = numpy.sqrt(numpy.sum(B**2.)) #module
225 bk[ih] = numpy.sum(u_rr*B)
226 alpha[ih] = numpy.arccos(bk[ih]/bfm[ih])*180/numpy.pi
227 rgc[ih,:] = numpy.array([cv_gclon, cv_gclat, cv_gcalt])
228
229 return bk, alpha, bfm, rr, rgc
230
231
232 def __igrfkudeki(self,heights,time,latitude,longitude,ae=6371.2):
233 """
234 __igrfkudeki calculates the International Geomagnetic Reference Field for given in-
235 put conditions based on IGRF2005 coefficients.
236
237 Parameters
238 ----------
239 heights = Scalar or vector giving the height above the Earth of the point in ques-
240 tion in kilometers.
241 time = Scalar or vector giving the decimal year of time in question (e.g. 1991.2).
242 latitude = Latitude of point in question in decimal degrees. Scalar or vector.
243 longitude = Longitude of point in question in decimal degrees. Scalar or vector.
244 ae =
245 first_time =
246
247 Return
248 ------
249 bn =
250 be =
251 bd =
252 bmod =
253 balpha =
254 first_time =
255
256 Modification History
257 --------------------
258 Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009.
259 """
260
261 global first_time
262 global gs, hs, nvec, mvec, maxcoef
263
264 heights = numpy.atleast_1d(heights)
265 time = numpy.atleast_1d(time)
266 latitude = numpy.atleast_1d(latitude)
267 longitude = numpy.atleast_1d(longitude)
268
269 if numpy.max(latitude)==90:
270 # print "Field calculations are not supported at geographic poles"
271 pass
272
273 # output arrays
274 bn = numpy.zeros(heights.size)
275 be = numpy.zeros(heights.size)
276 bd = numpy.zeros(heights.size)
277
278 if first_time==None:first_time=0
279
280 time0 = time[0]
281 if time!=first_time:
282 #print "Getting coefficients for", time0
283 [periods,g,h ] = self.__readIGRFcoeff()
284 top_year = numpy.max(periods)
285 nperiod = (top_year - 1900)/5 + 1
286
287 maxcoef = 10
288
289 if time0>=2000:maxcoef = 12
290
291
292 # Normalization array for Schmidt fucntions
293 multer = numpy.zeros((2+maxcoef,1+maxcoef)) + 1
294 for cn in (numpy.arange(maxcoef)+1):
295 for rm in (numpy.arange(cn)+1):
296 tmp = numpy.arange(2*rm) + cn - rm + 1.
297 multer[rm+1,cn] = ((-1.)**rm)*numpy.sqrt(2./tmp.prod())
298
299 schmidt = multer[1:,1:].transpose()
300
301 # n and m arrays
302 nvec = numpy.atleast_2d(numpy.arange(maxcoef)+2)
303 mvec = numpy.atleast_2d(numpy.arange(maxcoef+1)).transpose()
304
305 # Time adjusted igrf g and h with Schmidt normalization
306 # IGRF coefficient arrays: g0(n,m), n=1, maxcoeff,m=0, maxcoeff, ...
307 if time0<top_year:
308 dtime = (time0 - 1900) % 5
309 ntime = int((time0 - 1900 - dtime)/5)
310 else:
311 # Estimating coefficients for times > top_year
312 dtime = (time0 - top_year) + 5
313 ntime = int(g[:,0,0].size - 2)
314
315
316 g0 = g[ntime,1:maxcoef+1,:maxcoef+1]
317 h0 = h[ntime,1:maxcoef+1,:maxcoef+1]
318 gdot = g[ntime+1,1:maxcoef+1,:maxcoef+1]-g[ntime,1:maxcoef+1,:maxcoef+1]
319 hdot = h[ntime+1,1:maxcoef+1,:maxcoef+1]-h[ntime,1:maxcoef+1,:maxcoef+1]
320 gs = (g0 + dtime*(gdot/5.))*schmidt[:maxcoef,0:maxcoef+1]
321 hs = (h0 + dtime*(hdot/5.))*schmidt[:maxcoef,0:maxcoef+1]
322
323 first_time = time0
324
325 for ii in numpy.arange(heights.size):
326 # Height dependence array rad = (ae/(ae+height))**(n+3)
327 rad = numpy.atleast_2d((ae/(ae + heights[ii]))**(nvec+1))
328
329 # Sin and Cos of m times longitude phi arrays
330 mphi = mvec*longitude[ii]*numpy.pi/180.
331 cosmphi = numpy.atleast_2d(numpy.cos(mphi))
332 sinmphi = numpy.atleast_2d(numpy.sin(mphi))
333
334 # Cos of colatitude theta
335 c = numpy.cos((90 - latitude[ii])*numpy.pi/180.)
336
337 # Legendre functions p(n,m|c)
338 [p,dp]= lpmn(maxcoef+1,maxcoef+1,c)
339 p = p[:,:-1].transpose()
340 s = numpy.sqrt((1. - c)*(1 + c))
341
342 # Generate derivative array dpdtheta = -s*dpdc
343 dpdtheta = c*p/s
344 for m in numpy.arange(maxcoef+2): dpdtheta[:,m] = m*dpdtheta[:,m]
345 dpdtheta = dpdtheta + numpy.roll(p,-1,axis=1)
346
347 # Extracting arrays required for field calculations
348 p = p[1:maxcoef+1,:maxcoef+1]
349 dpdtheta = dpdtheta[1:maxcoef+1,:maxcoef+1]
350
351 # Weigh p and dpdtheta with gs and hs coefficients.
352 gp = gs*p
353 hp = hs*p
354 gdpdtheta = gs*dpdtheta
355 hdpdtheta = hs*dpdtheta
356 # Calcultate field components
357 matrix0 = numpy.dot(gdpdtheta,cosmphi)
358 matrix1 = numpy.dot(hdpdtheta,sinmphi)
359 bn[ii] = numpy.dot(rad,(matrix0 + matrix1))
360 matrix0 = numpy.dot(hp,(mvec*cosmphi))
361 matrix1 = numpy.dot(gp,(mvec*sinmphi))
362 be[ii] = numpy.dot((-1*rad),((matrix0 - matrix1)/s))
363 matrix0 = numpy.dot(gp,cosmphi)
364 matrix1 = numpy.dot(hp,sinmphi)
365 bd[ii] = numpy.dot((-1*nvec*rad),(matrix0 + matrix1))
366
367 bmod = numpy.sqrt(bn**2. + be**2. + bd**2.)
368 btheta = numpy.arctan(bd/numpy.sqrt(be**2. + bn**2.))*180/numpy.pi
369 balpha = numpy.arctan(be/bn)*180./numpy.pi
370
371 #bn : north
372 #be : east
373 #bn : radial
374 #bmod : module
375
376
377 return bn, be, bd, bmod, btheta, balpha
378
379 def str2num(self, datum):
380 try:
381 return int(datum)
382 except:
383 try:
384 return float(datum)
385 except:
386 return datum
387
388 def __readIGRFfile(self, filename):
389 list_years=[]
390 for i in range(1,26):
391 list_years.append(1895.0 + i*5)
392
393 epochs=list_years
394 epochs.append(epochs[-1]+5)
395 nepochs = numpy.shape(epochs)
396
397 gg = numpy.zeros((13,14,nepochs[0]),dtype=float)
398 hh = numpy.zeros((13,14,nepochs[0]),dtype=float)
399
400 coeffs_file=open(filename)
401 lines=coeffs_file.readlines()
402
403 coeffs_file.close()
404
405 for line in lines:
406 items = line.split()
407 g_h = items[0]
408 n = self.str2num(items[1])
409 m = self.str2num(items[2])
410
411 coeffs = items[3:]
412
413 for i in range(len(coeffs)-1):
414 coeffs[i] = self.str2num(coeffs[i])
415
416 #coeffs = numpy.array(coeffs)
417 ncoeffs = numpy.shape(coeffs)[0]
418
419 if g_h == 'g':
420 # print n," g ",m
421 gg[n-1,m,:]=coeffs
422 elif g_h=='h':
423 # print n," h ",m
424 hh[n-1,m,:]=coeffs
425 # else :
426 # continue
427
428 # Ultimo Reordenamiento para almacenar .
429 gg[:,:,nepochs[0]-1] = gg[:,:,nepochs[0]-2] + 5*gg[:,:,nepochs[0]-1]
430 hh[:,:,nepochs[0]-1] = hh[:,:,nepochs[0]-2] + 5*hh[:,:,nepochs[0]-1]
431
432 # return numpy.array([gg,hh])
433 periods = numpy.array(epochs)
434 g = gg
435 h = hh
436 return periods, g, h
437
438
439 def __readIGRFcoeff(self,filename="igrf10coeffs.dat"):
440 """
441 __readIGRFcoeff reads the coefficients from a binary file which is located in the
442 folder "resource."
443
444 Parameter
445 ---------
446 filename = A string to specify the name of the file which contains thec coeffs. The
447 default value is "igrf10coeffs.dat"
448
449 Return
450 ------
451 periods = A lineal array giving...
452 g1 =
453 h1 =
454
455 Modification History
456 --------------------
457 Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009.
458 """
459
460 base_path = os.path.dirname(os.path.abspath(__file__))
461 filename = os.path.join(base_path, "igrf13coeffs.txt")
462
463 period_v, g_v, h_v = self.__readIGRFfile(filename)
464 g2 = numpy.zeros((14,14,26))
465 h2 = numpy.zeros((14,14,26))
466 g2[1:14,:,:] = g_v
467 h2[1:14,:,:] = h_v
468
469 g = numpy.transpose(g2, (2,0,1))
470 h = numpy.transpose(h2, (2,0,1))
471 periods = period_v.copy()
472
473 return periods, g, h
474
475 def rotvector(self,vector,axis=1,ang=0):
476 """
477 rotvector function returns the new vector generated rotating the rectagular coords.
478
479 Parameters
480 ----------
481 vector = A lineal 3-elements array (x,y,z).
482 axis = A integer to specify the axis used to rotate the coord systems. The default
483 value is 1.
484 axis = 1 -> Around "x"
485 axis = 2 -> Around "y"
486 axis = 3 -> Around "z"
487 ang = Angle of rotation (in radians). The default value is zero.
488
489 Return
490 ------
491 rotvector = A lineal array of 3 elements giving the new coordinates.
492
493 Modification History
494 --------------------
495 Converted to Python by Freddy R. Galindo, ROJ, 01 October 2009.
496 """
497
498 if axis==1:
499 t = [[1,0,0],[0,numpy.cos(ang),numpy.sin(ang)],[0,-numpy.sin(ang),numpy.cos(ang)]]
500 elif axis==2:
501 t = [[numpy.cos(ang),0,-numpy.sin(ang)],[0,1,0],[numpy.sin(ang),0,numpy.cos(ang)]]
502 elif axis==3:
503 t = [[numpy.cos(ang),numpy.sin(ang),0],[-numpy.sin(ang),numpy.cos(ang),0],[0,0,1]]
504
505 rotvector = numpy.array(numpy.dot(numpy.array(t),numpy.array(vector)))
506
507 return rotvector
Show file before | Show file after | Hide comments
@@ -0,0 +1,61
1 """
2 The module MISC_ROUTINES gathers classes and functions which are useful for daily processing. As an
3 example we have conversion factor or universal constants.
4
5 MODULES CALLED:
6 NUMPY, SYS
7
8 MODIFICATION HISTORY:
9 Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ, 21 October 2009.
10 """
11
12 import numpy
13 import sys
14
15 class CoFactors():
16 """
17 CoFactor class used to call pre-defined conversion factor (e.g. degree to radian). The cu-
18 The current available factor are:
19
20 d2r = degree to radian.
21 s2r = seconds to radian?, degree to arcsecond.?
22 h2r = hour to radian.
23 h2d = hour to degree
24 """
25
26 d2r = numpy.pi/180.
27 s2r = numpy.pi/(180.*3600.)
28 h2r = numpy.pi/12.
29 h2d = 15.
30
31
32 class Vector:
33 """
34 direction = 0 Polar to rectangular; direction=1 rectangular to polar
35 """
36 def __init__(self,vect,direction=0):
37 nsize = numpy.size(vect)
38 if nsize <= 3:
39 vect = vect.reshape(1,nsize)
40
41 self.vect = vect
42 self.dirc = direction
43
44
45
46 def Polar2Rect(self):
47 if self.dirc == 0:
48 jvect = self.vect*numpy.pi/180.
49 mmx = numpy.cos(jvect[:,1])*numpy.sin(jvect[:,0])
50 mmy = numpy.cos(jvect[:,1])*numpy.cos(jvect[:,0])
51 mmz = numpy.sin(jvect[:,1])
52 mm = numpy.array([mmx,mmy,mmz]).transpose()
53
54 elif self.dirc == 1:
55 mm = [numpy.arctan2(self.vect[:,0],self.vect[:,1]),numpy.arcsin(self.vect[:,2])]
56 mm = numpy.array(mm)*180./numpy.pi
57
58 return mm
59
60
61 No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,430
1 """
2 The TIME_CONVERSIONS.py module gathers classes and functions for time system transformations
3 (e.g. between seconds from 1970 to datetime format).
4
5 MODULES CALLED:
6 NUMPY, TIME, DATETIME, CALENDAR
7
8 MODIFICATION HISTORY:
9 Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Aug 13, 2009.
10 """
11
12 import numpy
13 import time
14 from datetime import datetime as dt
15 import calendar
16
17 class Time:
18 """
19 time(year,month,dom,hour,min,secs)
20
21 An object represents a date and time of certain event..
22
23 Parameters
24 ----------
25 YEAR = Number of the desired year. Year must be valid values from the civil calendar.
26 Years B.C.E must be represented as negative integers. Years in the common era are repre-
27 sented as positive integers. In particular, note that there is no year 0 in the civil
28 calendar. 1 B.C.E. (-1) is followed by 1 C.E. (1).
29
30 MONTH = Number of desired month (1=Jan, ..., 12=December).
31
32 DOM = Number of day of the month.
33
34 HOUR = Number of the hour of the day. By default hour=0
35
36 MINS = Number of the minute of the hour. By default min=0
37
38 SECS = Number of the second of the minute. By default secs=0.
39
40 Examples
41 --------
42 time_info = time(2008,9,30,12,30,00)
43
44 time_info = time(2008,9,30)
45 """
46
47 def __init__(self, year=None, month=None, dom=None, hour=0, mins=0, secs=0):
48 # If one the first three inputs are not defined, it takes the current date.
49 date = time.localtime()
50 if year==None:year=date[0]
51 if month==None:month=date[1]
52 if dom==None:dom=date[2]
53
54 # Converting to arrays
55 year = numpy.array([year]); month = numpy.array([month]); dom = numpy.array([dom])
56 hour = numpy.array([hour]); mins = numpy.array([mins]); secs = numpy.array([secs])
57
58 # Defining time information object.
59 self.year = numpy.atleast_1d(year)
60 self.month = numpy.atleast_1d(month)
61 self.dom = numpy.atleast_1d(dom)
62 self.hour = numpy.atleast_1d(hour)
63 self.mins = numpy.atleast_1d(mins)
64 self.secs = numpy.atleast_1d(secs)
65
66 def change2julday(self):
67 """
68 Converts a datetime to Julian days.
69 """
70
71 # Defining constants
72 greg = 2299171 # incorrect Julian day for Oct, 25, 1582.
73 min_calendar = -4716
74 max_calendar = 5000000
75
76 min_year = numpy.nanmin(self.year)
77 max_year = numpy.nanmax(self.year)
78 if (min_year<min_calendar) or (max_year>max_calendar):
79 print ("Value of Julian date is out of allowed range")
80 return -1
81
82 noyear = numpy.sum(self.year==0)
83 if noyear>0:
84 print ("There is no year zero in the civil calendar")
85 return -1
86
87 # Knowing if the year is less than 0.
88 bc = self.year<0
89
90 # Knowing if the month is less than March.
91 inJanFeb = self.month<=2
92
93 jy = self.year + bc - inJanFeb
94 jm = self.month + (1 + 12*inJanFeb)
95
96 # Computing Julian days.
97 jul= numpy.floor(365.25*jy) + numpy.floor(30.6001*jm) + (self.dom+1720995.0)
98
99 # Test whether to change to Gregorian Calendar
100 if numpy.min(jul) >= greg:
101 ja = numpy.int32(0.01*jy)
102 jul = jul + 2 - ja + numpy.int32(0.25*ja)
103 else:
104 gregchange = numpy.where(jul >= greg)
105 if gregchange[0].size>0:
106 ja = numpy.int32(0.01 + jy[gregchange])
107 jy[gregchange] = jy[gregchange] + 2 - ja + numpy.int32(0.25*ja)
108
109 # Determining machine-specific parameters affecting floating-point.
110 eps = 0.0 # Replace this line for a function to get precision.
111 eps = abs(jul)*0.0 > eps
112
113 jul = jul + (self.hour/24. -0.5) + (self.mins/1440.) + (self.secs/86400.) + eps
114
115 return jul[0]
116
117 def change2secs(self):
118 """
119 Converts datetime to number of seconds respect to 1970.
120 """
121
122 dtime = dt(self.year[0], self.month[0], self.dom[0])
123 return (dtime-dt(1970, 1, 1)).total_seconds()
124
125 year = self.year
126 if year.size>1: year = year[0]
127
128 month = self.month
129 if month.size>1: month = month[0]
130
131 dom = self.dom
132 if dom.size>1: dom = dom[0]
133
134 # Resizing hour, mins and secs if it was necessary.
135 hour = self.hour
136 if hour.size>1:hour = hour[0]
137 if hour.size==1:hour = numpy.resize(hour,year.size)
138
139 mins = self.mins
140 if mins.size>1:mins = mins[0]
141 if mins.size==1:mins = numpy.resize(mins,year.size)
142
143 secs = self.secs
144 if secs.size>1:secs = secs[0]
145 if secs.size==1:secs = numpy.resize(secs,year.size)
146
147 # Using time.mktime to compute seconds respect to 1970.
148 secs1970 = numpy.zeros(year.size)
149 for ii in numpy.arange(year.size):
150 secs1970[ii] = time.mktime((int(year[ii]),int(month[ii]),int(dom[ii]),\
151 int(hour[ii]),int(mins[ii]),int(secs[ii]),0,0,0))
152
153 secs1970 = numpy.int32(secs1970 - time.timezone)
154
155 return secs1970
156
157 def change2strdate(self,mode=1):
158 """
159 change2strdate method converts a date and time of certain event to date string. The
160 string format is like localtime (e.g. Fri Oct 9 15:00:19 2009).
161
162 Parameters
163 ----------
164 None.
165
166 Return
167 ------
168
169 Modification History
170 --------------------
171 Created by Freddy R. Galindo, ROJ, 09 October 2009.
172
173 """
174
175 secs = numpy.atleast_1d(self.change2secs())
176 strdate = []
177 for ii in numpy.arange(numpy.size(secs)):
178 secs_tmp = time.localtime(secs[ii] + time.timezone)
179 if mode==1:
180 strdate.append(time.strftime("%d-%b-%Y (%j) %H:%M:%S",secs_tmp))
181 elif mode==2:
182 strdate.append(time.strftime("%d-%b-%Y (%j)",secs_tmp))
183
184 strdate = numpy.array(strdate)
185
186 return strdate
187
188
189 class Secs:
190 """
191 secs(secs):
192
193 An object represents the number of seconds respect to 1970.
194
195 Parameters
196 ----------
197
198 SECS = A scalar or array giving the number of seconds respect to 1970.
199
200 Example:
201 --------
202 secs_info = secs(1251241373)
203
204 secs_info = secs([1251241373,1251241383,1251241393])
205 """
206 def __init__(self,secs):
207 self.secs = secs
208
209 def change2julday(self):
210 """
211 Convert seconds from 1970 to Julian days.
212 """
213
214 secs_1970 = time(1970,1,1,0,0,0).change2julday()
215
216 julian = self.secs/86400.0 + secs_1970
217
218 return julian
219
220 def change2time(self):
221 """
222 Converts seconds from 1970 to datetime.
223 """
224
225 secs1970 = numpy.atleast_1d(self.secs)
226
227 datetime = numpy.zeros((9,secs1970.size))
228 for ii in numpy.arange(secs1970.size):
229 tuple = time.gmtime(secs1970[ii])
230 datetime[0,ii] = tuple[0]
231 datetime[1,ii] = tuple[1]
232 datetime[2,ii] = tuple[2]
233 datetime[3,ii] = tuple[3]
234 datetime[4,ii] = tuple[4]
235 datetime[5,ii] = tuple[5]
236 datetime[6,ii] = tuple[6]
237 datetime[7,ii] = tuple[7]
238 datetime[8,ii] = tuple[8]
239
240 datetime = numpy.int32(datetime)
241
242 return datetime
243
244
245 class Julian:
246 """
247 julian(julian):
248
249 An object represents julian days.
250
251 Parameters
252 ----------
253
254 JULIAN = A scalar or array giving the julina days.
255
256 Example:
257 --------
258 julian_info = julian(2454740)
259
260 julian_info = julian([2454740,2454760,2454780])
261 """
262 def __init__(self,julian):
263 self.julian = numpy.atleast_1d(julian)
264
265 def change2time(self):
266 """
267 change2time method converts from julian day to calendar date and time.
268
269 Return
270 ------
271 year = An array giving the year of the desired julian day.
272 month = An array giving the month of the desired julian day.
273 dom = An array giving the day of the desired julian day.
274 hour = An array giving the hour of the desired julian day.
275 mins = An array giving the minute of the desired julian day.
276 secs = An array giving the second of the desired julian day.
277
278 Examples
279 --------
280 >> jd = 2455119.0
281 >> [yy,mo,dd,hh,mi,ss] = TimeTools.julian(jd).change2time()
282 >> print [yy,mo,dd,hh,mi,ss]
283 [2009] [10] [ 14.] [ 12.] [ 0.] [ 0.]
284
285 Modification history
286 --------------------
287 Translated from "Numerical Recipies in C", by William H. Press, Brian P. Flannery,
288 Saul A. Teukolsky, and William T. Vetterling. Cambridge University Press, 1988.
289 Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
290 """
291
292 min_julian = -1095
293 max_julian = 1827933925
294 if (numpy.min(self.julian) < min_julian) or (numpy.max(self.julian) > max_julian):
295 print ('Value of Julian date is out of allowed range.')
296 return None
297
298 # Beginning of Gregorian calendar
299 igreg = 2299161
300 julLong = numpy.floor(self.julian + 0.5)
301 minJul = numpy.min(julLong)
302
303 if (minJul >= igreg):
304 # All are Gregorian
305 jalpha = numpy.int32(((julLong - 1867216) - 0.25)/36524.25)
306 ja = julLong + 1 + jalpha - numpy.int32(0.25*jalpha)
307 else:
308 ja = julLong
309 gregChange = numpy.where(julLong >= igreg)
310 if gregChange[0].size>0:
311 jalpha = numpy.int32(((julLong[gregChange]-1867216) - 0.25)/36524.25)
312 ja[gregChange] = julLong[gregChange]+1+jalpha-numpy.int32(0.25*jalpha)
313
314 # clear memory.
315 jalpha = -1
316
317 jb = ja + 1524
318 jc = numpy.int32(6680. + ((jb-2439870)-122.1)/365.25)
319 jd = numpy.int32(365.*jc + (0.25*jc))
320 je = numpy.int32((jb - jd)/30.6001)
321
322 dom = jb - jd - numpy.int32(30.6001*je)
323 month = je - 1
324 month = ((month - 1) % 12) + 1
325 month = numpy.atleast_1d(month)
326 year = jc - 4715
327 year = year - (month > 2)*1
328 year = year - (year <= 0)*1
329 year = numpy.atleast_1d(year)
330
331 # Getting hours, minutes, seconds
332 fraction = self.julian + 0.5 - julLong
333 eps_0 = dom*0.0 + 1.0e-12
334 eps_1 = 1.0e-12*numpy.abs(julLong)
335 eps = (eps_0>eps_1)*eps_0 + (eps_0<=eps_1)*eps_1
336
337 hour_0 = dom*0 + 23
338 hour_2 = dom*0 + 0
339 hour_1 = numpy.floor(fraction*24.0 + eps)
340 hour = ((hour_1>hour_0)*23) + ((hour_1<=hour_0)*hour_1)
341 hour = ((hour_1<hour_2)*0) + ((hour_1>=hour_2)*hour_1)
342
343 fraction = fraction - (hour/24.0)
344 mins_0 = dom*0 + 59
345 mins_2 = dom*0 + 0
346 mins_1 = numpy.floor(fraction*1440.0 + eps)
347 mins = ((mins_1>mins_0)*59) + ((mins_1<=mins_0)*mins_1)
348 mins = ((mins_1<mins_2)*0) + ((mins_1>=mins_2)*mins_1)
349
350 secs_2 = dom*0 + 0
351 secs_1 = (fraction - mins/1440.0)*86400.0
352 secs = ((secs_1<secs_2)*0) + ((secs_1>=secs_2)*secs_1)
353
354 return year,month,dom,hour,mins,secs
355
356 def change2secs(self):
357 """
358 Converts from Julian days to seconds from 1970.
359 """
360
361 jul_1970 = Time(1970,1,1,0,0,0).change2julday()
362
363 secs = numpy.int32((self.julian - jul_1970)*86400)
364
365 return secs
366
367 def change2lst(self,longitude=-76.874369):
368 """
369 CT2LST converts from local civil time to local mean sideral time
370
371 longitude = The longitude in degrees (east of Greenwich) of the place for which
372 the local sideral time is desired, scalar. The Greenwich mean sideral time (GMST)
373 can be found by setting longitude=0.
374 """
375
376 # Useful constants, see Meus, p. 84
377 c = numpy.array([280.46061837, 360.98564736629, 0.000387933, 38710000.0])
378 jd2000 = 2451545.0
379 t0 = self.julian - jd2000
380 t = t0/36525.
381
382 # Computing GST in seconds
383 theta = c[0] + (c[1]*t0) + (t**2)*(c[2]-t/c[3])
384
385 # Computing LST in hours
386 lst = (theta + longitude)/15.0
387 neg = numpy.where(lst < 0.0)
388 if neg[0].size>0:lst[neg] = 24.0 + (lst[neg] % 24)
389 lst = lst % 24.0
390
391 return lst
392
393
394 class date2doy:
395 def __init__(self,year,month,day):
396 self.year = year
397 self.month = month
398 self.day = day
399
400 def change2doy(self):
401 if calendar.isleap(self.year) == True:
402 tfactor = 1
403 else:
404 tfactor = 2
405
406 day = self.day
407 month = self.month
408
409 doy = numpy.floor((275*month)/9.0) - (tfactor*numpy.floor((month+9)/12.0)) + day - 30
410
411 return numpy.int32(doy)
412
413
414 class Doy2Date:
415 def __init__(self,year,doy):
416 self.year = year
417 self.doy = doy
418
419 def change2date(self):
420 months = numpy.arange(12) + 1
421
422 first_dem = date2doy(self.year,months,1)
423 first_dem = first_dem.change2doy()
424
425 imm = numpy.where((self.doy - first_dem) > 0)
426
427 month = imm[0].size
428 dom = self.doy -first_dem[month - 1] + 1
429
430 return month, dom
Show file before | Show file after | Hide comments
@@ -1,709 +1,716
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
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 schainpy.model.data.jrodata import PlotterData
37 37 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
38 38 from schainpy.utils import log
39 39
40 40 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
41 41 blu_values = matplotlib.pyplot.get_cmap(
42 42 'seismic_r', 20)(numpy.arange(20))[10:15]
43 43 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
44 44 'jro', numpy.vstack((blu_values, jet_values)))
45 45 matplotlib.pyplot.register_cmap(cmap=ncmap)
46 46
47 47 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
48 48 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
49 49
50 50 EARTH_RADIUS = 6.3710e3
51 51
52 52 def ll2xy(lat1, lon1, lat2, lon2):
53 53
54 54 p = 0.017453292519943295
55 55 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
56 56 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
57 57 r = 12742 * numpy.arcsin(numpy.sqrt(a))
58 58 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
59 59 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
60 60 theta = -theta + numpy.pi/2
61 61 return r*numpy.cos(theta), r*numpy.sin(theta)
62 62
63 63
64 64 def km2deg(km):
65 65 '''
66 66 Convert distance in km to degrees
67 67 '''
68 68
69 69 return numpy.rad2deg(km/EARTH_RADIUS)
70 70
71 71
72 72 def figpause(interval):
73 73 backend = plt.rcParams['backend']
74 74 if backend in matplotlib.rcsetup.interactive_bk:
75 75 figManager = matplotlib._pylab_helpers.Gcf.get_active()
76 76 if figManager is not None:
77 77 canvas = figManager.canvas
78 78 if canvas.figure.stale:
79 79 canvas.draw()
80 80 try:
81 81 canvas.start_event_loop(interval)
82 82 except:
83 83 pass
84 84 return
85 85
86 86 def popup(message):
87 87 '''
88 88 '''
89 89
90 90 fig = plt.figure(figsize=(12, 8), facecolor='r')
91 91 text = '\n'.join([s.strip() for s in message.split(':')])
92 92 fig.text(0.01, 0.5, text, ha='left', va='center',
93 93 size='20', weight='heavy', color='w')
94 94 fig.show()
95 95 figpause(1000)
96 96
97 97
98 98 class Throttle(object):
99 99 '''
100 100 Decorator that prevents a function from being called more than once every
101 101 time period.
102 102 To create a function that cannot be called more than once a minute, but
103 103 will sleep until it can be called:
104 104 @Throttle(minutes=1)
105 105 def foo():
106 106 pass
107 107
108 108 for i in range(10):
109 109 foo()
110 110 print "This function has run %s times." % i
111 111 '''
112 112
113 113 def __init__(self, seconds=0, minutes=0, hours=0):
114 114 self.throttle_period = datetime.timedelta(
115 115 seconds=seconds, minutes=minutes, hours=hours
116 116 )
117 117
118 118 self.time_of_last_call = datetime.datetime.min
119 119
120 120 def __call__(self, fn):
121 121 @wraps(fn)
122 122 def wrapper(*args, **kwargs):
123 123 coerce = kwargs.pop('coerce', None)
124 124 if coerce:
125 125 self.time_of_last_call = datetime.datetime.now()
126 126 return fn(*args, **kwargs)
127 127 else:
128 128 now = datetime.datetime.now()
129 129 time_since_last_call = now - self.time_of_last_call
130 130 time_left = self.throttle_period - time_since_last_call
131 131
132 132 if time_left > datetime.timedelta(seconds=0):
133 133 return
134 134
135 135 self.time_of_last_call = datetime.datetime.now()
136 136 return fn(*args, **kwargs)
137 137
138 138 return wrapper
139 139
140 140 def apply_throttle(value):
141 141
142 142 @Throttle(seconds=value)
143 143 def fnThrottled(fn):
144 144 fn()
145 145
146 146 return fnThrottled
147 147
148 148
149 149 @MPDecorator
150 150 class Plot(Operation):
151 151 """Base class for Schain plotting operations
152 152
153 153 This class should never be use directtly you must subclass a new operation,
154 154 children classes must be defined as follow:
155 155
156 156 ExamplePlot(Plot):
157 157
158 158 CODE = 'code'
159 159 colormap = 'jet'
160 160 plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
161 161
162 162 def setup(self):
163 163 pass
164 164
165 165 def plot(self):
166 166 pass
167 167
168 168 """
169 169
170 170 CODE = 'Figure'
171 171 colormap = 'jet'
172 172 bgcolor = 'white'
173 173 buffering = True
174 174 __missing = 1E30
175 175
176 176 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
177 177 'showprofile']
178 178
179 179 def __init__(self):
180 180
181 181 Operation.__init__(self)
182 182 self.isConfig = False
183 183 self.isPlotConfig = False
184 184 self.save_time = 0
185 185 self.sender_time = 0
186 186 self.data = None
187 187 self.firsttime = True
188 188 self.sender_queue = deque(maxlen=10)
189 189 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
190 190
191 191 def __fmtTime(self, x, pos):
192 192 '''
193 193 '''
194 194 if self.t_units == "h_m":
195 195 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
196 196 if self.t_units == "h":
197 197 return '{}'.format(self.getDateTime(x).strftime('%H'))
198 198
199 199 def __setup(self, **kwargs):
200 200 '''
201 201 Initialize variables
202 202 '''
203 203
204 204 self.figures = []
205 205 self.axes = []
206 206 self.cb_axes = []
207 207 self.pf_axes = []
208 208 self.localtime = kwargs.pop('localtime', True)
209 209 self.show = kwargs.get('show', True)
210 210 self.save = kwargs.get('save', False)
211 211 self.save_period = kwargs.get('save_period', 0)
212 212 self.colormap = kwargs.get('colormap', self.colormap)
213 213 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
214 214 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
215 215 self.colormaps = kwargs.get('colormaps', None)
216 216 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
217 217 self.showprofile = kwargs.get('showprofile', False)
218 218 self.title = kwargs.get('wintitle', self.CODE.upper())
219 219 self.cb_label = kwargs.get('cb_label', None)
220 220 self.cb_labels = kwargs.get('cb_labels', None)
221 221 self.labels = kwargs.get('labels', None)
222 222 self.xaxis = kwargs.get('xaxis', 'frequency')
223 223 self.zmin = kwargs.get('zmin', None)
224 224 self.zmax = kwargs.get('zmax', None)
225 225 self.zlimits = kwargs.get('zlimits', None)
226 226 self.xmin = kwargs.get('xmin', None)
227 227 self.xmax = kwargs.get('xmax', None)
228 228 self.xrange = kwargs.get('xrange', 12)
229 229 self.xscale = kwargs.get('xscale', None)
230 230 self.ymin = kwargs.get('ymin', None)
231 231 self.ymax = kwargs.get('ymax', None)
232 232 self.yscale = kwargs.get('yscale', None)
233 233 self.xlabel = kwargs.get('xlabel', None)
234 234 self.attr_time = kwargs.get('attr_time', 'utctime')
235 235 self.attr_data = kwargs.get('attr_data', 'data_param')
236 236 self.decimation = kwargs.get('decimation', None)
237 237 self.oneFigure = kwargs.get('oneFigure', True)
238 238 self.width = kwargs.get('width', None)
239 239 self.height = kwargs.get('height', None)
240 240 self.colorbar = kwargs.get('colorbar', True)
241 241 self.factors = kwargs.get('factors', range(18))
242 242 self.channels = kwargs.get('channels', None)
243 243 self.titles = kwargs.get('titles', [])
244 244 self.polar = False
245 245 self.type = kwargs.get('type', 'iq')
246 246 self.grid = kwargs.get('grid', False)
247 247 self.pause = kwargs.get('pause', False)
248 248 self.save_code = kwargs.get('save_code', self.CODE)
249 249 self.throttle = kwargs.get('throttle', 0)
250 250 self.exp_code = kwargs.get('exp_code', None)
251 251 self.server = kwargs.get('server', False)
252 252 self.sender_period = kwargs.get('sender_period', 60)
253 253 self.tag = kwargs.get('tag', '')
254 254 self.height_index = kwargs.get('height_index', [])
255 255 self.__throttle_plot = apply_throttle(self.throttle)
256 256 code = self.attr_data if self.attr_data else self.CODE
257 257 self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
258 258 self.tmin = kwargs.get('tmin', None)
259 259 self.t_units = kwargs.get('t_units', "h_m")
260 260 self.selectedHeightsList = kwargs.get('selectedHeightsList', [])
261 self.extFile = kwargs.get('filename', None)
262 self.bFieldList = kwargs.get('bField', [])
263 self.celestialList = kwargs.get('celestial', [])
264
265 if isinstance(self.bFieldList, int):
266 self.bFieldList = [self.bFieldList]
267
261 268 if isinstance(self.selectedHeightsList, int):
262 269 self.selectedHeightsList = [self.selectedHeightsList]
263 270
264 271 if self.server:
265 272 if not self.server.startswith('tcp://'):
266 273 self.server = 'tcp://{}'.format(self.server)
267 274 log.success(
268 275 'Sending to server: {}'.format(self.server),
269 276 self.name
270 277 )
271 278
272 279 if isinstance(self.attr_data, str):
273 280 self.attr_data = [self.attr_data]
274 281
275 282 def __setup_plot(self):
276 283 '''
277 284 Common setup for all figures, here figures and axes are created
278 285 '''
279 286
280 287 self.setup()
281 288
282 289 self.time_label = 'LT' if self.localtime else 'UTC'
283 290
284 291 if self.width is None:
285 292 self.width = 8
286 293
287 294 self.figures = []
288 295 self.axes = []
289 296 self.cb_axes = []
290 297 self.pf_axes = []
291 298 self.cmaps = []
292 299
293 300 size = '15%' if self.ncols == 1 else '30%'
294 301 pad = '4%' if self.ncols == 1 else '8%'
295 302
296 303 if self.oneFigure:
297 304 if self.height is None:
298 305 self.height = 1.4 * self.nrows + 1
299 306 fig = plt.figure(figsize=(self.width, self.height),
300 307 edgecolor='k',
301 308 facecolor='w')
302 309 self.figures.append(fig)
303 310 for n in range(self.nplots):
304 311 ax = fig.add_subplot(self.nrows, self.ncols,
305 312 n + 1, polar=self.polar)
306 313 ax.tick_params(labelsize=8)
307 314 ax.firsttime = True
308 315 ax.index = 0
309 316 ax.press = None
310 317 self.axes.append(ax)
311 318 if self.showprofile:
312 319 cax = self.__add_axes(ax, size=size, pad=pad)
313 320 cax.tick_params(labelsize=8)
314 321 self.pf_axes.append(cax)
315 322 else:
316 323 if self.height is None:
317 324 self.height = 3
318 325 for n in range(self.nplots):
319 326 fig = plt.figure(figsize=(self.width, self.height),
320 327 edgecolor='k',
321 328 facecolor='w')
322 329 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
323 330 ax.tick_params(labelsize=8)
324 331 ax.firsttime = True
325 332 ax.index = 0
326 333 ax.press = None
327 334 self.figures.append(fig)
328 335 self.axes.append(ax)
329 336 if self.showprofile:
330 337 cax = self.__add_axes(ax, size=size, pad=pad)
331 338 cax.tick_params(labelsize=8)
332 339 self.pf_axes.append(cax)
333 340
334 341 for n in range(self.nrows):
335 342 if self.colormaps is not None:
336 343 cmap = plt.get_cmap(self.colormaps[n])
337 344 else:
338 345 cmap = plt.get_cmap(self.colormap)
339 346 cmap.set_bad(self.bgcolor, 1.)
340 347 self.cmaps.append(cmap)
341 348
342 349 def __add_axes(self, ax, size='30%', pad='8%'):
343 350 '''
344 351 Add new axes to the given figure
345 352 '''
346 353 divider = make_axes_locatable(ax)
347 354 nax = divider.new_horizontal(size=size, pad=pad)
348 355 ax.figure.add_axes(nax)
349 356 return nax
350 357
351 358 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
352 359 '''
353 360 Create a masked array for missing data
354 361 '''
355 362 if x_buffer.shape[0] < 2:
356 363 return x_buffer, y_buffer, z_buffer
357 364
358 365 deltas = x_buffer[1:] - x_buffer[0:-1]
359 366 x_median = numpy.median(deltas)
360 367
361 368 index = numpy.where(deltas > 5 * x_median)
362 369
363 370 if len(index[0]) != 0:
364 371 z_buffer[::, index[0], ::] = self.__missing
365 372 z_buffer = numpy.ma.masked_inside(z_buffer,
366 373 0.99 * self.__missing,
367 374 1.01 * self.__missing)
368 375
369 376 return x_buffer, y_buffer, z_buffer
370 377
371 378 def decimate(self):
372 379
373 380 # dx = int(len(self.x)/self.__MAXNUMX) + 1
374 381 dy = int(len(self.y) / self.decimation) + 1
375 382
376 383 # x = self.x[::dx]
377 384 x = self.x
378 385 y = self.y[::dy]
379 386 z = self.z[::, ::, ::dy]
380 387
381 388 return x, y, z
382 389
383 390 def format(self):
384 391 '''
385 392 Set min and max values, labels, ticks and titles
386 393 '''
387 394
388 395 for n, ax in enumerate(self.axes):
389 396 if ax.firsttime:
390 397 if self.xaxis != 'time':
391 398 xmin = self.xmin
392 399 xmax = self.xmax
393 400 else:
394 401 xmin = self.tmin
395 402 xmax = self.tmin + self.xrange*60*60
396 403 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
397 404 if self.t_units == "h_m":
398 405 ax.xaxis.set_major_locator(LinearLocator(9))
399 406 if self.t_units == "h":
400 407 ax.xaxis.set_major_locator(LinearLocator(int((xmax-xmin)/3600)+1))
401 408 ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
402 409 ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
403 410 ax.set_facecolor(self.bgcolor)
404 411 if self.xscale:
405 412 ax.xaxis.set_major_formatter(FuncFormatter(
406 413 lambda x, pos: '{0:g}'.format(x*self.xscale)))
407 414 if self.yscale:
408 415 ax.yaxis.set_major_formatter(FuncFormatter(
409 416 lambda x, pos: '{0:g}'.format(x*self.yscale)))
410 417 if self.xlabel is not None:
411 418 ax.set_xlabel(self.xlabel)
412 419 if self.ylabel is not None:
413 420 ax.set_ylabel(self.ylabel)
414 421 if self.showprofile:
415 422 self.pf_axes[n].set_ylim(ymin, ymax)
416 423 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
417 424 self.pf_axes[n].set_xlabel('dB')
418 425 self.pf_axes[n].grid(b=True, axis='x')
419 426 [tick.set_visible(False)
420 427 for tick in self.pf_axes[n].get_yticklabels()]
421 428
422 429 if self.colorbar and not(hasattr(ax, 'cbar')) :
423 430 ax.cbar = plt.colorbar(
424 431 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
425 432 ax.cbar.ax.tick_params(labelsize=8)
426 433 ax.cbar.ax.press = None
427 434 if self.cb_label:
428 435 ax.cbar.set_label(self.cb_label, size=8)
429 436 elif self.cb_labels:
430 437 ax.cbar.set_label(self.cb_labels[n], size=8)
431 438 else:
432 439 ax.cbar = None
433 440 ax.set_xlim(xmin, xmax)
434 441 ax.set_ylim(ymin, ymax)
435 442 ax.firsttime = False
436 443 if self.grid:
437 444 ax.grid(True)
438 445
439 446 if not self.polar:
440 447 ax.set_title('{} {} {}'.format(
441 448 self.titles[n],
442 449 self.getDateTime(self.data.max_time).strftime(
443 450 '%Y-%m-%d %H:%M:%S'),
444 451 self.time_label),
445 452 size=8)
446 453 else:
447 454
448 455 ax.set_title('{}'.format(self.titles[n]), size=8)
449 456 ax.set_ylim(0, 90)
450 457 ax.set_yticks(numpy.arange(0, 90, 20))
451 458 ax.yaxis.labelpad = 40
452 459
453 460 if self.firsttime:
454 461 for n, fig in enumerate(self.figures):
455 462 fig.subplots_adjust(**self.plots_adjust)
456 463 self.firsttime = False
457 464
458 465 def clear_figures(self):
459 466 '''
460 467 Reset axes for redraw plots
461 468 '''
462 469
463 470 for ax in self.axes+self.pf_axes+self.cb_axes:
464 471 ax.clear()
465 472 ax.firsttime = True
466 473 # #if hasattr(ax, 'cbar') and ax.cbar:
467 474 # # ax.cbar.remove()
468 475
469 476
470 477 def __plot(self):
471 478 '''
472 479 Main function to plot, format and save figures
473 480 '''
474 481
475 482 self.plot()
476 483 self.format()
477 484
478 485 for n, fig in enumerate(self.figures):
479 486 if self.nrows == 0 or self.nplots == 0:
480 487 log.warning('No data', self.name)
481 488 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
482 489 fig.canvas.manager.set_window_title(self.CODE)
483 490 continue
484 491
485 492 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
486 493 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
487 494
488 495 fig.canvas.draw()
489 496 if self.show:
490 497 fig.show()
491 498 figpause(0.01)
492 499
493 500 if self.save:
494 501 self.save_figure(n)
495 502
496 503 if self.server:
497 504 self.send_to_server()
498 505
499 506 def __update(self, dataOut, timestamp):
500 507 '''
501 508 '''
502 509
503 510 metadata = {
504 511 'yrange': dataOut.heightList,
505 512 'interval': dataOut.timeInterval,
506 513 'channels': dataOut.channelList
507 514 }
508 515 data, meta = self.update(dataOut)
509 516 metadata.update(meta)
510 517 self.data.update(data, timestamp, metadata)
511 518
512 519 def save_figure(self, n):
513 520 '''
514 521 '''
515 522
516 523 if (self.data.max_time - self.save_time) <= self.save_period:
517 524 return
518 525
519 526 self.save_time = self.data.max_time
520 527
521 528 fig = self.figures[n]
522 529
523 530 if self.throttle == 0:
524 531 figname = os.path.join(
525 532 self.save,
526 533 self.save_code,
527 534 '{}_{}.png'.format(
528 535 self.save_code,
529 536 self.getDateTime(self.data.max_time).strftime(
530 537 '%Y%m%d_%H%M%S'
531 538 ),
532 539 )
533 540 )
534 541 log.log('Saving figure: {}'.format(figname), self.name)
535 542 if not os.path.isdir(os.path.dirname(figname)):
536 543 os.makedirs(os.path.dirname(figname))
537 544 fig.savefig(figname)
538 545
539 546 figname = os.path.join(
540 547 self.save,
541 548 '{}_{}.png'.format(
542 549 self.save_code,
543 550 self.getDateTime(self.data.min_time).strftime(
544 551 '%Y%m%d'
545 552 ),
546 553 )
547 554 )
548 555
549 556 log.log('Saving figure: {}'.format(figname), self.name)
550 557 if not os.path.isdir(os.path.dirname(figname)):
551 558 os.makedirs(os.path.dirname(figname))
552 559 fig.savefig(figname)
553 560
554 561 def send_to_server(self):
555 562 '''
556 563 '''
557 564
558 565 if self.exp_code == None:
559 566 log.warning('Missing `exp_code` skipping sending to server...')
560 567
561 568 last_time = self.data.max_time
562 569 interval = last_time - self.sender_time
563 570 if interval < self.sender_period:
564 571 return
565 572
566 573 self.sender_time = last_time
567 574
568 575 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
569 576 for attr in attrs:
570 577 value = getattr(self, attr)
571 578 if value:
572 579 if isinstance(value, (numpy.float32, numpy.float64)):
573 580 value = round(float(value), 2)
574 581 self.data.meta[attr] = value
575 582 if self.colormap == 'jet':
576 583 self.data.meta['colormap'] = 'Jet'
577 584 elif 'RdBu' in self.colormap:
578 585 self.data.meta['colormap'] = 'RdBu'
579 586 else:
580 587 self.data.meta['colormap'] = 'Viridis'
581 588 self.data.meta['interval'] = int(interval)
582 589
583 590 self.sender_queue.append(last_time)
584 591
585 592 while 1:
586 593 try:
587 594 tm = self.sender_queue.popleft()
588 595 except IndexError:
589 596 break
590 597 msg = self.data.jsonify(tm, self.save_code, self.plot_type)
591 598 self.socket.send_string(msg)
592 599 socks = dict(self.poll.poll(2000))
593 600 if socks.get(self.socket) == zmq.POLLIN:
594 601 reply = self.socket.recv_string()
595 602 if reply == 'ok':
596 603 log.log("Response from server ok", self.name)
597 604 time.sleep(0.1)
598 605 continue
599 606 else:
600 607 log.warning(
601 608 "Malformed reply from server: {}".format(reply), self.name)
602 609 else:
603 610 log.warning(
604 611 "No response from server, retrying...", self.name)
605 612 self.sender_queue.appendleft(tm)
606 613 self.socket.setsockopt(zmq.LINGER, 0)
607 614 self.socket.close()
608 615 self.poll.unregister(self.socket)
609 616 self.socket = self.context.socket(zmq.REQ)
610 617 self.socket.connect(self.server)
611 618 self.poll.register(self.socket, zmq.POLLIN)
612 619 break
613 620
614 621 def setup(self):
615 622 '''
616 623 This method should be implemented in the child class, the following
617 624 attributes should be set:
618 625
619 626 self.nrows: number of rows
620 627 self.ncols: number of cols
621 628 self.nplots: number of plots (channels or pairs)
622 629 self.ylabel: label for Y axes
623 630 self.titles: list of axes title
624 631
625 632 '''
626 633 raise NotImplementedError
627 634
628 635 def plot(self):
629 636 '''
630 637 Must be defined in the child class, the actual plotting method
631 638 '''
632 639 raise NotImplementedError
633 640
634 641 def update(self, dataOut):
635 642 '''
636 643 Must be defined in the child class, update self.data with new data
637 644 '''
638 645
639 646 data = {
640 647 self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
641 648 }
642 649 meta = {}
643 650
644 651 return data, meta
645 652
646 653 def run(self, dataOut, **kwargs):
647 654 '''
648 655 Main plotting routine
649 656 '''
650 657 if self.isConfig is False:
651 658 self.__setup(**kwargs)
652 659
653 660 if self.localtime:
654 661 self.getDateTime = datetime.datetime.fromtimestamp
655 662 else:
656 663 self.getDateTime = datetime.datetime.utcfromtimestamp
657 664
658 665 self.data.setup()
659 666 self.isConfig = True
660 667 if self.server:
661 668 self.context = zmq.Context()
662 669 self.socket = self.context.socket(zmq.REQ)
663 670 self.socket.connect(self.server)
664 671 self.poll = zmq.Poller()
665 672 self.poll.register(self.socket, zmq.POLLIN)
666 673
667 674 tm = getattr(dataOut, self.attr_time)
668 675
669 676 if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
670 677 self.clear_figures()
671 678 self.save_time = tm
672 679 self.tmin += self.xrange*60*60
673 680 self.__plot()
674 681 self.data.setup()
675 682
676 683
677 684
678 685 self.__update(dataOut, tm)
679 686
680 687 if self.isPlotConfig is False:
681 688 self.__setup_plot()
682 689 self.isPlotConfig = True
683 690 if self.xaxis == 'time':
684 691 dt = self.getDateTime(tm)
685 692 if self.xmin is None:
686 693 self.tmin = tm
687 694 self.xmin = dt.hour
688 695 minutes = (self.xmin-int(self.xmin)) * 60
689 696 seconds = (minutes - int(minutes)) * 60
690 697 self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
691 698 datetime.datetime(1970, 1, 1)).total_seconds()
692 699 if self.localtime:
693 700 self.tmin += time.timezone
694 701
695 702 if self.xmin is not None and self.xmax is not None:
696 703 self.xrange = self.xmax - self.xmin
697 704
698 705 if self.throttle == 0:
699 706 self.__plot()
700 707 else:
701 708 self.__throttle_plot(self.__plot)#, coerce=coerce)
702 709
703 710 def close(self):
704 711
705 712 if self.data and not self.data.flagNoData:
706 713 self.save_time = 0
707 714 self.__plot()
708 715 if self.data and not self.data.flagNoData and self.pause:
709 716 figpause(10)
Show file before | Show file after | Hide comments
@@ -1,1469 +1,1725
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 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13 from itertools import combinations
14 14 from matplotlib.ticker import LinearLocator
15 15
16 from schainpy.model.utils.BField import BField
17 from scipy.interpolate import splrep
18 from scipy.interpolate import splev
19
16 20 from matplotlib import __version__ as plt_version
17 21
18 22 if plt_version >='3.3.4':
19 23 EXTRA_POINTS = 0
20 24 else:
21 25 EXTRA_POINTS = 1
22 26
23 27 class SpectraPlot(Plot):
24 28 '''
25 29 Plot for Spectra data
26 30 '''
27 31
28 32 CODE = 'spc'
29 33 colormap = 'jet'
30 34 plot_type = 'pcolor'
31 35 buffering = False
32 36 channelList = []
33 37 elevationList = []
34 38 azimuthList = []
35 39
36 40 def setup(self):
37 41
38 42 self.nplots = len(self.data.channels)
39 43 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
40 44 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
41 45 self.height = 3.4 * self.nrows
42 46
43 47 self.cb_label = 'dB'
44 48 if self.showprofile:
45 49 self.width = 5.2 * self.ncols
46 50 else:
47 51 self.width = 4.2* self.ncols
48 52 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
49 53 self.ylabel = 'Range [km]'
50 54
51 55
52 56 def update_list(self,dataOut):
53 57 if len(self.channelList) == 0:
54 58 self.channelList = dataOut.channelList
55 59 if len(self.elevationList) == 0:
56 60 self.elevationList = dataOut.elevationList
57 61 if len(self.azimuthList) == 0:
58 62 self.azimuthList = dataOut.azimuthList
59 63
60 64 def update(self, dataOut):
61 65
62 66 self.update_list(dataOut)
63 67 data = {}
64 68 meta = {}
65 69
66 70 #data['rti'] = dataOut.getPower()
67 71 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
68 72 noise = 10*numpy.log10(dataOut.getNoise()/norm)
69 73
74 <<<<<<< HEAD
75 =======
76
77
78 >>>>>>> 37cccf17c7b80521b59b978cb30e4ab2e6f37fce
70 79 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
71 80 for ch in range(dataOut.nChannels):
72 81 if hasattr(dataOut.normFactor,'ndim'):
73 82 if dataOut.normFactor.ndim > 1:
74 83 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
75 84
76 85 else:
77 86 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
78 87 else:
79 88 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
80 89
81 90
82 91 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
83 92 spc = 10*numpy.log10(z)
84 93
85 94 data['spc'] = spc
86 95 #print(spc[0].shape)
87 96 data['rti'] = spc.mean(axis=1)
88 97 data['noise'] = noise
89 98 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
90 99 if self.CODE == 'spc_moments':
91 100 data['moments'] = dataOut.moments
92 101
93 102 return data, meta
94 103
95 104 def plot(self):
96 105 if self.xaxis == "frequency":
97 106 x = self.data.xrange[0]
98 107 self.xlabel = "Frequency (kHz)"
99 108 elif self.xaxis == "time":
100 109 x = self.data.xrange[1]
101 110 self.xlabel = "Time (ms)"
102 111 else:
103 112 x = self.data.xrange[2]
104 113 self.xlabel = "Velocity (m/s)"
105 114
106 115 if self.CODE == 'spc_moments':
107 116 x = self.data.xrange[2]
108 117 self.xlabel = "Velocity (m/s)"
109 118
110 119 self.titles = []
111 120 y = self.data.yrange
112 121 self.y = y
113 122
114 123 data = self.data[-1]
115 124 z = data['spc']
116 125 #print(z.shape, x.shape, y.shape)
117 126 for n, ax in enumerate(self.axes):
118 127 noise = self.data['noise'][n][0]
119 128 #print(noise)
120 129 if self.CODE == 'spc_moments':
121 130 mean = data['moments'][n, 1]
122 131 if ax.firsttime:
123 132 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
124 133 self.xmin = self.xmin if self.xmin else -self.xmax
125 134 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
126 135 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
127 136 ax.plt = ax.pcolormesh(x, y, z[n].T,
128 137 vmin=self.zmin,
129 138 vmax=self.zmax,
130 139 cmap=plt.get_cmap(self.colormap)
131 140 )
132 141
133 142 if self.showprofile:
134 143 ax.plt_profile = self.pf_axes[n].plot(
135 144 data['rti'][n], y)[0]
136 145 # ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
137 146 # color="k", linestyle="dashed", lw=1)[0]
138 147 if self.CODE == 'spc_moments':
139 148 ax.plt_mean = ax.plot(mean, y, color='k')[0]
140 149 else:
141 150 ax.plt.set_array(z[n].T.ravel())
142 151 if self.showprofile:
143 152 ax.plt_profile.set_data(data['rti'][n], y)
144 153 #ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
145 154 if self.CODE == 'spc_moments':
146 155 ax.plt_mean.set_data(mean, y)
147 156 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
148 157 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
149 158 else:
150 159 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
151 160
152 161
153 162 class CrossSpectraPlot(Plot):
154 163
155 164 CODE = 'cspc'
156 165 colormap = 'jet'
157 166 plot_type = 'pcolor'
158 167 zmin_coh = None
159 168 zmax_coh = None
160 169 zmin_phase = None
161 170 zmax_phase = None
162 171 realChannels = None
163 172 crossPairs = None
164 173
165 174 def setup(self):
166 175
167 176 self.ncols = 4
168 177 self.nplots = len(self.data.pairs) * 2
169 178 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
170 179 self.width = 3.1 * self.ncols
171 180 self.height = 2.6 * self.nrows
172 181 self.ylabel = 'Range [km]'
173 182 self.showprofile = False
174 183 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
175 184
176 185 def update(self, dataOut):
177 186
178 187 data = {}
179 188 meta = {}
180 189
181 190 spc = dataOut.data_spc
182 191 cspc = dataOut.data_cspc
183 192 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
184 193 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
185 194 meta['pairs'] = rawPairs
186 195
187 196 if self.crossPairs == None:
188 197 self.crossPairs = dataOut.pairsList
189 198
190 199 tmp = []
191 200
192 201 for n, pair in enumerate(meta['pairs']):
193 202
194 203 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
195 204 coh = numpy.abs(out)
196 205 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
197 206 tmp.append(coh)
198 207 tmp.append(phase)
199 208
200 209 data['cspc'] = numpy.array(tmp)
201 210
202 211 return data, meta
203 212
204 213 def plot(self):
205 214
206 215 if self.xaxis == "frequency":
207 216 x = self.data.xrange[0]
208 217 self.xlabel = "Frequency (kHz)"
209 218 elif self.xaxis == "time":
210 219 x = self.data.xrange[1]
211 220 self.xlabel = "Time (ms)"
212 221 else:
213 222 x = self.data.xrange[2]
214 223 self.xlabel = "Velocity (m/s)"
215 224
216 225 self.titles = []
217 226
218 227 y = self.data.yrange
219 228 self.y = y
220 229
221 230 data = self.data[-1]
222 231 cspc = data['cspc']
223 232
224 233 for n in range(len(self.data.pairs)):
225 234
226 235 pair = self.crossPairs[n]
227 236
228 237 coh = cspc[n*2]
229 238 phase = cspc[n*2+1]
230 239 ax = self.axes[2 * n]
231 240
232 241 if ax.firsttime:
233 242 ax.plt = ax.pcolormesh(x, y, coh.T,
234 243 vmin=self.zmin_coh,
235 244 vmax=self.zmax_coh,
236 245 cmap=plt.get_cmap(self.colormap_coh)
237 246 )
238 247 else:
239 248 ax.plt.set_array(coh.T.ravel())
240 249 self.titles.append(
241 250 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
242 251
243 252 ax = self.axes[2 * n + 1]
244 253 if ax.firsttime:
245 254 ax.plt = ax.pcolormesh(x, y, phase.T,
246 255 vmin=-180,
247 256 vmax=180,
248 257 cmap=plt.get_cmap(self.colormap_phase)
249 258 )
250 259 else:
251 260 ax.plt.set_array(phase.T.ravel())
252 261
253 262 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
254 263
255 264
256 265 class RTIPlot(Plot):
257 266 '''
258 267 Plot for RTI data
259 268 '''
260 269
261 270 CODE = 'rti'
262 271 colormap = 'jet'
263 272 plot_type = 'pcolorbuffer'
264 273 titles = None
265 274 channelList = []
266 275 elevationList = []
267 276 azimuthList = []
268 277
269 278 def setup(self):
270 279 self.xaxis = 'time'
271 280 self.ncols = 1
272 281 #print("dataChannels ",self.data.channels)
273 282 self.nrows = len(self.data.channels)
274 283 self.nplots = len(self.data.channels)
275 284 self.ylabel = 'Range [km]'
276 285 #self.xlabel = 'Time'
277 286 self.cb_label = 'dB'
278 287 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
279 288 self.titles = ['{} Channel {}'.format(
280 289 self.CODE.upper(), x) for x in range(self.nplots)]
281 290
282 291 def update_list(self,dataOut):
283 292
284 293 if len(self.channelList) == 0:
285 294 self.channelList = dataOut.channelList
286 295 if len(self.elevationList) == 0:
287 296 self.elevationList = dataOut.elevationList
288 297 if len(self.azimuthList) == 0:
289 298 self.azimuthList = dataOut.azimuthList
290 299
291 300
292 301 def update(self, dataOut):
293 302 if len(self.channelList) == 0:
294 303 self.update_list(dataOut)
295 304 data = {}
296 305 meta = {}
297 306
298 307 data['rti'] = dataOut.getPower()
299 308
300 309 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
301 310 noise = 10*numpy.log10(dataOut.getNoise()/norm)
302 311 data['noise'] = noise
303 312
304 313 return data, meta
305 314
306 315 def plot(self):
307 316
308 317 self.x = self.data.times
309 318 self.y = self.data.yrange
310 319 #print(" x, y: ",self.x, self.y)
311 320 self.z = self.data[self.CODE]
312 321 self.z = numpy.array(self.z, dtype=float)
313 322 self.z = numpy.ma.masked_invalid(self.z)
314 323
315 324 try:
316 325 if self.channelList != None:
317 326 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
318 327 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
319 328 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
320 329 else:
321 330 self.titles = ['{} Channel {}'.format(
322 331 self.CODE.upper(), x) for x in self.channelList]
323 332 except:
324 333 if self.channelList.any() != None:
325 334 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
326 335 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
327 336 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
328 337 else:
329 338 self.titles = ['{} Channel {}'.format(
330 339 self.CODE.upper(), x) for x in self.channelList]
331 340
332 341 if self.decimation is None:
333 342 x, y, z = self.fill_gaps(self.x, self.y, self.z)
334 343 else:
335 344 x, y, z = self.fill_gaps(*self.decimate())
336 345
337 346 #dummy_var = self.axes #Extrañamente esto actualiza el valor axes
338 347 for n, ax in enumerate(self.axes):
339 348 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
340 349 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
341 350 data = self.data[-1]
342 351
343 352 if ax.firsttime:
344 353 if (n+1) == len(self.channelList):
345 354 ax.set_xlabel('Time')
346 355 ax.plt = ax.pcolormesh(x, y, z[n].T,
347 356 vmin=self.zmin,
348 357 vmax=self.zmax,
349 358 cmap=plt.get_cmap(self.colormap)
350 359 )
351 360 if self.showprofile:
352 361 ax.plot_profile = self.pf_axes[n].plot(data[self.CODE][n], self.y)[0]
353 362 if "noise" in self.data:
354 363
355 364 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
356 365 color="k", linestyle="dashed", lw=1)[0]
357 366 else:
358 367 ax.collections.remove(ax.collections[0])
359 368 ax.plt = ax.pcolormesh(x, y, z[n].T,
360 369 vmin=self.zmin,
361 370 vmax=self.zmax,
362 371 cmap=plt.get_cmap(self.colormap)
363 372 )
364 373 if self.showprofile:
365 374 ax.plot_profile.set_data(data[self.CODE][n], self.y)
366 375 if "noise" in self.data:
367 376 ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
368 377
369 378
370 379 class CoherencePlot(RTIPlot):
371 380 '''
372 381 Plot for Coherence data
373 382 '''
374 383
375 384 CODE = 'coh'
376 385 titles = None
377 386
378 387 def setup(self):
379 388 self.xaxis = 'time'
380 389 self.ncols = 1
381 390 self.nrows = len(self.data.pairs)
382 391 self.nplots = len(self.data.pairs)
383 392 self.ylabel = 'Range [km]'
384 393 #self.xlabel = 'Time'
385 394 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
386 395 if self.CODE == 'coh':
387 396 self.cb_label = ''
388 397 self.titles = [
389 398 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
390 399 else:
391 400 self.cb_label = 'Degrees'
392 401 self.titles = [
393 402 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
394 403
395 404
396 405 def update(self, dataOut):
397 406
398 407 data = {}
399 408 meta = {}
400 409 data['coh'] = dataOut.getCoherence()
401 410 meta['pairs'] = dataOut.pairsList
402 411
403 412
404 413 return data, meta
405 414
406 415 def plot(self):
407 416
408 417 self.x = self.data.times
409 418 self.y = self.data.yrange
410 419 self.z = self.data[self.CODE]
411 420
412 421 self.z = numpy.ma.masked_invalid(self.z)
413 422
414 423 if self.decimation is None:
415 424 x, y, z = self.fill_gaps(self.x, self.y, self.z)
416 425 else:
417 426 x, y, z = self.fill_gaps(*self.decimate())
418 427
419 428 for n, ax in enumerate(self.axes):
420 429 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
421 430 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
422 431 if ax.firsttime:
423 432 if (n+1) == len(self.channelList):
424 433 ax.set_xlabel('Time')
425 434 ax.plt = ax.pcolormesh(x, y, z[n].T,
426 435 vmin=self.zmin,
427 436 vmax=self.zmax,
428 437 cmap=plt.get_cmap(self.colormap)
429 438 )
430 439 if self.showprofile:
431 440 ax.plot_profile = self.pf_axes[n].plot(
432 441 self.data[self.CODE][n][-1], self.y)[0]
433 442 # ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
434 443 # color="k", linestyle="dashed", lw=1)[0]
435 444 else:
436 445 ax.collections.remove(ax.collections[0])
437 446 ax.plt = ax.pcolormesh(x, y, z[n].T,
438 447 vmin=self.zmin,
439 448 vmax=self.zmax,
440 449 cmap=plt.get_cmap(self.colormap)
441 450 )
442 451 if self.showprofile:
443 452 ax.plot_profile.set_data(self.data[self.CODE][n][-1], self.y)
444 453 # ax.plot_noise.set_data(numpy.repeat(
445 454 # self.data['noise'][n][-1], len(self.y)), self.y)
446 455
447 456
448 457
449 458 class PhasePlot(CoherencePlot):
450 459 '''
451 460 Plot for Phase map data
452 461 '''
453 462
454 463 CODE = 'phase'
455 464 colormap = 'seismic'
456 465
457 466 def update(self, dataOut):
458 467
459 468 data = {}
460 469 meta = {}
461 470 data['phase'] = dataOut.getCoherence(phase=True)
462 471 meta['pairs'] = dataOut.pairsList
463 472
464 473 return data, meta
465 474
466 475 class NoisePlot(Plot):
467 476 '''
468 477 Plot for noise
469 478 '''
470 479
471 480 CODE = 'noise'
472 481 plot_type = 'scatterbuffer'
473 482
474 483 def setup(self):
475 484 self.xaxis = 'time'
476 485 self.ncols = 1
477 486 self.nrows = 1
478 487 self.nplots = 1
479 488 self.ylabel = 'Intensity [dB]'
480 489 self.xlabel = 'Time'
481 490 self.titles = ['Noise']
482 491 self.colorbar = False
483 492 self.plots_adjust.update({'right': 0.85 })
484 493 #if not self.titles:
485 494 self.titles = ['Noise Plot']
486 495
487 496 def update(self, dataOut):
488 497
489 498 data = {}
490 499 meta = {}
491 500 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
492 501 noise = 10*numpy.log10(dataOut.getNoise())
493 502 noise = noise.reshape(dataOut.nChannels, 1)
494 503 data['noise'] = noise
495 504 meta['yrange'] = numpy.array([])
496 505
497 506 return data, meta
498 507
499 508 def plot(self):
500 509
501 510 x = self.data.times
502 511 xmin = self.data.min_time
503 512 xmax = xmin + self.xrange * 60 * 60
504 513 Y = self.data['noise']
505 514
506 515 if self.axes[0].firsttime:
507 516 if self.ymin is None: self.ymin = numpy.nanmin(Y) - 5
508 517 if self.ymax is None: self.ymax = numpy.nanmax(Y) + 5
509 518 for ch in self.data.channels:
510 519 y = Y[ch]
511 520 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
512 521 plt.legend(bbox_to_anchor=(1.18, 1.0))
513 522 else:
514 523 for ch in self.data.channels:
515 524 y = Y[ch]
516 525 self.axes[0].lines[ch].set_data(x, y)
517 526
518 527
519 528 class PowerProfilePlot(Plot):
520 529
521 530 CODE = 'pow_profile'
522 531 plot_type = 'scatter'
523 532
524 533 def setup(self):
525 534
526 535 self.ncols = 1
527 536 self.nrows = 1
528 537 self.nplots = 1
529 538 self.height = 4
530 539 self.width = 3
531 540 self.ylabel = 'Range [km]'
532 541 self.xlabel = 'Intensity [dB]'
533 542 self.titles = ['Power Profile']
534 543 self.colorbar = False
535 544
536 545 def update(self, dataOut):
537 546
538 547 data = {}
539 548 meta = {}
540 549 data[self.CODE] = dataOut.getPower()
541 550
542 551 return data, meta
543 552
544 553 def plot(self):
545 554
546 555 y = self.data.yrange
547 556 self.y = y
548 557
549 558 x = self.data[-1][self.CODE]
550 559
551 560 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
552 561 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
553 562
554 563 if self.axes[0].firsttime:
555 564 for ch in self.data.channels:
556 565 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
557 566 plt.legend()
558 567 else:
559 568 for ch in self.data.channels:
560 569 self.axes[0].lines[ch].set_data(x[ch], y)
561 570
562 571
563 572 class SpectraCutPlot(Plot):
564 573
565 574 CODE = 'spc_cut'
566 575 plot_type = 'scatter'
567 576 buffering = False
568 577 heights = []
569 578 channelList = []
570 579 maintitle = "Spectra Cuts"
571 580 flag_setIndex = False
572 581
573 582 def setup(self):
574 583
575 584 self.nplots = len(self.data.channels)
576 585 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
577 586 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
578 587 self.width = 4.5 * self.ncols + 2.5
579 588 self.height = 4.8 * self.nrows
580 589 self.ylabel = 'Power [dB]'
581 590 self.colorbar = False
582 591 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
583 592
584 593 if len(self.selectedHeightsList) > 0:
585 594 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
586 595
587 596
588 597
589 598 def update(self, dataOut):
590 599 if len(self.channelList) == 0:
591 600 self.channelList = dataOut.channelList
592 601
593 602 self.heights = dataOut.heightList
594 603 #print("sels: ",self.selectedHeightsList)
595 604 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
596 605
597 606 for sel_height in self.selectedHeightsList:
598 607 index_list = numpy.where(self.heights >= sel_height)
599 608 index_list = index_list[0]
600 609 self.height_index.append(index_list[0])
601 610 #print("sels i:"", self.height_index)
602 611 self.flag_setIndex = True
603 612 #print(self.height_index)
604 613 data = {}
605 614 meta = {}
606 615
607 616 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
608 617 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
609 618 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
610 619
611 620
612 621 z = []
613 622 for ch in range(dataOut.nChannels):
614 623 if hasattr(dataOut.normFactor,'shape'):
615 624 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
616 625 else:
617 626 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
618 627
619 628 z = numpy.asarray(z)
620 629 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
621 630 spc = 10*numpy.log10(z)
622 631
623 632
624 633 data['spc'] = spc - noise
625 634 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
626 635
627 636 return data, meta
628 637
629 638 def plot(self):
630 639 if self.xaxis == "frequency":
631 640 x = self.data.xrange[0][0:]
632 641 self.xlabel = "Frequency (kHz)"
633 642 elif self.xaxis == "time":
634 643 x = self.data.xrange[1]
635 644 self.xlabel = "Time (ms)"
636 645 else:
637 646 x = self.data.xrange[2]
638 647 self.xlabel = "Velocity (m/s)"
639 648
640 649 self.titles = []
641 650
642 651 y = self.data.yrange
643 652 z = self.data[-1]['spc']
644 653 #print(z.shape)
645 654 if len(self.height_index) > 0:
646 655 index = self.height_index
647 656 else:
648 657 index = numpy.arange(0, len(y), int((len(y))/9))
649 658 #print("inde x ", index, self.axes)
650 659
651 660 for n, ax in enumerate(self.axes):
652 661
653 662 if ax.firsttime:
654 663
655 664
656 665 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
657 666 self.xmin = self.xmin if self.xmin else -self.xmax
658 667 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
659 668 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
660 669
661 670
662 671 ax.plt = ax.plot(x, z[n, :, index].T)
663 672 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
664 673 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
665 674 ax.minorticks_on()
666 675 ax.grid(which='major', axis='both')
667 676 ax.grid(which='minor', axis='x')
668 677 else:
669 678 for i, line in enumerate(ax.plt):
670 679 line.set_data(x, z[n, :, index[i]])
671 680
672 681
673 682 self.titles.append('CH {}'.format(self.channelList[n]))
674 683 plt.suptitle(self.maintitle, fontsize=10)
675 684
676 685
677 686 class BeaconPhase(Plot):
678 687
679 688 __isConfig = None
680 689 __nsubplots = None
681 690
682 691 PREFIX = 'beacon_phase'
683 692
684 693 def __init__(self):
685 694 Plot.__init__(self)
686 695 self.timerange = 24*60*60
687 696 self.isConfig = False
688 697 self.__nsubplots = 1
689 698 self.counter_imagwr = 0
690 699 self.WIDTH = 800
691 700 self.HEIGHT = 400
692 701 self.WIDTHPROF = 120
693 702 self.HEIGHTPROF = 0
694 703 self.xdata = None
695 704 self.ydata = None
696 705
697 706 self.PLOT_CODE = BEACON_CODE
698 707
699 708 self.FTP_WEI = None
700 709 self.EXP_CODE = None
701 710 self.SUB_EXP_CODE = None
702 711 self.PLOT_POS = None
703 712
704 713 self.filename_phase = None
705 714
706 715 self.figfile = None
707 716
708 717 self.xmin = None
709 718 self.xmax = None
710 719
711 720 def getSubplots(self):
712 721
713 722 ncol = 1
714 723 nrow = 1
715 724
716 725 return nrow, ncol
717 726
718 727 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
719 728
720 729 self.__showprofile = showprofile
721 730 self.nplots = nplots
722 731
723 732 ncolspan = 7
724 733 colspan = 6
725 734 self.__nsubplots = 2
726 735
727 736 self.createFigure(id = id,
728 737 wintitle = wintitle,
729 738 widthplot = self.WIDTH+self.WIDTHPROF,
730 739 heightplot = self.HEIGHT+self.HEIGHTPROF,
731 740 show=show)
732 741
733 742 nrow, ncol = self.getSubplots()
734 743
735 744 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
736 745
737 746 def save_phase(self, filename_phase):
738 747 f = open(filename_phase,'w+')
739 748 f.write('\n\n')
740 749 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
741 750 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
742 751 f.close()
743 752
744 753 def save_data(self, filename_phase, data, data_datetime):
745 754 f=open(filename_phase,'a')
746 755 timetuple_data = data_datetime.timetuple()
747 756 day = str(timetuple_data.tm_mday)
748 757 month = str(timetuple_data.tm_mon)
749 758 year = str(timetuple_data.tm_year)
750 759 hour = str(timetuple_data.tm_hour)
751 760 minute = str(timetuple_data.tm_min)
752 761 second = str(timetuple_data.tm_sec)
753 762 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
754 763 f.close()
755 764
756 765 def plot(self):
757 766 log.warning('TODO: Not yet implemented...')
758 767
759 768 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
760 769 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
761 770 timerange=None,
762 771 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
763 772 server=None, folder=None, username=None, password=None,
764 773 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
765 774
766 775 if dataOut.flagNoData:
767 776 return dataOut
768 777
769 778 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
770 779 return
771 780
772 781 if pairsList == None:
773 782 pairsIndexList = dataOut.pairsIndexList[:10]
774 783 else:
775 784 pairsIndexList = []
776 785 for pair in pairsList:
777 786 if pair not in dataOut.pairsList:
778 787 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
779 788 pairsIndexList.append(dataOut.pairsList.index(pair))
780 789
781 790 if pairsIndexList == []:
782 791 return
783 792
784 793 # if len(pairsIndexList) > 4:
785 794 # pairsIndexList = pairsIndexList[0:4]
786 795
787 796 hmin_index = None
788 797 hmax_index = None
789 798
790 799 if hmin != None and hmax != None:
791 800 indexes = numpy.arange(dataOut.nHeights)
792 801 hmin_list = indexes[dataOut.heightList >= hmin]
793 802 hmax_list = indexes[dataOut.heightList <= hmax]
794 803
795 804 if hmin_list.any():
796 805 hmin_index = hmin_list[0]
797 806
798 807 if hmax_list.any():
799 808 hmax_index = hmax_list[-1]+1
800 809
801 810 x = dataOut.getTimeRange()
802 811
803 812 thisDatetime = dataOut.datatime
804 813
805 814 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
806 815 xlabel = "Local Time"
807 816 ylabel = "Phase (degrees)"
808 817
809 818 update_figfile = False
810 819
811 820 nplots = len(pairsIndexList)
812 821 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
813 822 phase_beacon = numpy.zeros(len(pairsIndexList))
814 823 for i in range(nplots):
815 824 pair = dataOut.pairsList[pairsIndexList[i]]
816 825 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
817 826 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
818 827 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
819 828 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
820 829 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
821 830
822 831 if dataOut.beacon_heiIndexList:
823 832 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
824 833 else:
825 834 phase_beacon[i] = numpy.average(phase)
826 835
827 836 if not self.isConfig:
828 837
829 838 nplots = len(pairsIndexList)
830 839
831 840 self.setup(id=id,
832 841 nplots=nplots,
833 842 wintitle=wintitle,
834 843 showprofile=showprofile,
835 844 show=show)
836 845
837 846 if timerange != None:
838 847 self.timerange = timerange
839 848
840 849 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
841 850
842 851 if ymin == None: ymin = 0
843 852 if ymax == None: ymax = 360
844 853
845 854 self.FTP_WEI = ftp_wei
846 855 self.EXP_CODE = exp_code
847 856 self.SUB_EXP_CODE = sub_exp_code
848 857 self.PLOT_POS = plot_pos
849 858
850 859 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
851 860 self.isConfig = True
852 861 self.figfile = figfile
853 862 self.xdata = numpy.array([])
854 863 self.ydata = numpy.array([])
855 864
856 865 update_figfile = True
857 866
858 867 #open file beacon phase
859 868 path = '%s%03d' %(self.PREFIX, self.id)
860 869 beacon_file = os.path.join(path,'%s.txt'%self.name)
861 870 self.filename_phase = os.path.join(figpath,beacon_file)
862 871 #self.save_phase(self.filename_phase)
863 872
864 873
865 874 #store data beacon phase
866 875 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
867 876
868 877 self.setWinTitle(title)
869 878
870 879
871 880 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
872 881
873 882 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
874 883
875 884 axes = self.axesList[0]
876 885
877 886 self.xdata = numpy.hstack((self.xdata, x[0:1]))
878 887
879 888 if len(self.ydata)==0:
880 889 self.ydata = phase_beacon.reshape(-1,1)
881 890 else:
882 891 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
883 892
884 893
885 894 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
886 895 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
887 896 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
888 897 XAxisAsTime=True, grid='both'
889 898 )
890 899
891 900 self.draw()
892 901
893 902 if dataOut.ltctime >= self.xmax:
894 903 self.counter_imagwr = wr_period
895 904 self.isConfig = False
896 905 update_figfile = True
897 906
898 907 self.save(figpath=figpath,
899 908 figfile=figfile,
900 909 save=save,
901 910 ftp=ftp,
902 911 wr_period=wr_period,
903 912 thisDatetime=thisDatetime,
904 913 update_figfile=update_figfile)
905 914
906 915 return dataOut
907 916
908 917 class NoiselessSpectraPlot(Plot):
909 918 '''
910 919 Plot for Spectra data, subtracting
911 920 the noise in all channels, using for
912 921 amisr-14 data
913 922 '''
914 923
915 924 CODE = 'noiseless_spc'
916 925 colormap = 'jet'
917 926 plot_type = 'pcolor'
918 927 buffering = False
919 928 channelList = []
920 929 last_noise = None
921 930
922 931 def setup(self):
923 932
924 933 self.nplots = len(self.data.channels)
925 934 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
926 935 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
927 936 self.height = 3.5 * self.nrows
928 937
929 938 self.cb_label = 'dB'
930 939 if self.showprofile:
931 940 self.width = 5.8 * self.ncols
932 941 else:
933 942 self.width = 4.8* self.ncols
934 943 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
935 944
936 945 self.ylabel = 'Range [km]'
937 946
938 947
939 948 def update_list(self,dataOut):
940 949 if len(self.channelList) == 0:
941 950 self.channelList = dataOut.channelList
942 951
943 952 def update(self, dataOut):
944 953
945 954 self.update_list(dataOut)
946 955 data = {}
947 956 meta = {}
948 957
949 958 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
950 959 n0 = (dataOut.getNoise()/norm)
951 960 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
952 961 noise = 10*numpy.log10(noise)
953 962
954 963 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
955 964 for ch in range(dataOut.nChannels):
956 965 if hasattr(dataOut.normFactor,'ndim'):
957 966 if dataOut.normFactor.ndim > 1:
958 967 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
959 968 else:
960 969 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
961 970 else:
962 971 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
963 972
964 973 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
965 974 spc = 10*numpy.log10(z)
966 975
967 976
968 977 data['spc'] = spc - noise
969 978 #print(spc.shape)
970 979 data['rti'] = spc.mean(axis=1)
971 980 data['noise'] = noise
972 981
973 982
974 983
975 984 # data['noise'] = noise
976 985 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
977 986
978 987 return data, meta
979 988
980 989 def plot(self):
981 990 if self.xaxis == "frequency":
982 991 x = self.data.xrange[0]
983 992 self.xlabel = "Frequency (kHz)"
984 993 elif self.xaxis == "time":
985 994 x = self.data.xrange[1]
986 995 self.xlabel = "Time (ms)"
987 996 else:
988 997 x = self.data.xrange[2]
989 998 self.xlabel = "Velocity (m/s)"
990 999
991 1000 self.titles = []
992 1001 y = self.data.yrange
993 1002 self.y = y
994 1003
995 1004 data = self.data[-1]
996 1005 z = data['spc']
997 1006
998 1007 for n, ax in enumerate(self.axes):
999 1008 #noise = data['noise'][n]
1000 1009
1001 1010 if ax.firsttime:
1002 1011 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1003 1012 self.xmin = self.xmin if self.xmin else -self.xmax
1004 1013 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1005 1014 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1006 1015 ax.plt = ax.pcolormesh(x, y, z[n].T,
1007 1016 vmin=self.zmin,
1008 1017 vmax=self.zmax,
1009 1018 cmap=plt.get_cmap(self.colormap)
1010 1019 )
1011 1020
1012 1021 if self.showprofile:
1013 1022 ax.plt_profile = self.pf_axes[n].plot(
1014 1023 data['rti'][n], y)[0]
1015 1024
1016 1025
1017 1026 else:
1018 1027 ax.plt.set_array(z[n].T.ravel())
1019 1028 if self.showprofile:
1020 1029 ax.plt_profile.set_data(data['rti'][n], y)
1021 1030
1022 1031
1023 1032 self.titles.append('CH {}'.format(self.channelList[n]))
1024 1033
1025 1034
1026 1035 class NoiselessRTIPlot(RTIPlot):
1027 1036 '''
1028 1037 Plot for RTI data
1029 1038 '''
1030 1039
1031 1040 CODE = 'noiseless_rti'
1032 1041 colormap = 'jet'
1033 1042 plot_type = 'pcolorbuffer'
1034 1043 titles = None
1035 1044 channelList = []
1036 1045 elevationList = []
1037 1046 azimuthList = []
1038 1047 last_noise = None
1039 1048
1040 1049 def setup(self):
1041 1050 self.xaxis = 'time'
1042 1051 self.ncols = 1
1043 1052 #print("dataChannels ",self.data.channels)
1044 1053 self.nrows = len(self.data.channels)
1045 1054 self.nplots = len(self.data.channels)
1046 1055 self.ylabel = 'Range [km]'
1047 1056 #self.xlabel = 'Time'
1048 1057 self.cb_label = 'dB'
1049 1058 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1050 1059 self.titles = ['{} Channel {}'.format(
1051 1060 self.CODE.upper(), x) for x in range(self.nplots)]
1052 1061
1053 1062 def update_list(self,dataOut):
1054 1063 if len(self.channelList) == 0:
1055 1064 self.channelList = dataOut.channelList
1056 1065 if len(self.elevationList) == 0:
1057 1066 self.elevationList = dataOut.elevationList
1058 1067 if len(self.azimuthList) == 0:
1059 1068 self.azimuthList = dataOut.azimuthList
1060 1069
1061 1070 def update(self, dataOut):
1062 1071 if len(self.channelList) == 0:
1063 1072 self.update_list(dataOut)
1064 1073
1065 1074 data = {}
1066 1075 meta = {}
1067 1076 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1068 1077 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt
1069 1078 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1070 1079 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1071 1080 data['noise'] = n0
1072 1081 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1073 1082 noiseless_data = dataOut.getPower() - noise
1074 1083
1075 1084 #print("power, noise:", dataOut.getPower(), n0)
1076 1085 #print(noise)
1077 1086 #print(noiseless_data)
1078 1087
1079 1088 data['noiseless_rti'] = noiseless_data
1080 1089
1081 1090 return data, meta
1082 1091
1083 1092 def plot(self):
1084 1093 from matplotlib import pyplot as plt
1085 1094 self.x = self.data.times
1086 1095 self.y = self.data.yrange
1087 1096 self.z = self.data['noiseless_rti']
1088 1097 self.z = numpy.array(self.z, dtype=float)
1089 1098 self.z = numpy.ma.masked_invalid(self.z)
1090 1099
1091 1100
1092 1101 try:
1093 1102 if self.channelList != None:
1094 1103 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1095 1104 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1096 1105 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1097 1106 else:
1098 1107 self.titles = ['{} Channel {}'.format(
1099 1108 self.CODE.upper(), x) for x in self.channelList]
1100 1109 except:
1101 1110 if self.channelList.any() != None:
1102 1111 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1103 1112 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1104 1113 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1105 1114 else:
1106 1115 self.titles = ['{} Channel {}'.format(
1107 1116 self.CODE.upper(), x) for x in self.channelList]
1108 1117
1109 1118
1110 1119 if self.decimation is None:
1111 1120 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1112 1121 else:
1113 1122 x, y, z = self.fill_gaps(*self.decimate())
1114 1123
1115 1124 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
1116 1125 #print("plot shapes ", z.shape, x.shape, y.shape)
1117 1126 #print(self.axes)
1118 1127 for n, ax in enumerate(self.axes):
1119 1128
1120 1129
1121 1130 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1122 1131 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1123 1132 data = self.data[-1]
1124 1133 if ax.firsttime:
1125 1134 if (n+1) == len(self.channelList):
1126 1135 ax.set_xlabel('Time')
1127 1136 ax.plt = ax.pcolormesh(x, y, z[n].T,
1128 1137 vmin=self.zmin,
1129 1138 vmax=self.zmax,
1130 1139 cmap=plt.get_cmap(self.colormap)
1131 1140 )
1132 1141 if self.showprofile:
1133 1142 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1134 1143
1135 1144 else:
1136 1145 ax.collections.remove(ax.collections[0])
1137 1146 ax.plt = ax.pcolormesh(x, y, z[n].T,
1138 1147 vmin=self.zmin,
1139 1148 vmax=self.zmax,
1140 1149 cmap=plt.get_cmap(self.colormap)
1141 1150 )
1142 1151 if self.showprofile:
1143 1152 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1144 1153 # if "noise" in self.data:
1145 1154 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1146 1155 # ax.plot_noise.set_data(data['noise'][n], self.y)
1147 1156
1148 1157
1149 1158 class OutliersRTIPlot(Plot):
1150 1159 '''
1151 1160 Plot for data_xxxx object
1152 1161 '''
1153 1162
1154 1163 CODE = 'outlier_rtc' # Range Time Counts
1155 1164 colormap = 'cool'
1156 1165 plot_type = 'pcolorbuffer'
1157 1166
1158 1167 def setup(self):
1159 1168 self.xaxis = 'time'
1160 1169 self.ncols = 1
1161 1170 self.nrows = self.data.shape('outlier_rtc')[0]
1162 1171 self.nplots = self.nrows
1163 1172 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1164 1173
1165 1174
1166 1175 if not self.xlabel:
1167 1176 self.xlabel = 'Time'
1168 1177
1169 1178 self.ylabel = 'Height [km]'
1170 1179 if not self.titles:
1171 1180 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1172 1181
1173 1182 def update(self, dataOut):
1174 1183
1175 1184 data = {}
1176 1185 data['outlier_rtc'] = dataOut.data_outlier
1177 1186
1178 1187 meta = {}
1179 1188
1180 1189 return data, meta
1181 1190
1182 1191 def plot(self):
1183 1192 # self.data.normalize_heights()
1184 1193 self.x = self.data.times
1185 1194 self.y = self.data.yrange
1186 1195 self.z = self.data['outlier_rtc']
1187 1196
1188 1197 #self.z = numpy.ma.masked_invalid(self.z)
1189 1198
1190 1199 if self.decimation is None:
1191 1200 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1192 1201 else:
1193 1202 x, y, z = self.fill_gaps(*self.decimate())
1194 1203
1195 1204 for n, ax in enumerate(self.axes):
1196 1205
1197 1206 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1198 1207 self.z[n])
1199 1208 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1200 1209 self.z[n])
1201 1210 data = self.data[-1]
1202 1211 if ax.firsttime:
1203 1212 if self.zlimits is not None:
1204 1213 self.zmin, self.zmax = self.zlimits[n]
1205 1214
1206 1215 ax.plt = ax.pcolormesh(x, y, z[n].T,
1207 1216 vmin=self.zmin,
1208 1217 vmax=self.zmax,
1209 1218 cmap=self.cmaps[n]
1210 1219 )
1211 1220 if self.showprofile:
1212 1221 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1213 1222 self.pf_axes[n].set_xlabel('')
1214 1223 else:
1215 1224 if self.zlimits is not None:
1216 1225 self.zmin, self.zmax = self.zlimits[n]
1217 1226 ax.collections.remove(ax.collections[0])
1218 1227 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1219 1228 vmin=self.zmin,
1220 1229 vmax=self.zmax,
1221 1230 cmap=self.cmaps[n]
1222 1231 )
1223 1232 if self.showprofile:
1224 1233 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1225 1234 self.pf_axes[n].set_xlabel('')
1226 1235
1227 1236 class NIncohIntRTIPlot(Plot):
1228 1237 '''
1229 1238 Plot for data_xxxx object
1230 1239 '''
1231 1240
1232 1241 CODE = 'integrations_rtc' # Range Time Counts
1233 1242 colormap = 'BuGn'
1234 1243 plot_type = 'pcolorbuffer'
1235 1244
1236 1245 def setup(self):
1237 1246 self.xaxis = 'time'
1238 1247 self.ncols = 1
1239 1248 self.nrows = self.data.shape('integrations_rtc')[0]
1240 1249 self.nplots = self.nrows
1241 1250 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1242 1251
1243 1252
1244 1253 if not self.xlabel:
1245 1254 self.xlabel = 'Time'
1246 1255
1247 1256 self.ylabel = 'Height [km]'
1248 1257 if not self.titles:
1249 1258 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1250 1259
1251 1260 def update(self, dataOut):
1252 1261
1253 1262 data = {}
1254 1263 data['integrations_rtc'] = dataOut.nIncohInt
1255 1264
1256 1265 meta = {}
1257 1266
1258 1267 return data, meta
1259 1268
1260 1269 def plot(self):
1261 1270 # self.data.normalize_heights()
1262 1271 self.x = self.data.times
1263 1272 self.y = self.data.yrange
1264 1273 self.z = self.data['integrations_rtc']
1265 1274
1266 1275 #self.z = numpy.ma.masked_invalid(self.z)
1267 1276
1268 1277 if self.decimation is None:
1269 1278 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1270 1279 else:
1271 1280 x, y, z = self.fill_gaps(*self.decimate())
1272 1281
1273 1282 for n, ax in enumerate(self.axes):
1274 1283
1275 1284 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1276 1285 self.z[n])
1277 1286 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1278 1287 self.z[n])
1279 1288 data = self.data[-1]
1280 1289 if ax.firsttime:
1281 1290 if self.zlimits is not None:
1282 1291 self.zmin, self.zmax = self.zlimits[n]
1283 1292
1284 1293 ax.plt = ax.pcolormesh(x, y, z[n].T,
1285 1294 vmin=self.zmin,
1286 1295 vmax=self.zmax,
1287 1296 cmap=self.cmaps[n]
1288 1297 )
1289 1298 if self.showprofile:
1290 1299 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1291 1300 self.pf_axes[n].set_xlabel('')
1292 1301 else:
1293 1302 if self.zlimits is not None:
1294 1303 self.zmin, self.zmax = self.zlimits[n]
1295 1304 ax.collections.remove(ax.collections[0])
1296 1305 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1297 1306 vmin=self.zmin,
1298 1307 vmax=self.zmax,
1299 1308 cmap=self.cmaps[n]
1300 1309 )
1301 1310 if self.showprofile:
1302 1311 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1303 1312 self.pf_axes[n].set_xlabel('')
1304 1313
1305 1314
1306 1315 import datetime
1307 1316 class NoiselessRTILinePlot(Plot):
1308 1317 '''
1309 1318 Plot for RTI data
1310 1319 '''
1311 1320
1312 1321 CODE = 'noiseless_rtiLine'
1313 1322
1314 1323 plot_type = 'scatter'
1315 1324 titles = None
1316 1325 channelList = []
1317 1326 elevationList = []
1318 1327 azimuthList = []
1319 1328 last_noise = None
1320 1329
1321 1330 def setup(self):
1322 1331 self.xaxis = 'Range (Km)'
1323 1332 self.nplots = len(self.data.channels)
1324 1333 self.nrows = int(numpy.ceil(self.nplots/2))
1325 1334 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
1326 1335 self.ylabel = 'Intensity [dB]'
1327 1336 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
1328 1337 self.colorbar = False
1329 1338 self.width = 6
1330 1339 self.height = 4
1331 1340
1332 1341 def update_list(self,dataOut):
1333 1342 if len(self.channelList) == 0:
1334 1343 self.channelList = dataOut.channelList
1335 1344 if len(self.elevationList) == 0:
1336 1345 self.elevationList = dataOut.elevationList
1337 1346 if len(self.azimuthList) == 0:
1338 1347 self.azimuthList = dataOut.azimuthList
1339 1348
1340 1349 def update(self, dataOut):
1341 1350 if len(self.channelList) == 0:
1342 1351 self.update_list(dataOut)
1343 1352
1344 1353 data = {}
1345 1354 meta = {}
1346 1355 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1347 1356 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt)
1348 1357 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1349 1358
1350 1359 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1351 1360 data['noise'] = n0
1352 1361
1353 1362 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1354 1363 noiseless_data = dataOut.getPower() - noise
1355 1364
1356 1365 #print("power, noise:", dataOut.getPower(), n0)
1357 1366 #print(noise)
1358 1367 #print(noiseless_data)
1359 1368
1360 1369 data['noiseless_rtiLine'] = noiseless_data
1361 1370
1362 1371 #print(noiseless_data.shape, self.name)
1363 1372 data['time'] = dataOut.utctime
1364 1373
1365 1374 return data, meta
1366 1375
1367 1376 def plot(self):
1368 1377
1369 1378 self.x = self.data.times
1370 1379 self.y = self.data.yrange
1371 1380 #print(self.data['noiseless_rtiLine'].shape, self.y.shape, self.name)
1372 1381 #ts = self.data['time'][0].squeeze()
1373 1382 if len(self.data['noiseless_rtiLine'])>2 :
1374 1383 self.z = self.data['noiseless_rtiLine'][:, -1,:]
1375 1384 else:
1376 1385 self.z = self.data['noiseless_rtiLine']
1377 1386 #print(self.z.shape, self.y.shape, ts)
1378 1387 #thisDatetime = datetime.datetime.utcfromtimestamp(ts)
1379 1388
1380 1389 for i,ax in enumerate(self.axes):
1381 1390 #self.titles[i] = "Channel {} {}".format(i, thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1382 1391
1383 1392
1384 1393 if ax.firsttime:
1385 1394 #self.xmin = min(self.z)
1386 1395 #self.xmax = max(self.z)
1387 1396 ax.plt_r = ax.plot(self.z[i], self.y)[0]
1388 1397 else:
1389 1398 ax.plt_r.set_data(self.z[i], self.y)
1390 1399
1391 1400
1392 1401
1393 1402 class GeneralProfilePlot(Plot):
1394 1403 '''
1395 1404 Plot for RTI data
1396 1405 '''
1397 1406
1398 1407 CODE = 'general_profilePlot'
1399 1408
1400 1409 plot_type = 'scatter'
1401 1410 titles = None
1402 1411 channelList = []
1403 1412 elevationList = []
1404 1413 azimuthList = []
1405 1414 last_noise = None
1406 1415
1407 1416 def setup(self):
1408 1417 self.xaxis = 'Range (Km)'
1409 1418 self.nplots = len(self.data.channels)
1410 1419 self.nrows = int(numpy.ceil(self.nplots/2))
1411 1420 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
1412 1421 self.ylabel = 'Intensity [dB]'
1413 1422 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
1414 1423 self.colorbar = False
1415 1424 self.width = 6
1416 1425 self.height = 4
1417 1426
1418 1427 def update_list(self,dataOut):
1419 1428 if len(self.channelList) == 0:
1420 1429 self.channelList = dataOut.channelList
1421 1430 if len(self.elevationList) == 0:
1422 1431 self.elevationList = dataOut.elevationList
1423 1432 if len(self.azimuthList) == 0:
1424 1433 self.azimuthList = dataOut.azimuthList
1425 1434
1426 1435 def update(self, dataOut):
1427 1436 if len(self.channelList) == 0:
1428 1437 self.update_list(dataOut)
1429 1438
1430 1439 data = {}
1431 1440 meta = {}
1432 1441
1433 1442 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1434 1443 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1435 1444 data['noise'] = n0
1436 1445
1437 1446 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1438 1447 noiseless_data = dataOut.getPower() - noise
1439 1448
1440 1449 data['noiseless_rtiLine'] = noiseless_data
1441 1450
1442 1451 #print(noiseless_data.shape, self.name)
1443 1452 data['time'] = dataOut.utctime
1444 1453
1445 1454 return data, meta
1446 1455
1447 1456 def plot(self):
1448 1457
1449 1458 self.x = self.data.times
1450 1459 self.y = self.data.yrange
1451 1460 #print(self.data['noiseless_rtiLine'].shape, self.y.shape, self.name)
1452 1461 #ts = self.data['time'][0].squeeze()
1453 1462 if len(self.data['noiseless_rtiLine'])>2 :
1454 1463 self.z = self.data['noiseless_rtiLine'][:, -1,:]
1455 1464 else:
1456 1465 self.z = self.data['noiseless_rtiLine']
1457 1466 #print(self.z.shape, self.y.shape, ts)
1458 1467 #thisDatetime = datetime.datetime.utcfromtimestamp(ts)
1459 1468
1460 1469 for i,ax in enumerate(self.axes):
1461 1470 #self.titles[i] = "Channel {} {}".format(i, thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1462 1471
1463 1472
1464 1473 if ax.firsttime:
1465 1474 #self.xmin = min(self.z)
1466 1475 #self.xmax = max(self.z)
1467 1476 ax.plt_r = ax.plot(self.z[i], self.y)[0]
1468 1477 else:
1469 ax.plt_r.set_data(self.z[i], self.y) No newline at end of file
1478 ax.plt_r.set_data(self.z[i], self.y)
1479
1480
1481 ##########################################################################################################
1482 ########################################## AMISR_V4 ######################################################
1483
1484 class RTIMapPlot(Plot):
1485 '''
1486 Plot for RTI data
1487
1488 Example:
1489
1490 controllerObj = Project()
1491 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1492 ##.......................................................................................
1493 ##.......................................................................................
1494 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1495 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1496 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1497 syncronization=False,shiftChannels=0)
1498
1499 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1500
1501 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1502 opObj01.addParameter(name='code', value=code, format='floatlist')
1503 opObj01.addParameter(name='nCode', value=1, format='int')
1504 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1505 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1506
1507 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1508 opObj12.addParameter(name='minHei', value='90', format='float')
1509 opObj12.addParameter(name='maxHei', value='150', format='float')
1510
1511 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1512 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1513
1514 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1515 opObj11.addParameter(name='n', value='1', format='int')
1516
1517 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1518
1519 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1520 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1521 opObj12.addParameter(name='bField', value='100', format='int')
1522 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1523
1524 '''
1525
1526 CODE = 'rti_skymap'
1527
1528 plot_type = 'scatter'
1529 titles = None
1530 colormap = 'jet'
1531 channelList = []
1532 elevationList = []
1533 azimuthList = []
1534 last_noise = None
1535 flag_setIndex = False
1536 heights = []
1537 dcosx = []
1538 dcosy = []
1539 fullDcosy = None
1540 fullDcosy = None
1541 hindex = []
1542 mapFile = False
1543 ##### BField ####
1544 flagBField = False
1545 dcosxB = []
1546 dcosyB = []
1547 Bmarker = ['+','*','D','x','s','>','o','^']
1548
1549
1550
1551 def setup(self):
1552
1553 self.xaxis = 'Range (Km)'
1554 if len(self.selectedHeightsList) > 0:
1555 self.nplots = len(self.selectedHeightsList)
1556 else:
1557 self.nplots = 4
1558 self.ncols = int(numpy.ceil(self.nplots/2))
1559 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1560 self.ylabel = 'dcosy'
1561 self.xlabel = 'dcosx'
1562 self.colorbar = True
1563 self.width = 6 + 4.1*self.nrows
1564 self.height = 3 + 3.5*self.ncols
1565
1566
1567 if self.extFile!=None:
1568 try:
1569 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1570 full_azi = pointings[:,1]
1571 full_elev = pointings[:,2]
1572 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1573 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1574 mapFile = True
1575 except Exception as e:
1576 self.extFile = None
1577 print(e)
1578
1579
1580
1581
1582 def update_list(self,dataOut):
1583 if len(self.channelList) == 0:
1584 self.channelList = dataOut.channelList
1585 if len(self.elevationList) == 0:
1586 self.elevationList = dataOut.elevationList
1587 if len(self.azimuthList) == 0:
1588 self.azimuthList = dataOut.azimuthList
1589 a = numpy.radians(numpy.asarray(self.azimuthList))
1590 e = numpy.radians(numpy.asarray(self.elevationList))
1591 self.heights = dataOut.heightList
1592 self.dcosx = numpy.cos(e)*numpy.sin(a)
1593 self.dcosy = numpy.cos(e)*numpy.cos(a)
1594
1595 if len(self.bFieldList)>0:
1596 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1597 doy = datetObj.timetuple().tm_yday
1598 year = datetObj.year
1599 # self.dcosxB, self.dcosyB
1600 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1601 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1602
1603 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1604 for ih in range(len(self.bFieldList)):
1605 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1606 for ilon in numpy.arange(nlon):
1607 myx = (alpha[ilon,:,ih])[::-1]
1608 myy = (dcos[ilon,:,ih,0])[::-1]
1609 tck = splrep(myx,myy,s=0)
1610 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1611
1612 myx = (alpha[ilon,:,ih])[::-1]
1613 myy = (dcos[ilon,:,ih,1])[::-1]
1614 tck = splrep(myx,myy,s=0)
1615 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1616 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1617 self.dcosxB.append(alpha_location[:,0,ih])
1618 self.dcosyB.append(alpha_location[:,1,ih])
1619 self.flagBField = True
1620
1621 if len(self.celestialList)>0:
1622 #getBField(self.bFieldList, date)
1623 #pass = kwargs.get('celestial', [])
1624 pass
1625
1626
1627
1628 def update(self, dataOut):
1629
1630 if len(self.channelList) == 0:
1631 self.update_list(dataOut)
1632
1633 if not self.flag_setIndex:
1634 if len(self.selectedHeightsList)>0:
1635 for sel_height in self.selectedHeightsList:
1636 index_list = numpy.where(self.heights >= sel_height)
1637 index_list = index_list[0]
1638 self.hindex.append(index_list[0])
1639 # else:
1640 # k = len(self.heights)
1641 # self.hindex.append(int(k/2))
1642 self.flag_setIndex = True
1643
1644 data = {}
1645 meta = {}
1646
1647 data['rti_skymap'] = dataOut.getPower()
1648 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1649 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1650 data['noise'] = noise
1651
1652 return data, meta
1653
1654 def plot(self):
1655
1656 ######
1657 self.x = self.dcosx
1658 self.y = self.dcosy
1659 self.z = self.data[-1]['rti_skymap']
1660 self.z = numpy.array(self.z, dtype=float)
1661
1662 #print("inde x1 ", self.height_index)
1663 if len(self.hindex) > 0:
1664 index = self.hindex
1665 else:
1666 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1667
1668 #print(index)
1669 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1670 for n, ax in enumerate(self.axes):
1671
1672 if ax.firsttime:
1673
1674
1675 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1676 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1677
1678 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1679 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1680
1681 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1682 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1683
1684
1685 if self.extFile!=None:
1686 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1687 #print(self.fullDcosx)
1688 pass
1689
1690
1691 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1692 s=60, marker="s", vmax = self.zmax)
1693
1694
1695 ax.minorticks_on()
1696 ax.grid(which='major', axis='both')
1697 ax.grid(which='minor', axis='x')
1698
1699 if self.flagBField :
1700
1701 for ih in range(len(self.bFieldList)):
1702 label = str(self.bFieldList[ih]) + ' km'
1703 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1704 label=label, linestyle='--', ms=4.0,lw=0.5)
1705 handles, labels = ax.get_legend_handles_labels()
1706 a = -0.05
1707 b = 1.15 - 1.19*(self.nrows)
1708 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1709
1710 else:
1711
1712 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1713 s=80, marker="s", vmax = self.zmax)
1714
1715 if self.flagBField :
1716 for ih in range(len(self.bFieldList)):
1717 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1718 linestyle='--', ms=4.0,lw=0.5)
1719
1720 # handles, labels = ax.get_legend_handles_labels()
1721 # a = -0.05
1722 # b = 1.15 - 1.19*(self.nrows)
1723 # self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1724
1725
Show file before | Show file after | Hide comments
@@ -1,409 +1,418
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from schainpy.model.graphics.jroplot_base import Plot, plt
11 11
12 12
13 13 class ScopePlot(Plot):
14 14
15 15 '''
16 16 Plot for Scope
17 17 '''
18 18
19 19 CODE = 'scope'
20 20 plot_type = 'scatter'
21 21
22 22 def setup(self):
23 23
24 24 self.xaxis = 'Range (Km)'
25 25 self.nplots = len(self.data.channels)
26 26 self.nrows = int(numpy.ceil(self.nplots/2))
27 27 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
28 28 self.ylabel = 'Intensity [dB]'
29 29 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
30 30 self.colorbar = False
31 31 self.width = 6
32 32 self.height = 4
33 33
34 34 def update(self, dataOut):
35 35
36 36 data = {}
37 37 meta = {
38 38 'nProfiles': dataOut.nProfiles,
39 39 'flagDataAsBlock': dataOut.flagDataAsBlock,
40 40 'profileIndex': dataOut.profileIndex,
41 41 }
42 42 if self.CODE == 'scope':
43 43 data[self.CODE] = dataOut.data
44 44 elif self.CODE == 'pp_power':
45 45 data[self.CODE] = dataOut.dataPP_POWER
46 46 elif self.CODE == 'pp_signal':
47 47 data[self.CODE] = dataOut.dataPP_POW
48 48 elif self.CODE == 'pp_velocity':
49 49 data[self.CODE] = dataOut.dataPP_DOP
50 50 elif self.CODE == 'pp_specwidth':
51 51 data[self.CODE] = dataOut.dataPP_WIDTH
52 52
53 53 return data, meta
54 54
55 55 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
56 56
57 57 yreal = y[channelIndexList,:].real
58 58 yimag = y[channelIndexList,:].imag
59 59 Maintitle = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
60 60 self.xlabel = "Range (Km)"
61 61 self.ylabel = "Intensity - IQ"
62 62
63 63 self.y = yreal
64 64 self.x = x
65 65
66 66 for i,ax in enumerate(self.axes):
67 67 title = "Channel %d" %(i)
68 68 if ax.firsttime:
69 69 self.xmin = min(x)
70 70 self.xmax = max(x)
71 71 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
72 72 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
73 73 else:
74 74 ax.plt_r.set_data(x, yreal[i,:])
75 75 ax.plt_i.set_data(x, yimag[i,:])
76 76 plt.suptitle(Maintitle)
77 77
78 78 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
79 79 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
80 80 yreal = y.real
81 81 yreal = 10*numpy.log10(yreal)
82 82 self.y = yreal
83 83 title = wintitle + " Power: %s" %(thisDatetime.strftime("%d-%b-%Y"))
84 84 self.xlabel = "Range (Km)"
85 85 self.ylabel = "Intensity [dB]"
86 86
87 87
88 88 self.titles[0] = title
89 89
90 90 for i,ax in enumerate(self.axes):
91 91 title = "Channel %d" %(i)
92 92 ychannel = yreal[i,:]
93 93
94 94 if ax.firsttime:
95 95 self.xmin = min(x)
96 96 self.xmax = max(x)
97 97 ax.plt_r = ax.plot(x, ychannel)[0]
98 98 else:
99 99 ax.plt_r.set_data(x, ychannel)
100 100
101 101
102 102 def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):
103 103
104 104
105 105 y = y[channelIndexList,:]
106 106 yreal = y.real
107 107 yreal = 10*numpy.log10(yreal)
108 108 self.y = yreal
109 109 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
110 110 self.xlabel = "Range (Km)"
111 111 self.ylabel = "Intensity"
112 112 self.xmin = min(x)
113 113 self.xmax = max(x)
114 114
115 115 self.titles[0] =title
116 116 for i,ax in enumerate(self.axes):
117 117 title = "Channel %d" %(i)
118 118
119 119 ychannel = yreal[i,:]
120 120
121 121 if ax.firsttime:
122 122 ax.plt_r = ax.plot(x, ychannel)[0]
123 123 else:
124 124 #pass
125 125 ax.plt_r.set_data(x, ychannel)
126 126
127 127 def plot_weathervelocity(self, x, y, channelIndexList, thisDatetime, wintitle):
128 128
129 129 x = x[channelIndexList,:]
130 130 yreal = y
131 131 self.y = yreal
132 132 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
133 133 self.xlabel = "Velocity (m/s)"
134 134 self.ylabel = "Range (Km)"
135 135 self.xmin = numpy.min(x)
136 136 self.xmax = numpy.max(x)
137 137 self.titles[0] =title
138 138 for i,ax in enumerate(self.axes):
139 139 title = "Channel %d" %(i)
140 140 xchannel = x[i,:]
141 141 if ax.firsttime:
142 142 ax.plt_r = ax.plot(xchannel, yreal)[0]
143 143 else:
144 144 #pass
145 145 ax.plt_r.set_data(xchannel, yreal)
146 146
147 147 def plot_weatherspecwidth(self, x, y, channelIndexList, thisDatetime, wintitle):
148 148
149 149 x = x[channelIndexList,:]
150 150 yreal = y
151 151 self.y = yreal
152 152 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
153 153 self.xlabel = "width "
154 154 self.ylabel = "Range (Km)"
155 155 self.xmin = numpy.min(x)
156 156 self.xmax = numpy.max(x)
157 157 self.titles[0] =title
158 158 for i,ax in enumerate(self.axes):
159 159 title = "Channel %d" %(i)
160 160 xchannel = x[i,:]
161 161 if ax.firsttime:
162 162 ax.plt_r = ax.plot(xchannel, yreal)[0]
163 163 else:
164 164 #pass
165 165 ax.plt_r.set_data(xchannel, yreal)
166 166
167 167 def plot(self):
168 168 if self.channels:
169 169 channels = self.channels
170 170 else:
171 171 channels = self.data.channels
172 172
173 173 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
174 174
175 175 scope = self.data[-1][self.CODE]
176 176
177 177 if self.data.flagDataAsBlock:
178 178
179 179 for i in range(self.data.nProfiles):
180 180
181 181 wintitle1 = " [Profile = %d] " %i
182 182 if self.CODE =="scope":
183 183 if self.type == "power":
184 184 self.plot_power(self.data.yrange,
185 185 scope[:,i,:],
186 186 channels,
187 187 thisDatetime,
188 188 wintitle1
189 189 )
190 190
191 191 if self.type == "iq":
192 192 self.plot_iq(self.data.yrange,
193 193 scope[:,i,:],
194 194 channels,
195 195 thisDatetime,
196 196 wintitle1
197 197 )
198 198 if self.CODE=="pp_power":
199 199 self.plot_weatherpower(self.data.yrange,
200 200 scope[:,i,:],
201 201 channels,
202 202 thisDatetime,
203 203 wintitle
204 204 )
205 205 if self.CODE=="pp_signal":
206 206 self.plot_weatherpower(self.data.yrange,
207 207 scope[:,i,:],
208 208 channels,
209 209 thisDatetime,
210 210 wintitle
211 211 )
212 212 if self.CODE=="pp_velocity":
213 213 self.plot_weathervelocity(scope[:,i,:],
214 214 self.data.yrange,
215 215 channels,
216 216 thisDatetime,
217 217 wintitle
218 218 )
219 219 if self.CODE=="pp_spcwidth":
220 220 self.plot_weatherspecwidth(scope[:,i,:],
221 221 self.data.yrange,
222 222 channels,
223 223 thisDatetime,
224 224 wintitle
225 225 )
226 226 else:
227 227 wintitle = " [Profile = %d] " %self.data.profileIndex
228 228 if self.CODE== "scope":
229 229 if self.type == "power":
230 230 self.plot_power(self.data.yrange,
231 231 scope,
232 232 channels,
233 233 thisDatetime,
234 234 wintitle
235 235 )
236 236
237 237 if self.type == "iq":
238 238 self.plot_iq(self.data.yrange,
239 239 scope,
240 240 channels,
241 241 thisDatetime,
242 242 wintitle
243 243 )
244 244 if self.CODE=="pp_power":
245 245 self.plot_weatherpower(self.data.yrange,
246 246 scope,
247 247 channels,
248 248 thisDatetime,
249 249 wintitle
250 250 )
251 251 if self.CODE=="pp_signal":
252 252 self.plot_weatherpower(self.data.yrange,
253 253 scope,
254 254 channels,
255 255 thisDatetime,
256 256 wintitle
257 257 )
258 258 if self.CODE=="pp_velocity":
259 259 self.plot_weathervelocity(scope,
260 260 self.data.yrange,
261 261 channels,
262 262 thisDatetime,
263 263 wintitle
264 264 )
265 265 if self.CODE=="pp_specwidth":
266 266 self.plot_weatherspecwidth(scope,
267 267 self.data.yrange,
268 268 channels,
269 269 thisDatetime,
270 270 wintitle
271 271 )
272 272
273 273
274 274 class PulsepairPowerPlot(ScopePlot):
275 275 '''
276 276 Plot for P= S+N
277 277 '''
278 278
279 279 CODE = 'pp_power'
280 280 plot_type = 'scatter'
281 281
282 282 class PulsepairVelocityPlot(ScopePlot):
283 283 '''
284 284 Plot for VELOCITY
285 285 '''
286 286 CODE = 'pp_velocity'
287 287 plot_type = 'scatter'
288 288
289 289 class PulsepairSpecwidthPlot(ScopePlot):
290 290 '''
291 291 Plot for WIDTH
292 292 '''
293 293 CODE = 'pp_specwidth'
294 294 plot_type = 'scatter'
295 295
296 296 class PulsepairSignalPlot(ScopePlot):
297 297 '''
298 298 Plot for S
299 299 '''
300 300
301 301 CODE = 'pp_signal'
302 302 plot_type = 'scatter'
303 303
304 304
305 305 class Spectra2DPlot(Plot):
306 306 '''
307 307 Plot for 2D Spectra data
308 Necessary data as Block
309 you could use profiles2Block Operation
310
311 Example:
312 # opObj11 = volts_proc.addOperation(name='profiles2Block', optype='other')
313 # # opObj11.addParameter(name='n', value=10, format='int')
314 # opObj11.addParameter(name='timeInterval', value='2', format='int')
315
316 # opObj12 = volts_proc.addOperation(name='Spectra2DPlot', optype='external')
308 317 '''
309 318
310 319 CODE = 'spc'
311 320 colormap = 'jet'
312 321 plot_type = 'pcolor'
313 322 buffering = False
314 323 channelList = []
315 324 elevationList = []
316 325 azimuthList = []
317 326
318 327 def setup(self):
319 328
320 329 self.nplots = len(self.data.channels)
321 330 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
322 331 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
323 332 self.height = 3.4 * self.nrows
324 333
325 334 self.cb_label = 'dB'
326 335 if self.showprofile:
327 336 self.width = 5.2 * self.ncols
328 337 else:
329 338 self.width = 4.2* self.ncols
330 339 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
331 340 self.ylabel = 'Range [km]'
332 341
333 342
334 343 def update_list(self,dataOut):
335 344 if len(self.channelList) == 0:
336 345 self.channelList = dataOut.channelList
337 346 if len(self.elevationList) == 0:
338 347 self.elevationList = dataOut.elevationList
339 348 if len(self.azimuthList) == 0:
340 349 self.azimuthList = dataOut.azimuthList
341 350
342 351 def update(self, dataOut):
343 352
344 353 self.update_list(dataOut)
345 354 data = {}
346 355 meta = {}
347 356
348 357
349 358 spectrum = numpy.fft.fftshift(numpy.fft.fft2(dataOut.data, axes=(1,2)))
350 359 z = numpy.abs(spectrum)
351 360 phase = numpy.angle(spectrum)
352 361
353 362 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
354 363 spc = 10*numpy.log10(z)
355 364 dt1 = dataOut.ippSeconds
356 365 dt2 = dataOut.radarControllerHeaderObj.heightResolution/150000
357 366 data['spc'] = spc
358 367 #print(spc[0].shape)
359 368 data['phase'] = phase
360 369 f1 = numpy.fft.fftshift(numpy.fft.fftfreq(spectrum.shape[1],d=dt1)/1000)
361 370 f2 = numpy.fft.fftshift(numpy.fft.fftfreq(spectrum.shape[2],d=dt2)/1000)
362 371 meta['range'] = (f1, f2)
363 372
364 373 return data, meta
365 374
366 375 def plot(self):
367 376 x = self.data.range[0]
368 377 y = self.data.range[1]
369 378 self.xlabel = "Frequency (kHz)"
370 379 self.ylabel = "Frequency (kHz)"
371 380
372 381 # if self.xaxis == "frequency":
373 382 # x = self.data.range[1]
374 383 # y = self.data.range[2]
375 384 # self.xlabel = "Frequency (kHz)"
376 385 # self.ylabel = "Frequency (kHz)"
377 386 # else:
378 387 # x = self.data.xrange[2]
379 388 # self.xlabel = "Velocity (m/s)"
380 389
381 390
382 391 self.titles = []
383 392 self.y = y
384 393
385 394 data = self.data[-1]
386 395 z = data['spc']
387 396 #print(z.shape, x.shape, y.shape)
388 397 for n, ax in enumerate(self.axes):
389 398 # noise = self.data['noise'][n][0]
390 399 #print(noise)
391 400
392 401 if ax.firsttime:
393 402 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
394 403 self.xmin = self.xmin if self.xmin else -self.xmax
395 404 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
396 405 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
397 406 ax.plt = ax.pcolormesh(x, y, z[n].T,
398 407 vmin=self.zmin,
399 408 vmax=self.zmax,
400 409 cmap=plt.get_cmap(self.colormap)
401 410 )
402 411
403 412 else:
404 413 ax.plt.set_array(z[n].T.ravel())
405 414
406 415 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
407 416 self.titles.append('CH {}: {:2.1f}elv {:2.1f} az '.format(self.channelList[n], self.elevationList[n], self.azimuthList[n]))
408 417 else:
409 418 self.titles.append('CH {}: '.format(self.channelList[n])) No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,755 +1,783
1 1 ''''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5
6 6 @upgrade: 2021, Joab Apaza
7 7
8 8 '''
9 9
10 10 import os
11 11 import sys
12 12 import glob
13 13 import fnmatch
14 14 import datetime
15 15 import time
16 16 import re
17 17 import h5py
18 18 import numpy
19 19
20 20 try:
21 21 from gevent import sleep
22 22 except:
23 23 from time import sleep
24 24
25 25 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader,ProcessingHeader
26 26 from schainpy.model.data.jrodata import Voltage
27 27 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
28 28 from numpy import imag
29 29 from schainpy.utils import log
30 30
31 31
32 32 class AMISRReader(ProcessingUnit):
33 33 '''
34 34 classdocs
35 35 '''
36 36
37 37 def __init__(self):
38 38 '''
39 39 Constructor
40 40 '''
41 41
42 42 ProcessingUnit.__init__(self)
43 43
44 44 self.set = None
45 45 self.subset = None
46 46 self.extension_file = '.h5'
47 47 self.dtc_str = 'dtc'
48 48 self.dtc_id = 0
49 49 self.status = True
50 50 self.isConfig = False
51 51 self.dirnameList = []
52 52 self.filenameList = []
53 53 self.fileIndex = None
54 54 self.flagNoMoreFiles = False
55 55 self.flagIsNewFile = 0
56 56 self.filename = ''
57 57 self.amisrFilePointer = None
58 58
59 59 self.beamCodeMap = None
60 60 self.azimuthList = []
61 61 self.elevationList = []
62 62 self.dataShape = None
63 63 self.flag_old_beams = False
64
65
64
65 self.flagAsync = False #Use when the experiment has no syncronization
66 self.shiftChannels = 0
66 67 self.profileIndex = 0
67 68
68 69
69 70 self.beamCodeByFrame = None
70 71 self.radacTimeByFrame = None
71 72
72 73 self.dataset = None
73 74
74 75 self.__firstFile = True
75 76
76 77 self.buffer = None
77 78
78 79 self.timezone = 'ut'
79 80
80 81 self.__waitForNewFile = 20
81 82 self.__filename_online = None
82 83 #Is really necessary create the output object in the initializer
83 84 self.dataOut = Voltage()
84 85 self.dataOut.error=False
85 86 self.margin_days = 1
86 87 self.flag_ignoreFiles = False #to activate the ignoring Files flag
87 88 self.flag_standby = False # just keep waiting, use when ignoring files
88 89 self.ignStartDateTime=None
89 90 self.ignEndDateTime=None
90 91
91 92 def setup(self,path=None,
92 93 startDate=None,
93 94 endDate=None,
94 95 startTime=None,
95 96 endTime=None,
96 97 walk=True,
97 98 timezone='ut',
98 99 all=0,
99 100 code = 1,
100 101 nCode = 1,
101 102 nBaud = 0,
102 103 nOsamp = 0,
103 104 online=False,
104 105 old_beams=False,
105 106 margin_days=1,
106 107 nFFT = None,
107 108 nChannels = None,
108 109 ignStartDate=None,
109 110 ignEndDate=None,
110 111 ignStartTime=None,
111 112 ignEndTime=None,
113 syncronization=True,
114 shiftChannels=0
112 115 ):
113 116
114 117
115 118
116 119 self.timezone = timezone
117 120 self.all = all
118 121 self.online = online
119 122 self.flag_old_beams = old_beams
120 123 self.code = code
121 124 self.nCode = int(nCode)
122 125 self.nBaud = int(nBaud)
123 126 self.nOsamp = int(nOsamp)
124 127 self.margin_days = margin_days
125 128 self.__sampleRate = None
126
129 self.flagAsync = not syncronization
130 self.shiftChannels = shiftChannels
127 131 self.nFFT = nFFT
128 132 self.nChannels = nChannels
129 133 if ignStartTime!=None and ignEndTime!=None:
130 134 if ignStartDate!=None and ignEndDate!=None:
131 135 self.ignStartDateTime=datetime.datetime.combine(ignStartDate,ignStartTime)
132 136 self.ignEndDateTime=datetime.datetime.combine(ignEndDate,ignEndTime)
133 137 else:
134 138 self.ignStartDateTime=datetime.datetime.combine(startDate,ignStartTime)
135 139 self.ignEndDateTime=datetime.datetime.combine(endDate,ignEndTime)
136 140 self.flag_ignoreFiles = True
137 141
138 142 #self.findFiles()
139 143 if not(online):
140 144 #Busqueda de archivos offline
141 145 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk,)
142 146 else:
143 147 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
144 148
145 149 if not(self.filenameList):
146 150 raise schainpy.admin.SchainWarning("There is no files into the folder: %s"%(path))
147 151 #sys.exit(0)
148 152 self.dataOut.error = True
149 153
150 154 self.fileIndex = 0
151 155
152 156 self.readNextFile(online)
153 157
154 158 '''
155 159 Add code
156 160 '''
157 161 self.isConfig = True
158 162 # print("Setup Done")
159 163 pass
160 164
161 165
162 166 def readAMISRHeader(self,fp):
163 167
164 168 if self.isConfig and (not self.flagNoMoreFiles):
165 169 newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
166 170 if self.dataShape != newShape and newShape != None and not self.flag_standby:
167 171 raise schainpy.admin.SchainError("NEW FILE HAS A DIFFERENT SHAPE: ")
168 172 print(self.dataShape,newShape,"\n")
169 173 return 0
170 174 else:
171 175 self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
172 176
173 177
174 178 header = 'Raw11/Data/RadacHeader'
175 179 if self.nChannels == None:
176 180 expFile = fp['Setup/Experimentfile'][()].decode()
177 181 linesExp = expFile.split("\n")
178 182 a = [line for line in linesExp if "nbeamcodes" in line]
179 183 self.nChannels = int(a[0][11:])
180 184
185 if not self.flagAsync: #for experiments with no syncronization
186 self.shiftChannels = 0
187
188
189
181 190 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
191
192
182 193 if (self.startDate > datetime.date(2021, 7, 15)) or self.flag_old_beams: #Se cambió la forma de extracción de Apuntes el 17 o forzar con flag de reorganización
183 194 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
184 195 self.trueBeams = self.beamcodeFile.split("\n")
185 196 self.trueBeams.pop()#remove last
186 197 if self.nFFT == None:
187 198 log.error("FFT or number of repetitions per channels is needed",self.name)
188 199 beams_idx = [k*self.nFFT for k in range(self.nChannels)]
189 200 beams = [self.trueBeams[b] for b in beams_idx]
190 201 self.beamCode = [int(x, 16) for x in beams]
191 202
203 if(self.flagAsync and self.shiftChannels == 0):
204 initBeam = self.beamCodeByPulse[0, 0]
205 self.shiftChannels = numpy.argwhere(self.beamCode ==initBeam)[0,0]
206
192 207 else:
193 208 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
194 209 self.beamCode = _beamCode[0,:]
195 210
211
212
213
196 214 if self.beamCodeMap == None:
197 215 self.beamCodeMap = fp['Setup/BeamcodeMap']
198 216 for beam in self.beamCode:
199 217 beamAziElev = numpy.where(self.beamCodeMap[:,0]==beam)
200 218 beamAziElev = beamAziElev[0].squeeze()
201 219 self.azimuthList.append(self.beamCodeMap[beamAziElev,1])
202 220 self.elevationList.append(self.beamCodeMap[beamAziElev,2])
203 221 #print("Beamssss: ",self.beamCodeMap[beamAziElev,1],self.beamCodeMap[beamAziElev,2])
204 222 #print(self.beamCode)
205 223 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
206 224 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
207 225 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
208 226 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
209 227 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
210 228 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
211 229 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
212 230 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
213 231 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
214 232 self.frequency = fp.get('Rx/Frequency')
215 233 txAus = fp.get('Raw11/Data/Pulsewidth') #seconds
216 234 self.baud = fp.get('Raw11/Data/TxBaud')
217 235 sampleRate = fp.get('Rx/SampleRate')
218 236 self.__sampleRate = sampleRate[()]
219 237 self.nblocks = self.pulseCount.shape[0] #nblocks
220 238 self.profPerBlockRAW = self.pulseCount.shape[1] #profiles per block in raw data
221 239 self.nprofiles = self.pulseCount.shape[1] #nprofile
222 self.nsa = self.nsamplesPulse[0,0] #ngates
240 #self.nsa = self.nsamplesPulse[0,0] #ngates
241 self.nsa = len(self.rangeFromFile[0])
223 242 self.nchannels = len(self.beamCode)
224 243 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
225 244 #print("IPPS secs: ",self.ippSeconds)
226 245 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
227 246 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
228 247
229 248 #filling radar controller header parameters
230 249 self.__ippKm = self.ippSeconds *.15*1e6 # in km
231 250 #self.__txA = txAus[()]*.15 #(ipp[us]*.15km/1us) in km
232 251 self.__txA = txAus[()] #seconds
233 252 self.__txAKm = self.__txA*1e6*.15
234 253 self.__txB = 0
235 254 nWindows=1
236 255 self.__nSamples = self.nsa
237 256 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
238 257 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
239 258 #print("amisr-ipp:",self.ippSeconds, self.__ippKm)
240 259 #for now until understand why the code saved is different (code included even though code not in tuf file)
241 260 #self.__codeType = 0
242 261 # self.__nCode = None
243 262 # self.__nBaud = None
244 263 self.__code = self.code
245 264 self.__codeType = 0
246 265 if self.code != None:
247 266 self.__codeType = 1
248 267 self.__nCode = self.nCode
249 268 self.__nBaud = self.nBaud
250 269 self.__baudTX = self.__txA/(self.nBaud)
251 270 #self.__code = 0
252 271
253 272 #filling system header parameters
254 273 self.__nSamples = self.nsa
255 274 self.newProfiles = self.nprofiles/self.nchannels
256 275 self.__channelList = [n for n in range(self.nchannels)]
257 276
258 277 self.__frequency = self.frequency[0][0]
259 278
260 279
261 280 return 1
262 281
263 282
264 283 def createBuffers(self):
265 284
266 285 pass
267 286
268 287 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
269 288 self.path = path
270 289 self.startDate = startDate
271 290 self.endDate = endDate
272 291 self.startTime = startTime
273 292 self.endTime = endTime
274 293 self.walk = walk
275 294
276 295
277 296 def __checkPath(self):
278 297 if os.path.exists(self.path):
279 298 self.status = 1
280 299 else:
281 300 self.status = 0
282 301 print('Path:%s does not exists'%self.path)
283 302
284 303 return
285 304
286 305
287 306 def __selDates(self, amisr_dirname_format):
288 307 try:
289 308 year = int(amisr_dirname_format[0:4])
290 309 month = int(amisr_dirname_format[4:6])
291 310 dom = int(amisr_dirname_format[6:8])
292 311 thisDate = datetime.date(year,month,dom)
293 312 #margen de un día extra, igual luego se filtra for fecha y hora
294 313 if (thisDate>=(self.startDate - datetime.timedelta(days=self.margin_days)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)):
295 314 return amisr_dirname_format
296 315 except:
297 316 return None
298 317
299 318
300 319 def __findDataForDates(self,online=False):
301 320
302 321 if not(self.status):
303 322 return None
304 323
305 324 pat = '\d+.\d+'
306 325 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
307 326 dirnameList = [x for x in dirnameList if x!=None]
308 327 dirnameList = [x.string for x in dirnameList]
309 328 if not(online):
310 329 dirnameList = [self.__selDates(x) for x in dirnameList]
311 330 dirnameList = [x for x in dirnameList if x!=None]
312 331 if len(dirnameList)>0:
313 332 self.status = 1
314 333 self.dirnameList = dirnameList
315 334 self.dirnameList.sort()
316 335 else:
317 336 self.status = 0
318 337 return None
319 338
320 339 def __getTimeFromData(self):
321 340 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
322 341 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
323 342
324 343 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
325 344 print('........................................')
326 345 filter_filenameList = []
327 346 self.filenameList.sort()
328 347 total_files = len(self.filenameList)
329 348 #for i in range(len(self.filenameList)-1):
330 349 for i in range(total_files):
331 350 filename = self.filenameList[i]
332 351 #print("file-> ",filename)
333 352 try:
334 353 fp = h5py.File(filename,'r')
335 354 time_str = fp.get('Time/RadacTimeString')
336 355
337 356 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
338 357 #startDateTimeStr_File = "2019-12-16 09:21:11"
339 358 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
340 359 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
341 360
342 361 #endDateTimeStr_File = "2019-12-16 11:10:11"
343 362 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
344 363 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
345 364 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
346 365
347 366 fp.close()
348 367
349 368 #print("check time", startDateTime_File)
350 369 if self.timezone == 'lt':
351 370 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
352 371 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
353 372 if (startDateTime_File >=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
354 373 filter_filenameList.append(filename)
355 374
356 375 if (startDateTime_File>endDateTime_Reader):
357 376 break
358 377 except Exception as e:
359 378 log.warning("Error opening file {} -> {}".format(os.path.split(filename)[1],e))
360 379
361 380 filter_filenameList.sort()
362 381 self.filenameList = filter_filenameList
363 382
364 383 return 1
365 384
366 385 def __filterByGlob1(self, dirName):
367 386 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
368 387 filter_files.sort()
369 388 filterDict = {}
370 389 filterDict.setdefault(dirName)
371 390 filterDict[dirName] = filter_files
372 391 return filterDict
373 392
374 393 def __getFilenameList(self, fileListInKeys, dirList):
375 394 for value in fileListInKeys:
376 395 dirName = list(value.keys())[0]
377 396 for file in value[dirName]:
378 397 filename = os.path.join(dirName, file)
379 398 self.filenameList.append(filename)
380 399
381 400
382 401 def __selectDataForTimes(self, online=False):
383 402 #aun no esta implementado el filtro for tiempo-> implementado en readNextFile
384 403 if not(self.status):
385 404 return None
386 405
387 406 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
388 407 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
389 408 self.__getFilenameList(fileListInKeys, dirList)
390 409 if not(online):
391 410 #filtro por tiempo
392 411 if not(self.all):
393 412 self.__getTimeFromData()
394 413
395 414 if len(self.filenameList)>0:
396 415 self.status = 1
397 416 self.filenameList.sort()
398 417 else:
399 418 self.status = 0
400 419 return None
401 420
402 421 else:
403 422 #get the last file - 1
404 423 self.filenameList = [self.filenameList[-2]]
405 424 new_dirnameList = []
406 425 for dirname in self.dirnameList:
407 426 junk = numpy.array([dirname in x for x in self.filenameList])
408 427 junk_sum = junk.sum()
409 428 if junk_sum > 0:
410 429 new_dirnameList.append(dirname)
411 430 self.dirnameList = new_dirnameList
412 431 return 1
413 432
414 433 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
415 434 endTime=datetime.time(23,59,59),walk=True):
416 435
417 436 if endDate ==None:
418 437 startDate = datetime.datetime.utcnow().date()
419 438 endDate = datetime.datetime.utcnow().date()
420 439
421 440 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
422 441
423 442 self.__checkPath()
424 443
425 444 self.__findDataForDates(online=True)
426 445
427 446 self.dirnameList = [self.dirnameList[-1]]
428 447
429 448 self.__selectDataForTimes(online=True)
430 449
431 450 return
432 451
433 452
434 453 def searchFilesOffLine(self,
435 454 path,
436 455 startDate,
437 456 endDate,
438 457 startTime=datetime.time(0,0,0),
439 458 endTime=datetime.time(23,59,59),
440 459 walk=True):
441 460
442 461 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
443 462
444 463 self.__checkPath()
445 464
446 465 self.__findDataForDates()
447 466
448 467 self.__selectDataForTimes()
449 468
450 469 for i in range(len(self.filenameList)):
451 470 print("%s" %(self.filenameList[i]))
452 471
453 472 return
454 473
455 474 def __setNextFileOffline(self):
456 475
457 476 try:
458 477 self.filename = self.filenameList[self.fileIndex]
459 478 self.amisrFilePointer = h5py.File(self.filename,'r')
460 479 self.fileIndex += 1
461 480 except:
462 481 self.flagNoMoreFiles = 1
463 482 raise schainpy.admin.SchainError('No more files to read')
464 483 return 0
465 484
466 485 self.flagIsNewFile = 1
467 486 print("Setting the file: %s"%self.filename)
468 487
469 488 return 1
470 489
471 490
472 491 def __setNextFileOnline(self):
473 492 filename = self.filenameList[0]
474 493 if self.__filename_online != None:
475 494 self.__selectDataForTimes(online=True)
476 495 filename = self.filenameList[0]
477 496 wait = 0
478 497 self.__waitForNewFile=300 ## DEBUG:
479 498 while self.__filename_online == filename:
480 499 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
481 500 if wait == 5:
482 501 self.flagNoMoreFiles = 1
483 502 return 0
484 503 sleep(self.__waitForNewFile)
485 504 self.__selectDataForTimes(online=True)
486 505 filename = self.filenameList[0]
487 506 wait += 1
488 507
489 508 self.__filename_online = filename
490 509
491 510 self.amisrFilePointer = h5py.File(filename,'r')
492 511 self.flagIsNewFile = 1
493 512 self.filename = filename
494 513 print("Setting the file: %s"%self.filename)
495 514 return 1
496 515
497 516
498 517 def readData(self):
499 518 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
500 519 re = buffer[:,:,:,0]
501 520 im = buffer[:,:,:,1]
502 521 dataset = re + im*1j
503 522
504 523 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
505 524 timeset = self.radacTime[:,0]
506 525
507 526 return dataset,timeset
508 527
509 528 def reshapeData(self):
510 529 #print(self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa)
511 530 channels = self.beamCodeByPulse[0,:]
512 531 nchan = self.nchannels
513 532 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
514 533 nblocks = self.nblocks
515 534 nsamples = self.nsa
516 535 #print("Channels: ",self.nChannels)
517 536 #Dimensions : nChannels, nProfiles, nSamples
518 537 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
519 538 ############################################
520 539 profPerCH = int(self.profPerBlockRAW / (self.nFFT* self.nChannels))
521 540 #profPerCH = int(self.profPerBlockRAW / self.nChannels)
522 541 for thisChannel in range(nchan):
523 542
543 ich = thisChannel
524 544
525 idx_ch = [self.nFFT*(thisChannel + nchan*k) for k in range(profPerCH)]
545 idx_ch = [self.nFFT*(ich + nchan*k) for k in range(profPerCH)]
526 546 #print(idx_ch)
527 547 if self.nFFT > 1:
528 548 aux = [numpy.arange(i, i+self.nFFT) for i in idx_ch]
529 549 idx_ch = None
530 550 idx_ch =aux
531 551 idx_ch = numpy.array(idx_ch, dtype=int).flatten()
532 552 else:
533 553 idx_ch = numpy.array(idx_ch, dtype=int)
534 554
535 #print(thisChannel,profPerCH,idx_ch)
536 #print(numpy.where(channels==self.beamCode[thisChannel])[0])
537 #new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[thisChannel])[0],:]
538 new_block[:,thisChannel,:,:] = self.dataset[:,idx_ch,:]
555 #print(ich,profPerCH,idx_ch)
556 #print(numpy.where(channels==self.beamCode[ich])[0])
557 #new_block[:,ich,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[ich])[0],:]
558 new_block[:,ich,:,:] = self.dataset[:,idx_ch,:]
539 559
540 560 new_block = numpy.transpose(new_block, (1,0,2,3))
541 561 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
542
562 if self.flagAsync:
563 new_block = numpy.roll(new_block, self.shiftChannels, axis=0)
543 564 return new_block
544 565
545 566 def updateIndexes(self):
546 567
547 568 pass
548 569
549 570 def fillJROHeader(self):
550 571
551 572 #fill radar controller header
552 573
553 574 #fill system header
554 575 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
555 576 nProfiles=self.newProfiles,
556 577 nChannels=len(self.__channelList),
557 578 adcResolution=14,
558 579 pciDioBusWidth=32)
559 580
560 581 self.dataOut.type = "Voltage"
561 582 self.dataOut.data = None
562 583 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
563 584 # self.dataOut.nChannels = 0
564 585
565 586 # self.dataOut.nHeights = 0
566 587
567 588 self.dataOut.nProfiles = self.newProfiles*self.nblocks
568 589 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
569 590 ranges = numpy.reshape(self.rangeFromFile[()],(-1))
570 591 self.dataOut.heightList = ranges/1000.0 #km
571 592 self.dataOut.channelList = self.__channelList
572 593
573 594 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
574 595
575 596 # self.dataOut.channelIndexList = None
576 597
577 598
578 self.dataOut.azimuthList = numpy.array(self.azimuthList)
579 self.dataOut.elevationList = numpy.array(self.elevationList)
580 self.dataOut.codeList = numpy.array(self.beamCode)
581
599 #self.dataOut.azimuthList = numpy.roll( numpy.array(self.azimuthList) ,self.shiftChannels)
600 #self.dataOut.elevationList = numpy.roll(numpy.array(self.elevationList) ,self.shiftChannels)
601 #self.dataOut.codeList = numpy.roll(numpy.array(self.beamCode), self.shiftChannels)
602
603 self.dataOut.azimuthList = self.azimuthList
604 self.dataOut.elevationList = self.elevationList
605 self.dataOut.codeList = self.beamCode
582 606
583 607
584 608
585 609 #print(self.dataOut.elevationList)
586 610 self.dataOut.flagNoData = True
587 611
588 612 #Set to TRUE if the data is discontinuous
589 613 self.dataOut.flagDiscontinuousBlock = False
590 614
591 615 self.dataOut.utctime = None
592 616
593 617 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
594 618 if self.timezone == 'lt':
595 619 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
596 620 else:
597 621 self.dataOut.timeZone = 0 #by default time is UTC
598 622
599 623 self.dataOut.dstFlag = 0
600 624 self.dataOut.errorCount = 0
601 625 self.dataOut.nCohInt = 1
602 626 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
603 627 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
604 628 self.dataOut.flagShiftFFT = False
605 629 self.dataOut.ippSeconds = self.ippSeconds
606 630 self.dataOut.ipp = self.__ippKm
607 631 self.dataOut.nCode = self.__nCode
608 632 self.dataOut.code = self.__code
609 633 self.dataOut.nBaud = self.__nBaud
610 634
611 635
612 636 self.dataOut.frequency = self.__frequency
613 637 self.dataOut.realtime = self.online
614 638
615 639 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
616 640 txA=self.__txAKm,
617 641 txB=0,
618 642 nWindows=1,
619 643 nHeights=self.__nSamples,
620 644 firstHeight=self.__firstHeight,
621 645 codeType=self.__codeType,
622 646 nCode=self.__nCode, nBaud=self.__nBaud,
623 647 code = self.__code,
624 648 nOsamp=self.nOsamp,
625 649 frequency = self.__frequency,
626 650 sampleRate= self.__sampleRate,
627 651 fClock=self.__sampleRate)
628 652
629 653
630 654 self.dataOut.radarControllerHeaderObj.heightList = ranges/1000.0 #km
631 655 self.dataOut.radarControllerHeaderObj.heightResolution = self.__deltaHeight
632 656 self.dataOut.radarControllerHeaderObj.rangeIpp = self.__ippKm #km
633 657 self.dataOut.radarControllerHeaderObj.rangeTxA = self.__txA*1e6*.15 #km
634 658 self.dataOut.radarControllerHeaderObj.nChannels = self.nchannels
635 659 self.dataOut.radarControllerHeaderObj.channelList = self.__channelList
636 660 self.dataOut.radarControllerHeaderObj.azimuthList = self.azimuthList
637 661 self.dataOut.radarControllerHeaderObj.elevationList = self.elevationList
638 662 self.dataOut.radarControllerHeaderObj.dtype = "Voltage"
639 663 self.dataOut.ippSeconds = self.ippSeconds
664 <<<<<<< HEAD
640 665 self.dataOut.ippFactor = self.nchannels*self.nFFT
666 =======
667 self.dataOut.ippFactor = self.nFFT
668 >>>>>>> 37cccf17c7b80521b59b978cb30e4ab2e6f37fce
641 669 pass
642 670
643 671 def readNextFile(self,online=False):
644 672
645 673 if not(online):
646 674 newFile = self.__setNextFileOffline()
647 675 else:
648 676 newFile = self.__setNextFileOnline()
649 677
650 678 if not(newFile):
651 679 self.dataOut.error = True
652 680 return 0
653 681
654 682 if not self.readAMISRHeader(self.amisrFilePointer):
655 683 self.dataOut.error = True
656 684 return 0
657 685
658 686 #self.createBuffers()
659 687 self.fillJROHeader()
660 688
661 689 #self.__firstFile = False
662 690
663 691 self.dataset,self.timeset = self.readData()
664 692
665 693 if self.endDate!=None:
666 694 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
667 695 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
668 696 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
669 697 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
670 698 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
671 699 if self.timezone == 'lt':
672 700 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
673 701 if (startDateTime_File>endDateTime_Reader):
674 702 self.flag_standby = False
675 703 return 0
676 704 if self.flag_ignoreFiles and (startDateTime_File >= self.ignStartDateTime and startDateTime_File <= self.ignEndDateTime):
677 705 print("Ignoring...")
678 706 self.flag_standby = True
679 707 return 1
680 708 self.flag_standby = False
681 709
682 710 self.jrodataset = self.reshapeData()
683 711 #----self.updateIndexes()
684 712 self.profileIndex = 0
685 713
686 714 return 1
687 715
688 716
689 717 def __hasNotDataInBuffer(self):
690 718 if self.profileIndex >= (self.newProfiles*self.nblocks):
691 719 return 1
692 720 return 0
693 721
694 722
695 723 def getData(self):
696 724
697 725 if self.flagNoMoreFiles:
698 726 self.dataOut.flagNoData = True
699 727 return 0
700 728
701 729 if self.profileIndex >= (self.newProfiles*self.nblocks): #
702 730 #if self.__hasNotDataInBuffer():
703 731 if not (self.readNextFile(self.online)):
704 732 print("Profile Index break...")
705 733 return 0
706 734
707 735 if self.flag_standby: #Standby mode, if files are being ignoring, just return with no error flag
708 736 return 0
709 737
710 738 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
711 739 self.dataOut.flagNoData = True
712 740 print("No more data break...")
713 741 return 0
714 742
715 743 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
716 744
717 745 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
718 746
719 747 #print("R_t",self.timeset)
720 748
721 749 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
722 750 #verificar basic header de jro data y ver si es compatible con este valor
723 751 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
724 752 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
725 753 indexblock = self.profileIndex/self.newProfiles
726 754 #print (indexblock, indexprof)
727 755 diffUTC = 0
728 756 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
729 757
730 758 #print("utc :",indexblock," __ ",t_comp)
731 759 #print(numpy.shape(self.timeset))
732 760 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
733 761 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
734 762
735 763 self.dataOut.profileIndex = self.profileIndex
736 764 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
737 765 self.dataOut.flagNoData = False
738 766 # if indexprof == 0:
739 767 # print("kamisr: ",self.dataOut.utctime)
740 768
741 769 self.profileIndex += 1
742 770
743 771 return self.dataOut.data #retorno necesario??
744 772
745 773
746 774 def run(self, **kwargs):
747 775 '''
748 776 This method will be called many times so here you should put all your code
749 777 '''
750 778 #print("running kamisr")
751 779 if not self.isConfig:
752 780 self.setup(**kwargs)
753 781 self.isConfig = True
754 782
755 783 self.getData()
General Comments 0
You need to be logged in to leave comments. Login now