You can compute these dimensions using modular symbols (an auxiliary space which has the same Hecke action as modular forms, but is easier to compute). Here's a Sage example for weight 4 cusp forms of level Gamma0(17):
sage: S=ModularSymbols(Gamma0(17), weight=4, sign=1).cuspidal_submodule()
sage: S.atkin_lehner_operator().charpoly().factor()
(x + 17) * (x - 17)^3
Sage normalises the Atkin–Lehner operator so that its square is multiplication by $(-N)^{k-2}$, since that normalisation works better for odd weights (it avoids introducing square roots). If you want to normalise so that the operator becomes an involution, as is usual for $\Gamma_0$ levels, you need to scale by $N^{(k - 2) / 2}$. So, in my example, the scaling factor is 17, and the output means that the +1 eigenspace has dimension 3, and the -1 eigenspace has dimension 1.
You can compute the local AL operators for each prime $p \mid N$ similarly, and get the simultaneous eigenspaces by taking intersections of kernels – here's a code snippet which does this:
sage: def AL_eigenspace_dims(N, k):
....: facs = [p^r for (p, r) in N.factor()]
....: S = ModularSymbols(Gamma0(N), sign=1, weight=k).cuspidal_submodule()
....: Wmats = [S.atkin_lehner_operator(f).matrix() for f in facs]
....: for s in cartesian_product([ [-1, 1] for f in facs ]):
....: Vs = [(Wmats[i] - s[i]*facs[i]**((k-2)/2)).kernel() for i in range(len(facs))]
....: V = reduce(lambda u,v: u.intersection(v), Vs)
....: print(s, V.dimension())
....:
sage: AL_eigenspace_dims(21, 2)
(-1, -1) 0
(-1, 1) 1
(1, -1) 0
(1, 1) 0
sage: AL_eigenspace_dims(21, 12)
(-1, -1) 7
(-1, 1) 6
(1, -1) 7
(1, 1) 8
This modular symbol method will be slow if k or N is large, because it involves linear algebra with matrices having about $kN$ rows and columns. (For $k N$ around 500, it takes roughly ten seconds on my machine). In contrast, the algorithm in Lloyd Kilford's draft code from your link uses the Riemann–Roch theorem, and the running time is dominated by factorising N, so it should be practical for N up to about 100 digits; but his implementation only covers the case of weight 2. In principle it should be possible to extend the Riemann–Roch method to higher weights, but I don't know if anyone has implemented this in Sage.