How to find a unit vector in teh QD northward direction (along QD meridians) at any location, using apexpy:
import apexpy
import numpy as np
epoch = 2015 # year at which IGRF coefficients are evaluated
a = apexpy.Apex(epoch)
glat, glon = 60, 30 # set geodetic latitude and longitude
f1, f2 = a.basevectors_qd(glat, glon, 0, coords = 'geo') # last param is height - set to ground here
$\bf{f}_2$ is a vector that points along contours of constant longitude - i.e., along magnetic meridians. The apexpy function returns the east and north-components of this vector. However, it is not a unit vector, so the last step to find a QD northward unit vector is normalize $\bf{f}_2$:
qd_north = f2 / np.linalg.norm(f2)
print('QD north at (glat, glon) = (%.1f, %.1f) is %.3f east, %.3f north' % (glat, glon, qd_north[0], qd_north[1]))
QD north at (glat, glon) = (60.0, 30.0) is -0.177 east, 0.984 north