From 20fd882fbe8ce710c9bd722eb2bfed8bb20bec18 Mon Sep 17 00:00:00 2001 From: Malgosia Sobolewska Date: Fri, 21 Jan 2022 14:34:40 -0500 Subject: [PATCH 1/4] Major update of script building ACE html, remove alerts into a separate future script. --- ACE/Scripts/create_ace_html_page.py | 898 ++++++---------------------- 1 file changed, 187 insertions(+), 711 deletions(-) diff --git a/ACE/Scripts/create_ace_html_page.py b/ACE/Scripts/create_ace_html_page.py index e8ecbdf..f61b95f 100755 --- a/ACE/Scripts/create_ace_html_page.py +++ b/ACE/Scripts/create_ace_html_page.py @@ -1,769 +1,245 @@ #!/usr/bin/env /data/mta4/Script/Python3.8/envs/ska3-shiny/bin/python -##################################################################################################### -# # -# create_ace_html_page.py: read ace data and update html page # -# # -# author: t. isobe (tisobe@cfa.harvard.edu) # -# # -# Last update: Nov 04, 2021 # -# # -##################################################################################################### - import os -import sys import re -import string -import random -import operator -import time -import datetime -import numpy -import Chandra.Time -# -#--- reading directory list -# -path = '/data/mta4/Space_Weather/house_keeping/dir_list' - -with open(path, 'r') as f: - data = [line.strip() for line in f.readlines()] - -for ent in data: - atemp = re.split(':', ent) - var = atemp[1].strip() - line = atemp[0].strip() - exec("%s = %s" %(var, line)) -# -#--- append pathes to private folders to a python directory -# -sys.path.append('/data/mta4/Script/Python3.8/MTA/') -# -#--- import several functions -# -import mta_common_functions as mcf -# -#--- temp writing file name -# -rtail = int(time.time() * random.random()) -zspace = '/tmp/zspace' + str(rtail) -# -#--- set direcotries -# -data_dir = ace_dir + 'Data/' -templ_dir = ace_dir + 'Scripts/Template/' -web_dir = html_dir + 'ACE/' -plot_dir = web_dir + 'Plots/' -# -#--- ftp, html site setting -# -http_epam = 'http://' + noaa_site + 'images/ace-epam-7-day.gif' -http_mag = 'http://' + noaa_site + 'images/ace-mag-swepam-7-day.gif' -# -#--- set parameters -# -p5_p3_scale = 7. #--- scale P5 to P3 values, while P3 is broke -p6_p3_scale = 36. #--- scale P6 to P3 values, while P3 is broke -p7_p3_scale = 110. #--- scale P7 to P3 values, while P3 is broke +import argparse +import numpy as np +from Chandra.Time import DateTime +from astropy.table import Table, unique -#--------------------------------------------------------------------------------------------------- -#-- create_ace_html_page: read ace data and update html page --- -#--------------------------------------------------------------------------------------------------- -def create_ace_html_page(): - """ - read ace data and update html page - input: none, but read from: - http://services.swpc.noaa.gov/images/ace-epam-7-day.gif - http://services.swpc.noaa.gov/images/images/ace-mag-swepam-7-day.gif - /Data/ace_12h_archive - output: /ACE/ace.html - """ -# -#---- read 12h_archive data -# - ifile = data_dir + 'ace_12h_archive' - cdata = mcf.read_data_file(ifile) -# -#--- cdata: a list of lists of electron/proton flux data -#--- l_vals: a list of the 'last' entries of those electron/proton flux data -#--- data is also trimmed to the last 2 hours -# - try: - cdata_cols, l_vals = convert_to_col_data(cdata) - except: - exit(1) -# -#--- create data table for the html page -# - ace_table = create_ace_data_table(cdata_cols, l_vals) -# -#--- download images and reverse the color: -# - download_img(http_epam) - download_img(http_mag) -# -#--- update ace html page -# - ifile = templ_dir + 'header' - line = read_file(ifile) - line = line + '\n' - - ifile = templ_dir + 'header1' - line = line + read_file(ifile) - - line = line + ace_table - line = line + '\n' - - ifile = templ_dir + 'image2' - line = line + read_file(ifile) - line = line + '\n' - - ifile = templ_dir + 'footer' - line = line + read_file(ifile) - line = line + '\n' - - out = web_dir + 'ace.html' - with open(out, 'w') as fo: - fo.write(line) +# ftp, html sites with ACE data -#--------------------------------------------------------------------------------------------------- -#-- create_ace_data_table: craete data tables from a given data list -- -#--------------------------------------------------------------------------------------------------- +HTTP_EPAM = "http://services.swpc.noaa.gov/images/ace-epam-7-day.gif" +HTTP_MAG = "http://services.swpc.noaa.gov/images/ace-mag-swepam-7-day.gif" -def create_ace_data_table(cdata, l_vals): - """ - craete data tables from a given data list - input: cdata --- a list of lists of data - [atime, jtime, echk, ech1, ech2, pchk, pch2, pch3, pch5, pch6, pch7] - l_vals --- a list of the 'last'entry values - [ech1_last, ech2_last, pch2_last, pch3_last, pch5_last, pch6_last, pch7_last]] - output: line --- table +# Set ACE P3 scaling parameters - """ - c_len = len(cdata[0]) -# -#--- last valid entry records --- these are values just outside of the current data period -# - p2_last = l_vals[2] - p3_last = l_vals[3] - p5_last = l_vals[4] - p6_last = l_vals[5] - p7_last = l_vals[6] -# -#--- initialize scaled lists -# - p5_p3_scaled = [] - p6_p3_scaled = [] - p7_p3_scaled = [] - for k in range(0, c_len): - p5_p3_scaled.append(-10000) - p6_p3_scaled.append(-10000) - p7_p3_scaled.append(-10000) -# -#--- for ease of read... -# - ctime = cdata[0] #--- Chandra Time - jtime = cdata[1] #--- display time - echk = numpy.array(cdata[2]) #--- whether electron data is good - de1 = numpy.array(cdata[3]) - de4 = numpy.array(cdata[4]) - pchk = numpy.array(cdata[5]) #--- whether proton data is good - p2dat = numpy.array(cdata[6]) - p3dat = numpy.array(cdata[7]) - p5dat = numpy.array(cdata[8]) - p6dat = numpy.array(cdata[9]) - p7dat = numpy.array(cdata[10]) -# -#--- go through the data -# - line = '' - for k in range(0, c_len): -# -#--- pchk and p3 -# - if pchk[k] == 0: - plast = p3_last - p3_diff = p3dat[k] - p3_last - if p3dat[k] > 0 and p3_diff > 500000: - p3dat[k] = -999 - pchk[k] = 1 - - if p3dat[k] > 0: - p3_last = p3dat[k] -# -#---- p2 -# - p2_diff = p2dat[k] - p2_last - if p2dat[k] > 0 and p2_diff > 1000000: - p2dat[k] = -999 - pchk[k] = 1 - - if p2dat[k] > 0: - p2_last = p2dat[k] -# -#--- p5 -# - p5_diff = p5dat[k] - p5_last - if p5dat[k] and p5_diff > 40000: - p5dat[k] = -999 - pchk[k] = 1 - p5_p3_scaled[k] = -999 - - if p5dat[k] > 0: - p5_last = p5dat[k] -# -#--- p6 -# - p6_diff = p6dat[k] - p6_last - if p6dat[k] and p6_diff > 20000: - p6dat[k] = -999 - pchk[k] = 1 - p6_p3_scaled[k] = -999 - - if p6dat[k] > 0: - p6_last = p6dat[k] -# -#--- p7 -# - p7_diff = p7dat[k] - p7_last - if p7dat[k] > 0 and p7_diff> 10000: - p7dat[k] = -999 - pchk[k] = 1 - p7_p3_scaled[k] = -999 - - if p7dat[k] > 0: - p7_last = p7dat[k] -# -#--- replace the scaled data with valid data if the original values are good -# - if p5dat[k] > 0: - p5_p3_scaled[k] = p5dat[k] * p5_p3_scale - - if p6dat[k] > 0: - p6_p3_scaled[k] = p6dat[k] * p6_p3_scale - - if p7dat[k] > 0: - p7_p3_scaled[k] = p7dat[k] * p7_p3_scale -# -#--- if data is ok, add the line to the table data -# - if pchk[k] == 0: - aline = "%s%11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f\n"\ - % (cdata[1][k], de1[k], de4[k], p2dat[k], p3dat[k],\ - p5_p3_scaled[k], p6_p3_scaled[k], p5dat[k], p6dat[k], p7dat[k]) - - line = line + aline -# -#--- the table part is done, compute other entries -# -# -#--- compute mean and find min of p5, p6, and p7 -# - ind = p5dat > 0 - cdat = p5dat[ind] - if len(cdat) > 0: - p337a = numpy.mean(cdat) - p337m = min(cdat) - else: - p337a = 0 - p337m = 0 - - ind = p6dat > 0 - cdat = p6dat[ind] - if len(cdat) > 0: - p761a = numpy.mean(cdat) - p761m = min(cdat) - else: - p761a = 0 - p761m = 0 - - ind = p7dat > 0 - cdat = p7dat[ind] - if len(cdat) > 0: - p1073a = numpy.mean(cdat) - p1073m = min(cdat) - else: - p1073a = 0 - p1073m = 0 -# -#--- do the same for the scaled data -# - p5_p3_scaled = numpy.array(p5_p3_scaled) - p6_p3_scaled = numpy.array(p6_p3_scaled) - p7_p3_scaled = numpy.array(p5_p3_scaled) - p5_p3_scaled = p5_p3_scaled[p5_p3_scaled >0] - p6_p3_scaled = p6_p3_scaled[p6_p3_scaled >0] - p7_p3_scaled = p7_p3_scaled[p7_p3_scaled >0] - - if len(p5_p3_scaled) > 0: - p5_p3a = numpy.mean(p5_p3_scaled) - p5_p3m = numpy.min(p5_p3_scaled) - else: - p5_p3a = 0 - p5_p3m = 0 +P5_P3_SCALE = 7. # scale P5 to P3 values, while P3 is broke +P6_P3_SCALE = 36. # scale P6 to P3 values, while P3 is broke +P7_P3_SCALE = 110. # scale P7 to P3 values, while P3 is broke - if len(p6_p3_scaled) > 0: - p6_p3a = numpy.mean(numpy.array(p6_p3_scaled)) - p6_p3m = numpy.min(numpy.array(p6_p3_scaled)) - else: - p6_p3a = 0 - p6_p3m = 0 - if len(p7_p3_scaled) > 0: - p7_p3a = numpy.mean(numpy.array(p7_p3_scaled)) - p7_p3m = numpy.min(numpy.array(p7_p3_scaled)) - else: - p7_p3a = 0 - p7_p3m = 0 -# -#--- index to find good data -# - ind = pchk == 0 - ind2 = echk == 0 - -# -#--- check whether there are any good data left after removing bad ones -# - chk = len(echk[ind]) - chk2 = len(echk[ind2]) - - if (chk == 0) or (chk2 == 0): - line = line + '

' - line = line + " No Valid data for last 2 hours." - line = line + '

\n' - return line -# -#--- there are good data -# - p56 = p2dat[ind] - p56 = p56[p56 >0] - p130 = p3dat[ind] - p130 = p130[p130 > 0] - if len(p56) > 0: - p56a = numpy.mean(p56) - p56m = min(p56) - else: - p56a = 0 - p56m = 0 +def get_options(): + parser = argparse.ArgumentParser(description='ACE data and fluence') + parser.add_argument('--ace_dir', type=str, + default='/data/mta4/Space_Weather/ACE/', + help='Root directory for ACE data and scripts') + parser.add_argument('--html_dir', type=str, + default='.', + help='Root directory for ACE html output') + args = parser.parse_args() + return args - if len(p130) > 0: - p130a = numpy.mean(p130) - p130m = min(p56) - else: - p130a = 0 - p130m = 0 - e38 = de1[ind2] - e38 = e38[e38 > 0] - e175 = de4[ind2] - e175 = e175[e175 >0] +def create_ace_html_page(): + """ + Read 2h ACE data from local MTA 12h archive + and update MTA's ACE html page at /ACE/ace.html. + Download ACE images from: + http://services.swpc.noaa.gov/images/ace-epam-7-day.gif + http://services.swpc.noaa.gov/images/images/ace-mag-swepam-7-day.gif + """ - if len(e38) > 0: - e38a = numpy.mean(e38) - e38m = min(e38) - else: - e38a = 0 - e38m = 0 - if len(e175) > 0: - e175a = numpy.mean(e175) - e175m = min(e175) - else: - e175a = 0 - e175m = 0 -# -#--- compute fluence of 2 hours -# - e38f = e38a * 7200. - e175f = e175a * 7200. - p56f = p56a * 7200. - p130f = p130a * 7200. - p337f = p337a * 7200. - p761f = p761a * 7200. - p1073f = p1073a * 7200. - p5_p3f = p5_p3a * 7200. - p6_p3f = p6_p3a * 7200. - p7_p3f = p7_p3a * 7200. -# -#--- compute spectral indcies -# - if p337a > 0: - p3_p5 = p130a / p337a - else: - p3_p5 = 0 - if p761a > 0: - p3_p6 = p130a / p761a - else: - p3_p6 = 0 - if p761a > 0: - p5_p6 = p337a / p761a - else: - p5_p6 = 0 - if p1073a > 0: - p6_p7 = p761a / p1073a - else: - p6_p7 = 0 - if p1073a > 0: - p5_p7 = p337a / p1073a - else: - p337a = 0 -# -#--- violation check -# - if p130f > 360000000: - val = "%.4e" % p130f - ace_violation_protons(val) + args = get_options() - p5_p6_lim = 1.0 #----????? this is what given in the original, but use below - p5_p6_lim = 1.0e10 + # read 12h local MTA ACE archive data + ifile = f"{args.ace_dir}/Data/ace_12h_archive" - if (p5_p6 > p5_p6_lim) or (p5_p6 < 1): - speci = "%12.1f" % p5_p6 - speci_lim = "%8.1f" % p5_p6_lim + dat = read_2h_ace_data_from_12h_archive(ifile) - ace_invalid_spec(speci, speci_lim) -# -#--- crate a summary table -# - hfile = templ_dir + 'header2' - line = line + "\n" + read_file(hfile) + ace_table = '' - line = line + "%7s %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f\n"\ - % ("AVERAGE ", e38a, e175a, p56a, p130a, p5_p3a, p6_p3a, p337a, p761a, p1073a) + if len(dat) > 0: + # There are nominal quality data in the last 2 hours, + # so create a string, ace_table, with data for html display. + # However, some channels may still contain bad data + # indicated by -1.e5 values. These are going to be + # displayed but filtered out before computing the stats. - line = line + "%7s %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f %11.3f\n"\ - % ("MINIMUM ", e38m, e175m, p56m, p130m, p5_p3m, p6_p3m, p337m, p761m, p1073m) + # Add scaled P5 and P6 + dat['p5_scaled_p3'] = dat['p5'] * P5_P3_SCALE + dat['p6_scaled_p3'] = dat['p6'] * P6_P3_SCALE - line = line + "%7s %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n\n"\ - % ("FLUENCE ", e38f, e175f, p56f, p130f, p5_p3f, p6_p3f, p337f, p761f, p1073f) + # ace_table with data for html display - line = line + "%7s %11s %11.3f %11s %11.3f %11s %11.3f %11s %11.3f \n\n"\ - % ("SPECTRA ", "p3/p5", p3_p5, "p3/p6", p3_p6, "p5/p6", p5_p6, "p6/p7", p6_p7) + cols = ("yr", "mo", "da", "hhmm", "de1", "de4", + "p2", "p3", "p5", "p6", + "p5_scaled_p3", "p6_scaled_p3", "p7") + row_format = ' '.join(["{:s} " * 4, "{:11.3f} " * 9, "\n"]) - line = line + "%62s %4.1f\n"\ - % ("* This P3 channel is currently scaled from P5 data. P3* = P5 X ", p5_p3_scale) + for row in dat: + ace_table = ace_table + row_format.format(*row[cols]) - line = line + "%62s %4.1f\n"\ - % ("** This P3 channel is currently scaled from P6 data. P3** = P6 X ", p6_p3_scale) + # Add a header for the table with stats - line = line + "%62s %4.1f\n"\ - % ("*** This P3 channel (not shown) is currently scaled from P7 data. P3*** = P7 X ", p7_p3_scale) + with open(f"{args.ace_dir}/Scripts/Template/header2") as fp: + header = fp.read() - return line + ace_table = ace_table + "\n" + header -#--------------------------------------------------------------------------------------------------- -#-- ace_violation_protons: send out a warning email -- -#--------------------------------------------------------------------------------------------------- + # Compute electron/proton stats (mean, min, fluence) + means = {} + mins = {} + fluences_2h = {} + cols = ("de1", "de4", "p2", "p3", "p5", "p6", + "p5_scaled_p3", "p6_scaled_p3", "p7") -def ace_violation_protons(val): - """ - send out a warning email - input: val --- value - output: email sent out + for col in cols: + # Filter out negative (bad) values + ok = dat[col] > 0 + if len(dat[col][ok]) > 2: + # At least 2 measurements of nominal quality + # in the last 2 hours + means[col] = np.mean(dat[col][ok]) + mins[col] = np.min(dat[col][ok]) + fluences_2h[col] = np.mean(dat[col][ok]) * 7200 + else: + means[col] = -1e5 + mins[col] = -1e5 + fluences_2h[col] = -1e5 + + # Spectral indexes + spectra = {'p3_p5': means['p3'] / means['p5'], + 'p3_p6': means['p3'] / means['p6'], + 'p5_p6': means['p5'] / means['p6'], + 'p6_p7': means['p6'] / means['p7']} + + # Add table with stats for the display + + for stat, name, fmt in zip((means, mins, fluences_2h), + ('AVERAGE', 'MINIMUM', 'FLUENCE'), + ('{:11.3f}', '{:11.3f}', '{:11.4e}')): + vals = stat.values() + row_format = ' '.join([f"{name.ljust(15)}", f"{fmt} " * len(vals), "\n"]) + ace_table = ace_table + row_format.format(*list(vals)) + + + ace_table = ace_table + "\n" + "%7s %11s %11.3f %11s %11.3f %11s %11.3f %11s %11.3f \n\n"\ + % ("SPECTRA ", "p3/p5", spectra['p3_p5'], "p3/p6", spectra['p3_p6'], + "p5/p6", spectra['p5_p6'], "p6/p7", spectra['p6_p7']) + + ace_table = ace_table + "%62s %4.1f\n"\ + % ("* This P3 channel is currently scaled from P5 data. P3* = P5 X ", + P5_P3_SCALE) + + ace_table = ace_table + "%62s %4.1f\n"\ + % ("** This P3 channel is currently scaled from P6 data. P3** = P6 X ", + P6_P3_SCALE) + + ace_table = ace_table + "%62s %4.1f\n"\ + % ("*** This P3 channel (not shown) is currently scaled from P7 data. P3*** = P7 X ", + P7_P3_SCALE) + + else: + # No valid data in the last 2h + ace_table = """ +

+ No valid data for last 2 hours. +

+""" + + # download images and reverse the color: + download_img(img=HTTP_EPAM, + html_dir=args.html_dir, + ace_dir=args.ace_dir) + download_img(img=HTTP_MAG, + html_dir=args.html_dir, + ace_dir=args.ace_dir) + + # read in headers, footers and image html from template files + # and insert ace_table string with current data and stats + contents = [] + files = ['header', 'header1', 'image2', 'footer'] + for file in files: + with open(f"{args.ace_dir}/Scripts/Template/{file}") as fp: + contents.append(fp.read()) + + contents.insert(2, ace_table) + + # update ACE html page + with open(f"{args.html_dir}/ACE/ace.html", "w") as fo: + fo.write("\n".join(contents)) + + +def read_2h_ace_data_from_12h_archive(fname): """ - line = 'A Radiation violation of P3 (130KeV) has been observed by ACE\n' - line = line + 'Observed = ' + val + '\n' - line = line + '(limit = fluence of 3.6e8 particles/cm2-ster-MeV within 2 hours)\n' - line = line + 'see http://cxc.harvard.edu/mta/ace.html\n' -# -#--- add current spacecraft info -# - line = line + curr_state() -# -#--- read alert case -# - if os.path.isfile('/data/mta4/www/Snapshot/.scs107alert'): - line = line + 'The ACIS on-call person should review the data and call a telecon if necessary.\n' - line = line + 'This message sent to sot_ace_alert\n' - send_mail('ACE_p3', line, 'sot_ace_alert@cfa.harvard.edu') -# -#--- yellow alert case -# - else: - line = line + 'This message sent to sot_yellow_alert\n' - send_mail('ACE_p3', line, 'sot_yellow_alert@cfa.harvard.edu') - -#--------------------------------------------------------------------------------------------------- -#-- ace_invalid_spec: sending out a warning email -- -#--------------------------------------------------------------------------------------------------- + Reads 12h ACE archive and keeps only the last 2h + of nominal quality (both the e_status and p_status + flags equal 0). -def ace_invalid_spec(speci, speci_lim): + :param fname: ACE 12h local MTA archive + :returns: astropy Table with ACE data """ - sending out a warning email - input: speci --- the value reported - speci_lim --- the limit of the vluae - output: email sent out - """ -# -#--- check whether the mail is recently sent out -# - out = '/tmp/mta/prot_spec_violate' - if os.path.isfile(out): - cmd = 'date >> ' + out - os.system(cmd) -# -#--- if not, send the warning email -# - else: - line = " A spectral index violation of P5/P6 has been observed by ACE, " - line = line + "indicating a possibly invalid P5 channel.\n" - line = line + 'Observed = ' + speci - line = line + '(limit = ' + speci_lim + ')\n' - line = line + 'see http://cxc.harvard.edu/mta/ace.html\n' - line = line + 'This message sent to mtadude\n' - - send_mail('ACE_p5/p6', line, 'mtadude@cfa.harvard.edu') -# -#--- open a file indicating that the mail was sent -# - with open(out, 'w') as fo: - fo.write(line) -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- + # TODO: update scripts that create the 12h archive to + # include the column names + names = ("yr", "mo", "da", "hhmm", "mjd", "mjdsec", + "e_status", "de1", "de4", + "p_status", "p2", "p3", "p5", "p6", "p7", "anis_index") -def send_mail(subject, content, address): + formats = ('S4', 'S2', 'S2', 'S4', 'f8', 'f8', + 'i1', 'f8', 'f8', + 'i1', 'f8', 'f8', 'f8', 'f8', 'f8', 'i1') - fo = open(zspace, 'w') - fo.write(content) - fo.close() + dat = np.loadtxt(fname, dtype={'names': names, 'formats' : formats}) + dat = Table(dat) - cmd = 'cat ' + zspace + '|mailx -s "' + subject + '" ' + address - os.system(cmd) - - cmd = 'rm ' + zspace - os.system(cmd) + # remove duplicate lines; TODO check why they appear + dat = unique(dat, keys=("yr", "mo", "da", "hhmm")) -#--------------------------------------------------------------------------------------------------- -#-- curr_state: extract some satellite related information -- -#--------------------------------------------------------------------------------------------------- + # keep only the most recent two hours of nominal quality data + fits_fmt = [f"{yr}-{mo.zfill(2)}-{da}T{hhmm[:2]}:{hhmm[-2:]}:00.000" + for yr, mo, da, hhmm in zip(dat['yr'], dat['mo'], + dat['da'], dat['hhmm'])] -def curr_state(): - """ - extract some satellite related information - input: none - output: line --- lines which display the satellite information - """ -# -#--- read CRM file to get some information -# - sfile = crm_dir + 'CRMsummary.dat' - data = mcf.read_data_file(sfile) - - dirct = 'NA' - dist = 'NA' - if len(data) < 1: - dline = '' - else: - dline = data[0] - for ent in data: - mc = re.search('Geocentric Distance', ent) - if mc is not None: - atemp = re.split(':', ent) - btemp = re.split('\s+', atemp[1].strip()) - dist = btemp[0] - if btemp[1] == 'A': - dirct = 'Ascending' - else: - dirct = 'Descending' - - line = 'Currently ' + dirct + ' through ' + dist + 'km with ' + dline + '\n\n' -# -#-- get DSN contact info from ephem data -# - ctime = time.strftime("%Y:%j:%H:%M:%S", time.gmtime()) - start = Chandra.Time.DateTime(ctime).secs - stop = start + 3.0 * 86400. - - ifile = ephem_dir + 'dsn_summary.dat' - data = mcf.read_data_file(ifile) - data = data[2:] - line = line + data[0] + '\n' + data[1] + '\n' - for ent in data: - atemp = re.split('\s+', ent) - year = atemp[10] - yday = atemp[11] - stime = convert_to_stime(year, yday) - if stime >= start and stime < stop: - line = line + ent + '\n' - - line = line + '\n' - - return line + tt = DateTime(np.array(fits_fmt), format='fits').date + now_minus_2h = DateTime(DateTime() - 2 / 24).date -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- + # status: 0 = nominal, 4,6,7,8 = bad data, unable to process, 9 = no data + ok = (dat['e_status'] == 0) & (dat['p_status'] == 0) & (tt > now_minus_2h) + dat = dat[ok] -def convert_to_stime(year, yday): + return dat - atemp = re.split('\.', yday) - frac = float(atemp[1]) - val = 24. * frac - hh = int(val) - diff = val - hh - val = 60 * diff - mm = int(val) - diff = val - mm - ss = int(60 *diff) - htime = str(year) + ':' + atemp[0] + ':' + str(hh) + ':' + str(mm) + ':' + str(ss) - - stime = Chandra.Time.DateTime(htime).secs +def download_img(img, html_dir, ace_dir, chg=1): + """ + Download image from web site. Save image in {args.html_dir}/ACE/Plots - return stime + :param img: web address of the image + :param html_dir: root dir of space weather html products + :param ace_dir: space weather ACE directory + :param chg: if > 0, reverse the color + """ -#--------------------------------------------------------------------------------------------------- -#-- download_img: down load an image from web site -- -#--------------------------------------------------------------------------------------------------- + # get the name of output img file name -def download_img(file, chg=1): - """ - down load an image from web site - input: file --- image file address - chg --- if >0, reverse the color - output: / - """ -# -#--- get the name of output img file name -# - atemp = re.split('\/', file) + atemp = re.split('\/', img) ofile = atemp[-1] - oimg = plot_dir + ofile -# -#--- download the img -# - #cmd = 'lynx -source ' + file + '>' + oimg + oimg = f"{html_dir}/ACE/Plots/{ofile}" + + # download the img + try: - cmd = 'wget -q -O' + oimg + ' ' + file + cmd = f'wget -q -O {oimg} {img}' os.system(cmd) except: - mc = re.search('gif', oimg) + mc = re.search('gif', oimg) if mc is not None: - cmd = ' cp ' + house_keeping + 'no_plot.gif ' + oimg + cmd = f"cp {ace_dir}/Scripts/Template/no_plot.gif {oimg}" else: - cmd = ' cp ' + house_keeping + 'no_data.png ' + oimg + cmd = f"cp {ace_dir}/Scripts/Template/no_data.png {oimg}" os.system(cmd) return -# -#--- reverse the color of the image -# - if chg == 1: - cmd = 'convert -negate ' + oimg + ' ' + oimg - os.system(cmd) + # reverse the color of the image -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- -#--------------------------------------------------------------------------------------------------- - -def read_file(ifile): - - f = open(ifile, 'r') - data = f.read() - f.close() - - return data - - -#----------------------------------------------------------------------------- -#-- convert_to_col_data: read data into a list of lists -- -#----------------------------------------------------------------------------- + if chg == 1: + cmd = f'convert -negate {oimg} {oimg}' + os.system(cmd) -def convert_to_col_data(data): - """ - convert the list of row data into column data - input: a list of lists - output: a list of lists of: - atime --- a time in seconds from 1998.1.1 - jtime --- a string time - echk --- data quality of electron data - ech1 --- electron 38-53 - ech2 --- electron 175-315 - pchk --- data quality of proton - pch2 --- proton 47-68 - pch3 --- proton 115-195 - pch5 --- proton 310-580 - pch6 --- proton 795-1193 - pch7 --- proton 1060-1900 - - [ech1_last, ech2_last, pch2_last, pch3_last, pch5_last, pch6_last, pch7_last] - --- this is a list of the last valid flux data before the current data set - """ -# -#--- find the most recent entry time and set the cutting time to 2 hrs before that -# - atemp = re.split('\s+', data[-1]) - ltime = atemp[0] + ':' + atemp[1] + ':' + atemp[2] + ':' + atemp[3][0] + atemp[3][1] + ':' - ltime = ltime + atemp[3][2] + atemp[3][3] + ':00' - ltime = time.strftime('%Y:%j:%H:%M:%S', time.strptime(ltime, '%Y:%m:%d:%H:%M:%S')) - cut = int(Chandra.Time.DateTime(ltime).secs) - 2 * 3600.0 - 60.0 - - atime = [] - jtime = [] - echk = [] - ech1 = [] - ech2 = [] - pchk = [] - pch2 = [] - pch3 = [] - pch5 = [] - pch6 = [] - pch7 = [] - ptime = 0 - for ent in data: - atemp = re.split('\s+', ent) - clen = len(atemp) -# -#--- convert time in Chandra Time -# - ltime = atemp[0] + ':' + atemp[1] + ':' + atemp[2] + ':' + atemp[3][0] + atemp[3][1] + ':' - ltime = ltime + atemp[3][2] + atemp[3][3] + ':00' - ltime = time.strftime('%Y:%j:%H:%M:%S', time.strptime(ltime, '%Y:%m:%d:%H:%M:%S')) - stime = int(Chandra.Time.DateTime(ltime).secs) -# -#--- sometime, there are double entries; so remove those -# - if stime == ptime: - continue - else: - ptime = stime -# -#--- keep the record of 'previous' entry; only valid data -# - if stime < cut: - chk1 = int(float(atemp[6])) - chk2 = int(float(atemp[9])) - if chk1 == 0 and chk2 == 0: - ech1_last = float(atemp[7]) - ech2_last = float(atemp[8]) - pch2_last = float(atemp[10]) - pch3_last = float(atemp[11]) - pch5_last = float(atemp[12]) - pch6_last = float(atemp[13]) - pch7_last = float(atemp[14]) - continue -# -#--- save time part in a string format (YR MO DA HHMM Day Day) -# - ftime = atemp[0] + ' ' + atemp[1] + ' ' + atemp[2] + ' ' + atemp[3] -# ftime = ftime + '%8d%8d' % (float(atemp[4]), float(atemp[5])) - - atime.append(stime) - jtime.append(ftime) - echk.append(int(float(atemp[6]))) - ech1.append(float(atemp[7])) - ech2.append(float(atemp[8])) - pchk.append(int(float(atemp[9]))) - pch2.append(float(atemp[10])) - pch3.append(float(atemp[11])) - pch5.append(float(atemp[12])) - pch6.append(float(atemp[13])) - pch7.append(float(atemp[14])) - - return [atime, jtime, echk, ech1, ech2, pchk, pch2, pch3, pch5, pch6, pch7],\ - [ech1_last, ech2_last, pch2_last, pch3_last, pch5_last, pch6_last, pch7_last] #--------------------------------------------------------------------------------------------------- From 9045af13c3cf96ca7a2442cd45a1956c9cc63de1 Mon Sep 17 00:00:00 2001 From: Malgosia Sobolewska Date: Fri, 21 Jan 2022 21:00:59 -0500 Subject: [PATCH 2/4] Move stat computations to a separate function; unify names of scaled columns and constants. --- ACE/Scripts/create_ace_html_page.py | 95 +++++++++++++++++------------ 1 file changed, 56 insertions(+), 39 deletions(-) diff --git a/ACE/Scripts/create_ace_html_page.py b/ACE/Scripts/create_ace_html_page.py index f61b95f..9f19093 100755 --- a/ACE/Scripts/create_ace_html_page.py +++ b/ACE/Scripts/create_ace_html_page.py @@ -15,9 +15,9 @@ # Set ACE P3 scaling parameters -P5_P3_SCALE = 7. # scale P5 to P3 values, while P3 is broke -P6_P3_SCALE = 36. # scale P6 to P3 values, while P3 is broke -P7_P3_SCALE = 110. # scale P7 to P3 values, while P3 is broke +P3_P5_SCALE = 7. # scale P5 to P3 values, while P3 is broke +P3_P6_SCALE = 36. # scale P6 to P3 values, while P3 is broke +P3_P7_SCALE = 110. # scale P7 to P3 values, while P3 is broke def get_options(): @@ -58,14 +58,14 @@ def create_ace_html_page(): # displayed but filtered out before computing the stats. # Add scaled P5 and P6 - dat['p5_scaled_p3'] = dat['p5'] * P5_P3_SCALE - dat['p6_scaled_p3'] = dat['p6'] * P6_P3_SCALE + dat['p3_scaled_p5'] = dat['p5'] * P3_P5_SCALE + dat['p3_scaled_p6'] = dat['p6'] * P3_P6_SCALE # ace_table with data for html display cols = ("yr", "mo", "da", "hhmm", "de1", "de4", - "p2", "p3", "p5", "p6", - "p5_scaled_p3", "p6_scaled_p3", "p7") + "p2", "p3", "p3_scaled_p5", "p3_scaled_p6", + "p5", "p6", "p7") row_format = ' '.join(["{:s} " * 4, "{:11.3f} " * 9, "\n"]) for row in dat: @@ -78,58 +78,38 @@ def create_ace_html_page(): ace_table = ace_table + "\n" + header - # Compute electron/proton stats (mean, min, fluence) - means = {} - mins = {} - fluences_2h = {} - cols = ("de1", "de4", "p2", "p3", "p5", "p6", - "p5_scaled_p3", "p6_scaled_p3", "p7") - - for col in cols: - # Filter out negative (bad) values - ok = dat[col] > 0 - if len(dat[col][ok]) > 2: - # At least 2 measurements of nominal quality - # in the last 2 hours - means[col] = np.mean(dat[col][ok]) - mins[col] = np.min(dat[col][ok]) - fluences_2h[col] = np.mean(dat[col][ok]) * 7200 - else: - means[col] = -1e5 - mins[col] = -1e5 - fluences_2h[col] = -1e5 - - # Spectral indexes - spectra = {'p3_p5': means['p3'] / means['p5'], - 'p3_p6': means['p3'] / means['p6'], - 'p5_p6': means['p5'] / means['p6'], - 'p6_p7': means['p6'] / means['p7']} + # Calculate electron/proton stats (mean, min, 2h fluence, spectra) + + stats = calculate_stats(dat) # Add table with stats for the display - for stat, name, fmt in zip((means, mins, fluences_2h), + for stat, name, fmt in zip((stats['means'], stats['mins'], stats['fluences_2h']), ('AVERAGE', 'MINIMUM', 'FLUENCE'), ('{:11.3f}', '{:11.3f}', '{:11.4e}')): vals = stat.values() row_format = ' '.join([f"{name.ljust(15)}", f"{fmt} " * len(vals), "\n"]) ace_table = ace_table + row_format.format(*list(vals)) + # Add info on spectral indexes ace_table = ace_table + "\n" + "%7s %11s %11.3f %11s %11.3f %11s %11.3f %11s %11.3f \n\n"\ - % ("SPECTRA ", "p3/p5", spectra['p3_p5'], "p3/p6", spectra['p3_p6'], - "p5/p6", spectra['p5_p6'], "p6/p7", spectra['p6_p7']) + % ("SPECTRA ", "p3/p5", stats['spectra']['p3_p5'], + "p3/p6", stats['spectra']['p3_p6'], + "p5/p6", stats['spectra']['p5_p6'], + "p6/p7", stats['spectra']['p6_p7']) ace_table = ace_table + "%62s %4.1f\n"\ % ("* This P3 channel is currently scaled from P5 data. P3* = P5 X ", - P5_P3_SCALE) + P3_P5_SCALE) ace_table = ace_table + "%62s %4.1f\n"\ % ("** This P3 channel is currently scaled from P6 data. P3** = P6 X ", - P6_P3_SCALE) + P3_P6_SCALE) ace_table = ace_table + "%62s %4.1f\n"\ % ("*** This P3 channel (not shown) is currently scaled from P7 data. P3*** = P7 X ", - P7_P3_SCALE) + P3_P7_SCALE) else: # No valid data in the last 2h @@ -203,6 +183,43 @@ def read_2h_ace_data_from_12h_archive(fname): return dat +def calculate_stats(dat): + """ + Calculate statistics: means, mins, 2h fluences in each channel, + and spectral indexes. + :param dat: astropy Table with ACE data + :returns stats: dictionary with values being dictionaries of stats + for ACE channels + """ + + stats = {'means': {}, 'mins': {}, 'fluences_2h': {}, 'spectra': {}} + + cols = ("de1", "de4", "p2", "p3", "p3_scaled_p5", "p3_scaled_p5", + "p5", "p6", "p7") + + for col in cols: + # Filter out negative (bad) values + ok = dat[col] > 0 + if len(dat[col][ok]) > 2: + # At least 2 measurements of nominal quality + # in the last 2 hours + stats['means'][col] = np.mean(dat[col][ok]) + stats['mins'][col] = np.min(dat[col][ok]) + stats['fluences_2h'][col] = np.mean(dat[col][ok]) * 7200 + else: + stats['means'][col] = -1e5 + stats['mins'][col] = -1e5 + stats['fluences_2h'][col] = -1e5 + + # Spectral indexes + stats['spectra'] = {'p3_p5': stats['means']['p3'] / stats['means']['p5'], + 'p3_p6': stats['means']['p3'] / stats['means']['p6'], + 'p5_p6': stats['means']['p5'] / stats['means']['p6'], + 'p6_p7': stats['means']['p6'] / stats['means']['p7']} + + return stats + + def download_img(img, html_dir, ace_dir, chg=1): """ Download image from web site. Save image in {args.html_dir}/ACE/Plots From 70b5c3bf660f344372124d03cec4e02d0be3979f Mon Sep 17 00:00:00 2001 From: Malgosia Sobolewska Date: Tue, 1 Feb 2022 17:19:57 -0500 Subject: [PATCH 3/4] Used median instead of mean, ensured that spectral index is not negative --- ACE/Scripts/create_ace_html_page.py | 32 ++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/ACE/Scripts/create_ace_html_page.py b/ACE/Scripts/create_ace_html_page.py index 9f19093..166a115 100755 --- a/ACE/Scripts/create_ace_html_page.py +++ b/ACE/Scripts/create_ace_html_page.py @@ -84,7 +84,7 @@ def create_ace_html_page(): # Add table with stats for the display - for stat, name, fmt in zip((stats['means'], stats['mins'], stats['fluences_2h']), + for stat, name, fmt in zip((stats['medians'], stats['mins'], stats['fluences_2h']), ('AVERAGE', 'MINIMUM', 'FLUENCE'), ('{:11.3f}', '{:11.3f}', '{:11.4e}')): vals = stat.values() @@ -185,14 +185,18 @@ def read_2h_ace_data_from_12h_archive(fname): def calculate_stats(dat): """ - Calculate statistics: means, mins, 2h fluences in each channel, + Calculate statistics: medians, mins, 2h fluences in each channel, and spectral indexes. :param dat: astropy Table with ACE data :returns stats: dictionary with values being dictionaries of stats for ACE channels """ - stats = {'means': {}, 'mins': {}, 'fluences_2h': {}, 'spectra': {}} + stats = {'medians': {}, + 'mins': {}, + 'fluences_2h': {}, + 'spectra': {'p3_p5': -1e5, 'p3_p6': -1e5, + 'p5_p6': -1e5, 'p6_p7': -1e5}} cols = ("de1", "de4", "p2", "p3", "p3_scaled_p5", "p3_scaled_p5", "p5", "p6", "p7") @@ -203,19 +207,27 @@ def calculate_stats(dat): if len(dat[col][ok]) > 2: # At least 2 measurements of nominal quality # in the last 2 hours - stats['means'][col] = np.mean(dat[col][ok]) + stats['medians'][col] = np.median(dat[col][ok]) stats['mins'][col] = np.min(dat[col][ok]) - stats['fluences_2h'][col] = np.mean(dat[col][ok]) * 7200 + stats['fluences_2h'][col] = np.median(dat[col][ok]) * 7200 else: - stats['means'][col] = -1e5 + stats['medians'][col] = -1e5 stats['mins'][col] = -1e5 stats['fluences_2h'][col] = -1e5 # Spectral indexes - stats['spectra'] = {'p3_p5': stats['means']['p3'] / stats['means']['p5'], - 'p3_p6': stats['means']['p3'] / stats['means']['p6'], - 'p5_p6': stats['means']['p5'] / stats['means']['p6'], - 'p6_p7': stats['means']['p6'] / stats['means']['p7']} + meds = stats['medians'] + if all(x > 0 for x in [meds['p3'], meds['p5']]): + stats['spectra']['p3_p5'] = meds['p3'] / meds['p5'] + + if all(x > 0 for x in [meds['p3'], meds['p6']]): + stats['spectra']['p3_p6'] = meds['p3'] / meds['p6'] + + if all(x > 0 for x in [meds['p5'], meds['p6']]): + stats['spectra']['p5_p6'] = meds['p5'] / meds['p6'] + + if all(x > 0 for x in [meds['p6'], meds['p7']]): + stats['spectra']['p6_p7'] = meds['p6'] / meds['p7'] return stats From b21011e298730ee5245b751723b0d5757eca0925 Mon Sep 17 00:00:00 2001 From: Malgosia Sobolewska Date: Tue, 8 Feb 2022 16:35:18 -0500 Subject: [PATCH 4/4] Adjusted 2h filtering to account for processing time --- ACE/Scripts/create_ace_html_page.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ACE/Scripts/create_ace_html_page.py b/ACE/Scripts/create_ace_html_page.py index 166a115..7dc1cb1 100755 --- a/ACE/Scripts/create_ace_html_page.py +++ b/ACE/Scripts/create_ace_html_page.py @@ -168,13 +168,15 @@ def read_2h_ace_data_from_12h_archive(fname): # remove duplicate lines; TODO check why they appear dat = unique(dat, keys=("yr", "mo", "da", "hhmm")) - # keep only the most recent two hours of nominal quality data fits_fmt = [f"{yr}-{mo.zfill(2)}-{da}T{hhmm[:2]}:{hhmm[-2:]}:00.000" for yr, mo, da, hhmm in zip(dat['yr'], dat['mo'], dat['da'], dat['hhmm'])] tt = DateTime(np.array(fits_fmt), format='fits').date - now_minus_2h = DateTime(DateTime() - 2 / 24).date + + # keep only the most recent 2h of nominal quality data + # subtract 2.25h to account for cronjob processing every 5 min + now_minus_2h = DateTime(DateTime() - 2.25 / 24).date # status: 0 = nominal, 4,6,7,8 = bad data, unable to process, 9 = no data ok = (dat['e_status'] == 0) & (dat['p_status'] == 0) & (tt > now_minus_2h)