This answer doesn't talk about how to do it Python at all: rather how to deal with the rotation. I think once you can do that then turning the maths into Python is simple. Initially I'll assume that you are computing positions in terms of three orthogonal axes, and the positions look like $(x, y, z)$, and you're just projecting these down onto the $(x, y)$ plane (so, no perspective). So the question is given a set of rotated axes, $(x', y', z')$, how do you convert from one set to the other. Once you have $x' = x'(x, y, z)$ and similarly for the rest, you can compute the coordinates in the $(x', y', z')$ coordinates & then project down onto the $(x', y')$ plane.
The way to do this is to define rotation matrices: there are three angles you need to know. To see why, consider the new axes: there are two angles which define where the $z'$ axis is, and then you can spin the whole coordinate system around that axis, which is another angle.
And now I am going to make a mess of this, because I always get confused between rotation of vectors and rotation of coordinate systems: there are sign differences.
The rotation matrices are simply 3 by 3 matrices $R$ such that $R^TR = RR^T = I$ and $\det{R} = 1$, where $R^T$ means the transpose of $R$. These in fact are elements of a representation of the special orthogonal group in 3 dimensions, $SO(3)$: it's worth looking that up.
We can define these things as products of three basic rotations, about the $x, y, z$ axes respectively:
$$
\begin{align}
R_x(\theta) &=
\begin{bmatrix}
1 & 0 & 0\\
0 & \cos\theta & -\sin\theta\\
0 & \sin\theta & \cos\theta
\end{bmatrix}\\
R_y(\theta) &=
\begin{bmatrix}
\cos\theta & 0 & \sin\theta\\
0 & 1 & 0\\
-\sin\theta & 0 & \cos\theta\\
\end{bmatrix}\\
R_z(\theta) &=
\begin{bmatrix}
\cos\theta & -\sin\theta & 0\\
\sin\theta & \cos\theta & 0\\
0 & 0 & 1
\end{bmatrix}
\end{align}
$$
Then a general rotation about three angles $\alpha, \beta, \gamma$ is
$$
R(\alpha,\beta,\gamma) = R_x(\alpha)R_y(\beta)R_z(\gamma)
$$
Note that the multiplications here are of course matrix multiplication: the thing that in Numpy is np.matmul, and in particular they're not element-wise multiplication. Then, finally, to compute the new coordinates you do
$$[x',y',z'] = R(\alpha,\beta,\gamma)[x,y,z]^T$$
Where, again, this is matrix multiplication of course.
So here's an example: if we start off with a point at $(x, 0, 0)$, then what are its coordinates in a set of axes rotated by $\theta$ about $z$? Well, you can do the multiplications (or get your tame algebra system to do them for you) and the answer is $(x', y', z') = (x\cos\theta, x\sin\theta,0)$, and it's clear that the rotation of the axes is $\theta$, clockwise. I think this means I've botched a sign somewhere, but it doesn't really matter. And after a rotation about the $x$ axis you'll get $(x',y',z') = (x, 0, 0)$ which is obviously geometrically correct.
Combining rotations becomes less intuitive I think, but the maths just works.
I find by far the best way to work out what's going on is to write a program which implements the transformations, and then take some plots you understand and transform them, which will fairly quickly show you what the various angles mean and what mistakes you've made.
