plots.py
188 lines
| 5.9 KiB
| text/x-python
|
PythonLexer
/ utils / plots.py
|
r30 | |||
import os | ||||
import time | ||||
import numpy | ||||
import scipy | ||||
import base64 | ||||
from io import BytesIO | ||||
import scipy.interpolate | ||||
from matplotlib.figure import Figure | ||||
from utils import TimeTools | ||||
def skyNoise(jd, ut=-5.0, longitude=-76.87, filename='/app/utils/galaxy.txt'): | ||||
""" | ||||
hydrapos returns RA and Dec provided by Bill Coles (Oct 2003). | ||||
Parameters | ||||
---------- | ||||
jd = The julian date of the day (and time), scalar or vector. | ||||
dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica- | ||||
marca. The default value is -11.95. | ||||
filename = A string to specify name the skynoise file. The default value is skynoi- | ||||
se_jro.dat | ||||
Return | ||||
------ | ||||
maxra = The maximum right ascension to the declination used to get a cut. | ||||
ra = The right ascension. | ||||
Examples | ||||
-------- | ||||
>> [maxra,ra] = skynoise_jro() | ||||
>> print maxra, ra | ||||
139.45 -12.0833333333 | ||||
Modification history | ||||
-------------------- | ||||
Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | ||||
""" | ||||
# Defining date to compute SkyNoise. | ||||
[year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time() | ||||
is_dom = (month==9) & (dom==21) | ||||
if is_dom: | ||||
tmp = jd | ||||
jd = TimeTools.Time(year,9,22).change2julian() | ||||
dom = 22 | ||||
# Reading SkyNoise | ||||
g = os.getcwd() | ||||
f = open(filename) | ||||
lines = f.read() | ||||
f.close() | ||||
nlines = 99 | ||||
lines = lines.split('\n') | ||||
data = numpy.zeros((2,nlines))*numpy.float32(0.) | ||||
for ii in numpy.arange(nlines): | ||||
line = numpy.array([lines[ii][0:6],lines[ii][6:]]) | ||||
data[:,ii] = numpy.float32(line) | ||||
# Getting SkyNoise to the date desired. | ||||
otime = data[0,:]*60.0 | ||||
opowr = data[1,:] | ||||
hour = numpy.array([0,23]) | ||||
mins = numpy.array([0,59]) | ||||
secs = numpy.array([0,59]) | ||||
LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday() | ||||
LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1)) | ||||
lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst() | ||||
ipowr = lst*0.0 | ||||
# Interpolating using scipy (inside max and min "x") | ||||
otime = otime/3600. | ||||
val = numpy.where((lst>numpy.min(otime)) & (lst<numpy.max(otime))); val = val[0] | ||||
tck = scipy.interpolate.interp1d(otime,opowr) | ||||
ipowr[val] = tck(lst[val]) | ||||
# Extrapolating above maximum time data (23.75). | ||||
uval = numpy.where(lst>numpy.max(otime)) | ||||
if uval[0].size>0: | ||||
ii = numpy.min(uval[0]) | ||||
m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2]) | ||||
b = ipowr[ii-1] - m*lst[ii-1] | ||||
ipowr[uval] = m*lst[uval] + b | ||||
if is_dom: | ||||
lst = numpy.roll(lst,4) | ||||
ipowr = numpy.roll(ipowr,4) | ||||
new_lst = numpy.int32(lst*3600.) | ||||
new_pow = ipowr | ||||
return ipowr, LTtime, lst | ||||
def skynoise_plot(year, month, day): | ||||
""" | ||||
getPlot method creates a skynoise map over Jicamarca for a desired day. Additionally | ||||
save a PNG file of this plot. | ||||
Examples | ||||
-------- | ||||
>> SkyNoisePlot(skypower,skytime,skytime_lst).getPlot() | ||||
Modification history | ||||
-------------------- | ||||
Created by Freddy R. Galindo, ROJ, 18 October 2009. | ||||
""" | ||||
# Working with the time before to plot the SkyNoise | ||||
julian = TimeTools.Time(year, month, day).change2julday() | ||||
power, tm, time_lst = skyNoise(julian) | ||||
secs = TimeTools.Julian(tm).change2secs() | ||||
secs_lst = time_lst*1. | ||||
if time_lst.argmin()>0: | ||||
secs_lst[time_lst.argmin():] = secs_lst[time_lst.argmin():] + 24. | ||||
secs_lst = secs_lst*3600. | ||||
label_secs = time.localtime(secs[power.argmax()] + time.timezone) | ||||
label_secs_lst = time.localtime(secs_lst[power.argmax()] + time.timezone) | ||||
# Creating labels for main x-labelticks (Local Time): | ||||
snow = TimeTools.Time(year, month, day, 0, 0, 0).change2secs() | ||||
today = secs - snow | ||||
xticks_dpos = [] | ||||
xticks_dval = [] | ||||
for ix in [0,120,240,360,480,600,720,840,960,1080,1200,1320,1439]: | ||||
xticks_dpos.append(today[ix]) | ||||
time_tuple = time.localtime(today[ix] + time.timezone) | ||||
xticks_dval.append(time.strftime('%H:%M',time_tuple)) | ||||
if ix==1439:xticks_dval[12] = '' | ||||
# Creating labels for secondary x-labelticks (Sidereal Time): | ||||
xticks_upos = [secs_lst[0],secs_lst[-1]] | ||||
xticks_uval = ['',''] | ||||
for ix in [0,120,240,360,480,600,720,840,960,1080,1200,1320]: | ||||
index = numpy.argmin(numpy.abs(time_lst - (ix/60. + 1.))) | ||||
xticks_upos.append(secs_lst[index]) | ||||
time_tuple = time.localtime(secs_lst[index] + time.timezone) | ||||
if time_tuple[4]==59: | ||||
tmp_list = list(time_tuple) | ||||
tmp_list[4] = 0 | ||||
tmp_list[3] = tmp_list[3] + 1 | ||||
time_tuple = tuple(tmp_list) | ||||
xticks_uval.append(time.strftime('%H:%M',time_tuple)) | ||||
# Creating SkyNoise Map. | ||||
fig = Figure() | ||||
fig.subplots_adjust(top=0.80) | ||||
main_mesg = "SKY BRIGHTNESS AT 50Mhz - Date: " | ||||
main_mesg = main_mesg + time.strftime("%d-%b-%Y (%j)",label_secs) | ||||
main_mesg = main_mesg + "\nGalaxy Pass at " + time.strftime("%H:%M:%S LT",label_secs) | ||||
main_mesg = main_mesg + time.strftime(" (%H:%M:%S LST)",label_secs_lst) | ||||
fig.suptitle(main_mesg, fontsize=12) | ||||
axl = fig.add_subplot(1,1,1) | ||||
axl.plot(today,power,'b-') | ||||
axr = axl.twinx() | ||||
axl.set_xlabel('Local Time (Hours)') | ||||
axl.set_ylabel('Signal Strengh (mA)') | ||||
axr.set_ylabel('Temperature [K]x10^3') | ||||
axl.set_xlim(0,24) | ||||
axl.set_ylim(0,115) | ||||
axr.set_ylim(0,50) | ||||
axl.set_xticks(xticks_dpos) | ||||
axl.set_xticklabels(xticks_dval, fontsize=8) | ||||
axl.grid() | ||||
axt = axl.twiny() | ||||
axt.set_xlabel('Local Sidereal Time (Hours)') | ||||
axt.set_xlim(secs_lst[0],secs_lst[-1]) | ||||
axt.set_xticks(xticks_upos) | ||||
axt.set_xticklabels(xticks_uval, fontsize=8) | ||||
buf = BytesIO() | ||||
fig.savefig(buf, format="png") | ||||
return buf | ||||