Skip to content

Instantly share code, notes, and snippets.

@arodland
Created January 30, 2020 02:57
Show Gist options
  • Save arodland/52c70c706c279f67ae5b750aef988039 to your computer and use it in GitHub Desktop.
Save arodland/52c70c706c279f67ae5b750aef988039 to your computer and use it in GitHub Desktop.
from datetime import datetime, timezone, timedelta
import igrf12
import json
import numpy as np
import os
import re
from scipy.optimize import minimize_scalar
import subprocess
import urllib.request
import sys
import multiprocessing
now = datetime.utcnow()
data = None
def station_weight(station, tm):
hour = float(tm.hour) + float(tm.minute) / 60. + float(tm.second) / 3600.
local_time = (hour + float(station['longitude']) / 15.) % 24.
abs_dip = abs(station['dip_angle'])
time_weight = np.sqrt(2 + np.cos((local_time - 14.5) * np.pi / 12.)) - 0.73
if abs_dip < 25.:
return 0.75 * time_weight
elif abs_dip < 50.:
return 1. * time_weight
elif abs_dip < 56.:
return 0.75 * time_weight
else:
return 0.5 * time_weight
def recency_weight(tm):
delta = now - tm
return np.power(2.0, -(delta / timedelta(hours=6)))
def cs_to_stdev(cs):
if cs == -1 or cs > 100:
cs = 69
if cs == 100:
cs = 81
return 0.191 - 0.00147 * cs
iri = None
pred_pool = None
def get_pred(station, ssn, tm):
global iri
if iri is None:
iri = subprocess.Popen(
['/src/build/iri_opt'],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
universal_newlines=True,
)
iri.stdin.write("%f %f %f %s\n" % (station['latitude'], station['longitude'], ssn, tm.strftime("%Y %m %d %H %M %S")))
iri.stdin.flush()
data = [float(x) for x in iri.stdout.readline().split()]
return (data[3], data[6], data[5]) # fof2, hmf2, mufd
def station_err(x):
station, ssn = x
station_total = 0.0
station_tw = 0.0
num_records = 0
last_rw = 0.0
if len(station['history']) < 3:
return (0,0)
last_tm = station['history'][-1][0]
if (now - last_tm).total_seconds() > 7200:
return (0,0)
for record in station['history']:
tm, cs, fof2, mufd, hmf2 = record
if cs <= 25:
continue
fof2_pred, hmf2_pred, mufd_pred = get_pred(station, ssn, tm)
stdev = cs_to_stdev(cs)
err = 0.6 * np.square(np.log(fof2_pred) - np.log(fof2))
err += 0.3 * np.square(np.log(mufd_pred) - np.log(mufd))
err += 0.1 * np.square(np.log(hmf2_pred) - np.log(hmf2))
err = np.sqrt(err) / stdev
sw = station_weight(station, tm)
rw = recency_weight(tm)
weight = sw * rw
last_rw = rw
station_total += err * weight
station_tw += weight
num_records += 1
if num_records >= 3:
return (last_rw * station_total / float(num_records), last_rw * station_tw / float(num_records))
else:
return (0,0)
def err(ssn):
total = 0.0
total_weight = 0.0
# print(ssn, file=sys.stderr)
results = pred_pool.imap_unordered(station_err, [ (station, ssn) for station in data ])
for result in results:
total += result[0]
total_weight += result[1]
return total / total_weight
if __name__ == '__main__':
with urllib.request.urlopen(os.getenv("HISTORY_URI")) as res:
data = json.loads(res.read().decode())
for station in data:
mag = igrf12.igrf(
'2019-12-31',
glat=station['latitude'],
glon=station['longitude'],
alt_km=1.
)
station['dip_angle'] = mag['incl'].item()
for record in station['history']:
record[0] = re.sub(r"\.\d+$", "", record[0])
record[0] = datetime.strptime(record[0], '%Y-%m-%d %H:%M:%S')
pred_pool = multiprocessing.Pool(4)
res = minimize_scalar(err, bounds=(-20.0, 200.0), method='Bounded', options={'xatol':0.1})
ssn = res.x
sfi = 63.75 + ssn * (0.725 + ssn*0.000089)
print("%s\t%f\t%f\t%f\t%f" % (now.strftime("%s"), sfi, ssn, res.fun, err(-100.0)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment