I will show how to calculate magnetic apex (QD) coordinates, QD vector compontents, and magnetic local time. I will use apexpy for the coordinates and pyamps for MLT, since the MLT code in pyamps is faster than in apexpy.
Start by defining some dummy coordinates and measurements:
import numpy as np
glat, glon = np.array([70, 30, -60, 80]), np.zeros(4) # geodetic latitude and longitude
B = np.random.random((len(glat), 10000, 2)) # 10000 measurements, 2 components, from each station
set up Apex object - epoch 2010, and calculate coordinates:
import apexpy
a = apexpy.Apex(2010)
qdlat, qdlon = a.geo2qd(glat, glon, 0) # last parameter is height
# get QD base vectors:
f1, f2 = a.basevectors_qd(glat, glon, 0, coords = 'geo')
# calculate F parameters:
F = f1[0]*f2[1] - f1[1]*f2[0]
# calcualte QD components:
Eqd = f1[0].reshape((-1, 1)) * B[:, :, 0] / F.reshape((-1, 1)) + f1[1].reshape((-1, 1)) * B[:, :, 1] / F.reshape((-1, 1))
Nqd = f2[0].reshape((-1, 1)) * B[:, :, 0] / F.reshape((-1, 1)) + f2[1].reshape((-1, 1)) * B[:, :, 1] / F.reshape((-1, 1))
Calcualte magnetic local time using pyamps:
from pyamps.mlt_utils import mlon_to_mlt
import pandas as pd # need to set up some dummy times:
t = pd.date_range(start = '2010-01-01 00:00', end = '2012-01-01 00:00', periods = B.shape[1])
# loop through each location and calculate the MLT for each time (same length as )
mlt = np.array([mlon_to_mlt(lon, t, epoch = 2010) for lon in qdlon])
done!