From d43c87bf928825c22310af53eb1e243bb8fd08de 2019-06-23 22:02:34 From: Juan C. Espinoza Date: 2019-06-23 22:02:34 Subject: [PATCH] Add skynoise with python --- diff --git a/plotter/templates/tools.html b/plotter/templates/tools.html index 92b0000..4f26989 100644 --- a/plotter/templates/tools.html +++ b/plotter/templates/tools.html @@ -18,12 +18,12 @@ aria-hidden="true"> + aria-describedby="validationTooltipSkynoiseDate" value="{% now 'd-m-Y' %}" required>
Please enter a valid date.
Go + data-image="{% url 'url_skynoise' %}">Go @@ -69,8 +69,8 @@ var image = $(e.relatedTarget).data('image'); if (image.indexOf('skynoise') > 0) { - var dt = $('#skynoise-date').val().split('/') - image += 'month=' + dt[0] + '&dom=' + dt[1] + '&year=' + dt[2]; + var dt = $('#skynoise-date').val(); + image += '?date=' + dt; } //populate values diff --git a/plotter/urls.py b/plotter/urls.py index 110b9a7..691f1e5 100644 --- a/plotter/urls.py +++ b/plotter/urls.py @@ -1,10 +1,11 @@ from django.conf.urls import url -from .views import main, plot, tools, reports +from .views import main, plot, tools, plot_skynoise, reports urlpatterns = [ url(r'^$', main, name='url_main'), url(r'^tools/$', tools, name='url_tools'), + url(r'^tools/skynoise.png$', plot_skynoise, name='url_skynoise'), url(r'^reports/$', reports, name='url_reports'), url(r'^(?P[0-9]+)/(?P[-\w]+)/$', plot, name='url-plot'), ] diff --git a/plotter/views.py b/plotter/views.py index 86b7f08..d36103c 100644 --- a/plotter/views.py +++ b/plotter/views.py @@ -10,11 +10,14 @@ from django import forms from django.contrib import messages from django.utils.safestring import mark_safe from django.shortcuts import render +from django.http import HttpResponse import mongoengine from plotter.models import Experiment, ExpDetail, PlotMeta, PlotData +from utils.plots import skynoise_plot + host = os.environ.get('HOST_MONGO', 'localhost') mongoengine.connect('dbplots', host=host, port=27017) @@ -168,4 +171,15 @@ def plot(request, code=None, plot=None): else: return render(request, 'home.html', {}) - \ No newline at end of file +def plot_skynoise(request): + + date = request.GET.get('date', None) + if date is None: + date = datetime.now() + else: + date = datetime.strptime(date, '%d-%m-%Y') + + data = skynoise_plot(date.year, date.month, date.day) + response = HttpResponse(data.getvalue(), content_type='image/png') + + return response diff --git a/requirements.txt b/requirements.txt index 3a62801..33d092d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,4 +7,7 @@ pymongo pyzmq redis requests -simplejson \ No newline at end of file +simplejson +numpy +matplotlib +scipy \ No newline at end of file diff --git a/utils/TimeTools.py b/utils/TimeTools.py new file mode 100755 index 0000000..e74d553 --- /dev/null +++ b/utils/TimeTools.py @@ -0,0 +1,430 @@ +""" +The TIME_CONVERSIONS.py module gathers classes and functions for time system transformations +(e.g. between seconds from 1970 to datetime format). + +MODULES CALLED: +NUMPY, TIME, DATETIME, CALENDAR + +MODIFICATION HISTORY: +Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Aug 13, 2009. +""" + +import numpy +import time +from datetime import datetime as dt +import calendar + +class Time: + """ + time(year,month,dom,hour,min,secs) + + An object represents a date and time of certain event.. + + Parameters + ---------- + YEAR = Number of the desired year. Year must be valid values from the civil calendar. + Years B.C.E must be represented as negative integers. Years in the common era are repre- + sented as positive integers. In particular, note that there is no year 0 in the civil + calendar. 1 B.C.E. (-1) is followed by 1 C.E. (1). + + MONTH = Number of desired month (1=Jan, ..., 12=December). + + DOM = Number of day of the month. + + HOUR = Number of the hour of the day. By default hour=0 + + MINS = Number of the minute of the hour. By default min=0 + + SECS = Number of the second of the minute. By default secs=0. + + Examples + -------- + time_info = time(2008,9,30,12,30,00) + + time_info = time(2008,9,30) + """ + + def __init__(self, year=None, month=None, dom=None, hour=0, mins=0, secs=0): + # If one the first three inputs are not defined, it takes the current date. + date = time.localtime() + if year==None:year=date[0] + if month==None:month=date[1] + if dom==None:dom=date[2] + + # Converting to arrays + year = numpy.array([year]); month = numpy.array([month]); dom = numpy.array([dom]) + hour = numpy.array([hour]); mins = numpy.array([mins]); secs = numpy.array([secs]) + + # Defining time information object. + self.year = numpy.atleast_1d(year) + self.month = numpy.atleast_1d(month) + self.dom = numpy.atleast_1d(dom) + self.hour = numpy.atleast_1d(hour) + self.mins = numpy.atleast_1d(mins) + self.secs = numpy.atleast_1d(secs) + + def change2julday(self): + """ + Converts a datetime to Julian days. + """ + + # Defining constants + greg = 2299171 # incorrect Julian day for Oct, 25, 1582. + min_calendar = -4716 + max_calendar = 5000000 + + min_year = numpy.nanmin(self.year) + max_year = numpy.nanmax(self.year) + if (min_yearmax_calendar): + print ("Value of Julian date is out of allowed range") + return -1 + + noyear = numpy.sum(self.year==0) + if noyear>0: + print ("There is no year zero in the civil calendar") + return -1 + + # Knowing if the year is less than 0. + bc = self.year<0 + + # Knowing if the month is less than March. + inJanFeb = self.month<=2 + + jy = self.year + bc - inJanFeb + jm = self.month + (1 + 12*inJanFeb) + + # Computing Julian days. + jul= numpy.floor(365.25*jy) + numpy.floor(30.6001*jm) + (self.dom+1720995.0) + + # Test whether to change to Gregorian Calendar + if numpy.min(jul) >= greg: + ja = numpy.int32(0.01*jy) + jul = jul + 2 - ja + numpy.int32(0.25*ja) + else: + gregchange = numpy.where(jul >= greg) + if gregchange[0].size>0: + ja = numpy.int32(0.01 + jy[gregchange]) + jy[gregchange] = jy[gregchange] + 2 - ja + numpy.int32(0.25*ja) + + # Determining machine-specific parameters affecting floating-point. + eps = 0.0 # Replace this line for a function to get precision. + eps = abs(jul)*0.0 > eps + + jul = jul + (self.hour/24. -0.5) + (self.mins/1440.) + (self.secs/86400.) + eps + + return jul[0] + + def change2secs(self): + """ + Converts datetime to number of seconds respect to 1970. + """ + + dtime = dt(self.year, self.month, self.dom) + return (dtime-dt(1970, 1, 1)).total_seconds() + + year = self.year + if year.size>1: year = year[0] + + month = self.month + if month.size>1: month = month[0] + + dom = self.dom + if dom.size>1: dom = dom[0] + + # Resizing hour, mins and secs if it was necessary. + hour = self.hour + if hour.size>1:hour = hour[0] + if hour.size==1:hour = numpy.resize(hour,year.size) + + mins = self.mins + if mins.size>1:mins = mins[0] + if mins.size==1:mins = numpy.resize(mins,year.size) + + secs = self.secs + if secs.size>1:secs = secs[0] + if secs.size==1:secs = numpy.resize(secs,year.size) + + # Using time.mktime to compute seconds respect to 1970. + secs1970 = numpy.zeros(year.size) + for ii in numpy.arange(year.size): + secs1970[ii] = time.mktime((int(year[ii]),int(month[ii]),int(dom[ii]),\ + int(hour[ii]),int(mins[ii]),int(secs[ii]),0,0,0)) + + secs1970 = numpy.int32(secs1970 - time.timezone) + + return secs1970 + + def change2strdate(self,mode=1): + """ + change2strdate method converts a date and time of certain event to date string. The + string format is like localtime (e.g. Fri Oct 9 15:00:19 2009). + + Parameters + ---------- + None. + + Return + ------ + + Modification History + -------------------- + Created by Freddy R. Galindo, ROJ, 09 October 2009. + + """ + + secs = numpy.atleast_1d(self.change2secs()) + strdate = [] + for ii in numpy.arange(numpy.size(secs)): + secs_tmp = time.localtime(secs[ii] + time.timezone) + if mode==1: + strdate.append(time.strftime("%d-%b-%Y (%j) %H:%M:%S",secs_tmp)) + elif mode==2: + strdate.append(time.strftime("%d-%b-%Y (%j)",secs_tmp)) + + strdate = numpy.array(strdate) + + return strdate + + +class Secs: + """ + secs(secs): + + An object represents the number of seconds respect to 1970. + + Parameters + ---------- + + SECS = A scalar or array giving the number of seconds respect to 1970. + + Example: + -------- + secs_info = secs(1251241373) + + secs_info = secs([1251241373,1251241383,1251241393]) + """ + def __init__(self,secs): + self.secs = secs + + def change2julday(self): + """ + Convert seconds from 1970 to Julian days. + """ + + secs_1970 = time(1970,1,1,0,0,0).change2julday() + + julian = self.secs/86400.0 + secs_1970 + + return julian + + def change2time(self): + """ + Converts seconds from 1970 to datetime. + """ + + secs1970 = numpy.atleast_1d(self.secs) + + datetime = numpy.zeros((9,secs1970.size)) + for ii in numpy.arange(secs1970.size): + tuple = time.gmtime(secs1970[ii]) + datetime[0,ii] = tuple[0] + datetime[1,ii] = tuple[1] + datetime[2,ii] = tuple[2] + datetime[3,ii] = tuple[3] + datetime[4,ii] = tuple[4] + datetime[5,ii] = tuple[5] + datetime[6,ii] = tuple[6] + datetime[7,ii] = tuple[7] + datetime[8,ii] = tuple[8] + + datetime = numpy.int32(datetime) + + return datetime + + +class Julian: + """ + julian(julian): + + An object represents julian days. + + Parameters + ---------- + + JULIAN = A scalar or array giving the julina days. + + Example: + -------- + julian_info = julian(2454740) + + julian_info = julian([2454740,2454760,2454780]) + """ + def __init__(self,julian): + self.julian = numpy.atleast_1d(julian) + + def change2time(self): + """ + change2time method converts from julian day to calendar date and time. + + Return + ------ + year = An array giving the year of the desired julian day. + month = An array giving the month of the desired julian day. + dom = An array giving the day of the desired julian day. + hour = An array giving the hour of the desired julian day. + mins = An array giving the minute of the desired julian day. + secs = An array giving the second of the desired julian day. + + Examples + -------- + >> jd = 2455119.0 + >> [yy,mo,dd,hh,mi,ss] = TimeTools.julian(jd).change2time() + >> print [yy,mo,dd,hh,mi,ss] + [2009] [10] [ 14.] [ 12.] [ 0.] [ 0.] + + Modification history + -------------------- + Translated from "Numerical Recipies in C", by William H. Press, Brian P. Flannery, + Saul A. Teukolsky, and William T. Vetterling. Cambridge University Press, 1988. + Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. + """ + + min_julian = -1095 + max_julian = 1827933925 + if (numpy.min(self.julian) < min_julian) or (numpy.max(self.julian) > max_julian): + print ('Value of Julian date is out of allowed range.') + return None + + # Beginning of Gregorian calendar + igreg = 2299161 + julLong = numpy.floor(self.julian + 0.5) + minJul = numpy.min(julLong) + + if (minJul >= igreg): + # All are Gregorian + jalpha = numpy.int32(((julLong - 1867216) - 0.25)/36524.25) + ja = julLong + 1 + jalpha - numpy.int32(0.25*jalpha) + else: + ja = julLong + gregChange = numpy.where(julLong >= igreg) + if gregChange[0].size>0: + jalpha = numpy.int32(((julLong[gregChange]-1867216) - 0.25)/36524.25) + ja[gregChange] = julLong[gregChange]+1+jalpha-numpy.int32(0.25*jalpha) + + # clear memory. + jalpha = -1 + + jb = ja + 1524 + jc = numpy.int32(6680. + ((jb-2439870)-122.1)/365.25) + jd = numpy.int32(365.*jc + (0.25*jc)) + je = numpy.int32((jb - jd)/30.6001) + + dom = jb - jd - numpy.int32(30.6001*je) + month = je - 1 + month = ((month - 1) % 12) + 1 + month = numpy.atleast_1d(month) + year = jc - 4715 + year = year - (month > 2)*1 + year = year - (year <= 0)*1 + year = numpy.atleast_1d(year) + + # Getting hours, minutes, seconds + fraction = self.julian + 0.5 - julLong + eps_0 = dom*0.0 + 1.0e-12 + eps_1 = 1.0e-12*numpy.abs(julLong) + eps = (eps_0>eps_1)*eps_0 + (eps_0<=eps_1)*eps_1 + + hour_0 = dom*0 + 23 + hour_2 = dom*0 + 0 + hour_1 = numpy.floor(fraction*24.0 + eps) + hour = ((hour_1>hour_0)*23) + ((hour_1<=hour_0)*hour_1) + hour = ((hour_1=hour_2)*hour_1) + + fraction = fraction - (hour/24.0) + mins_0 = dom*0 + 59 + mins_2 = dom*0 + 0 + mins_1 = numpy.floor(fraction*1440.0 + eps) + mins = ((mins_1>mins_0)*59) + ((mins_1<=mins_0)*mins_1) + mins = ((mins_1=mins_2)*mins_1) + + secs_2 = dom*0 + 0 + secs_1 = (fraction - mins/1440.0)*86400.0 + secs = ((secs_1=secs_2)*secs_1) + + return year,month,dom,hour,mins,secs + + def change2secs(self): + """ + Converts from Julian days to seconds from 1970. + """ + + jul_1970 = Time(1970,1,1,0,0,0).change2julday() + + secs = numpy.int32((self.julian - jul_1970)*86400) + + return secs + + def change2lst(self,longitude=-76.8667): + """ + CT2LST converts from local civil time to local mean sideral time + + longitude = The longitude in degrees (east of Greenwich) of the place for which + the local sideral time is desired, scalar. The Greenwich mean sideral time (GMST) + can be found by setting longitude=0. + """ + + # Useful constants, see Meus, p. 84 + c = numpy.array([280.46061837, 360.98564736629, 0.000387933, 38710000.0]) + jd2000 = 2451545.0 + t0 = self.julian - jd2000 + t = t0/36525. + + # Computing GST in seconds + theta = c[0] + (c[1]*t0) + (t**2)*(c[2]-t/c[3]) + + # Computing LST in hours + lst = (theta + longitude)/15.0 + neg = numpy.where(lst < 0.0) + if neg[0].size>0:lst[neg] = 24.0 + (lst[neg] % 24) + lst = lst % 24.0 + + return lst + + +class date2doy: + def __init__(self,year,month,day): + self.year = year + self.month = month + self.day = day + + def change2doy(self): + if calendar.isleap(self.year) == True: + tfactor = 1 + else: + tfactor = 2 + + day = self.day + month = self.month + + doy = numpy.floor((275*month)/9.0) - (tfactor*numpy.floor((month+9)/12.0)) + day - 30 + + return numpy.int32(doy) + + +class Doy2Date: + def __init__(self,year,doy): + self.year = year + self.doy = doy + + def change2date(self): + months = numpy.arange(12) + 1 + + first_dem = date2doy(self.year,months,1) + first_dem = first_dem.change2doy() + + imm = numpy.where((self.doy - first_dem) > 0) + + month = imm[0].size + dom = self.doy -first_dem[month - 1] + 1 + + return month, dom diff --git a/utils/__init__.py b/utils/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/utils/__init__.py diff --git a/utils/galaxy.txt b/utils/galaxy.txt new file mode 100755 index 0000000..9ba726c --- /dev/null +++ b/utils/galaxy.txt @@ -0,0 +1,99 @@ +0000 13.00 +0015 13.00 +0030 13.00 +0045 13.00 +0060 13.00 +0075 12.50 +0090 12.50 +0105 12.50 +0120 11.50 +0135 11.00 +0150 10.50 +0165 10.00 +0180 10.00 +0195 10.00 +0210 9.00 +0225 8.50 +0240 8.50 +0255 8.50 +0270 8.50 +0285 8.00 +0300 8.00 +0315 8.50 +0330 9.00 +0345 10.00 +0360 11.00 +0375 12.50 +0390 14.50 +0405 17.00 +0420 18.00 +0435 17.00 +0450 15.50 +0465 15.00 +0480 14.00 +0495 12.50 +0510 11.00 +0525 10.00 +0540 9.50 +0555 9.00 +0570 23.00 +0585 8.00 +0600 10.00 +0615 10.50 +0630 10.00 +0645 9.00 +0660 8.50 +0675 9.00 +0690 10.00 +0705 11.00 +0720 12.00 +0735 12.50 +0750 13.50 +0765 13.00 +0780 13.00 +0795 13.00 +0810 13.00 +0825 12.50 +0840 12.00 +0855 12.50 +0870 13.00 +0885 14.00 +0900 15.00 +0915 17.00 +0930 18.00 +0945 17.50 +0960 16.50 +0975 16.50 +0990 17.00 +0990 17.00 +1005 17.00 +1020 20.00 +1035 26.00 +1050 30.00 +1065 36.00 +1080 47.00 +1095 71.00 +1102 60.00 +1110 77.00 +1115 87.00 +1120 83.00 +1130 60.00 +1140 50.00 +1155 35.00 +1170 28.00 +1185 21.00 +1200 18.00 +1215 16.00 +1237 15.50 +1260 15.00 +1275 15.50 +1290 16.00 +1305 15.50 +1320 15.00 +1335 14.50 +1350 14.00 +1365 13.00 +1380 12.00 +1395 12.50 +1410 13.00 +1425 12.50 diff --git a/utils/plots.py b/utils/plots.py new file mode 100644 index 0000000..73d8e90 --- /dev/null +++ b/utils/plots.py @@ -0,0 +1,188 @@ + +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)) & (lstnumpy.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 +