xpatch/plot_midas_vipir_hf.py at master · alexchartier/xpatch · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
import numpy as np
import pdb
import scipy.io
import datetime as dt
import h5py
import pickle
from plt_hf_vipir_nmf2 import load_vipir_nmf2(starttime, endtime, step, in_fname_fmt):


def main(
    tec_fname = 'data/midas_profs.mat'
    hf_fname_fmt = 'data/prc_analysis/no_badfrq/daily/data/%Y%b%d_analysis.pkl'
    vipir_fname_fmt = '/Users/chartat1/xpatch/data/vipir/Result/result_%Y_%m_%d.nc',
    start_time = dt.datetime(2019, 3, 1)
    end_time = dt.datetime(2019, 3, 13)
):
    """
    Plot NmF2 from MIDAS against RadioICE and VIPIR
    (need to make two different midas_profs.mat to cover both locations)
    """


    vipir_out = load_vipir_nmf2(starttime, endtime, step, in_fname_fmt)
    # Set up storage
    bins = []
    time = start_time
    while time < end_time:
        bins.append(time)
        time += dt.timedelta(minutes=15)
    bins = np.array(bins)
    tec_cts = np.zeros(bins.shape)
    hf_cts = np.zeros(bins.shape)

    # Load MIDAS
    Midas = {}
    with h5py.File(tec_fname, 'r') as f:
        for k, v in f.items():
            Midas[k] = v[...]
    Midas['dt'] = []
    for dn in Midas['Time']:
        Midas['dt'].append(datenum_to_datetime(dn))
    Midas['dt'] = np.array(Midas['dt'])
    tind = np.logical_and(Midas['dt'] >= start_time, Midas['dt'] <= end_time)
    for k, v in Midas.items():
        Midas[k] = v[tind]

    TEC_cutoff = 3.0
    for ind, t in enumerate(Midas['dt']):
        if Midas['TEC'][ind] > TEC_cutoff:
            delta_t = [td.total_seconds() for td in (bins - t)]
            tecind = np.argmin(np.abs(delta_t))
            tec_cts[tecind] = 1
            assert delta_t[tecind] < 1, 'Should be within 1 sec of a bin'

    time = start_time
    while time < end_time:
        # HF counts
        with open(time.strftime(hf_fname_fmt), 'rb') as f:
            HF = pickle.load(f)
        for t in HF[7.2]['time']:
            delta_t = [td.total_seconds() for td in (bins - t)]
            hfind = np.argmin(np.abs(delta_t))
            assert delta_t[hfind] < 15 * 60, 'Should be within 15 minutes of a bin'
            hf_cts[hfind] = 1

        time += dt.timedelta(days=1)

    notec_cts = np.logical_not(tec_cts)
    nohf_cts = np.logical_not(hf_cts)
    TEC_HF = np.float(np.sum(np.logical_and(tec_cts, hf_cts)))
    TEC_noHF = np.float(np.sum(np.logical_and(tec_cts, nohf_cts)))
    noTEC_HF = np.float(np.sum(np.logical_and(notec_cts, hf_cts)))
    noTEC_noHF = np.sum(np.logical_and(notec_cts, nohf_cts))
    print('         TEC > %1.1f,   TEC < %1.1f' % (TEC_cutoff, TEC_cutoff))
    print('7.2 MHz    %i          %i'  % (TEC_HF, noTEC_HF))
    print('No 7.2     %i          %i'  % (TEC_noHF, noTEC_noHF))
    print('True positive prediction: %2.0f pct' % (100. * TEC_HF / (TEC_HF + TEC_noHF)))
    print('True negative prediction: %2.0f pct' % (100. * noTEC_noHF / (noTEC_noHF + noTEC_HF)))


def datenum_to_datetime(datenum):
    """
    Convert Matlab datenum into Python datetime.

    :param datenum: Date in datenum format
    :return:        Datetime object corresponding to datenum.
    """
    days = datenum % 1
    hours = days % 1 * 24
    minutes = hours % 1 * 60
    seconds = minutes % 1 * 60
    return dt.datetime.fromordinal(int(datenum)) \
           + dt.timedelta(days=int(days)) \
           + dt.timedelta(hours=int(hours)) \
           + dt.timedelta(minutes=int(minutes)) \
           + dt.timedelta(seconds=round(seconds)) \
           - dt.timedelta(days=366)

if __name__ == '__main__':
    main()