Equation for non-orthogonal projection of a point onto two vectors representing the isometric axis?

Let's call your two "axis" vectors $U$ and $V$. The nice folks around here would probably call them "basis" vectors. And suppose you have some other vector $W$, lying in the plane of $U$ and $V$, that you want to decompose.

Essentially, you want to find two numbers $h$ and $k$ such that $$W = hU + kV$$ The vector $hU$ will be the projection of $W$ onto the $U$-axis in a direction parallel to $V$, and $kV$ will be the projection of $W$ onto the $V$-axis in a direction parallel to $U$. Draw a picture and you'll see why this is so.

Take dot products of the equation above with $U$ and $V$ in turn, getting $$h(U \cdot U) + k(U \cdot V) = W \cdot U$$ $$h(V \cdot U) + k(V \cdot V) = W \cdot V$$ Now you have two simultaneous equations for $h$ and $k$, which are easy to solve (using Cramer's rule, for example).

Of course, all of this is much easier if $U$ and $V$ are perpendicular, because then
$U \cdot V$ is zero. If $U$ and $V$ have unit length, it gets even easier, because then $U \cdot U$ and $V \cdot V$ are both equal to 1.

This projection onto two non-orthogonal vectors can be calculated more simply using the cross product. We know that the cross product of a vector with itself results in the zero vector. This is used to filter out both coefficients.

Suppose we have three vectors $\textbf{v}_1$, $\textbf{v}_2$ and $\textbf{v}_3$ all lying in the plane defined by the perpendicular unit vector $\textbf{n}_z$. Now we wish to find the projection of vector $\textbf{v}_1$ onto $\textbf{v}_2$ and $\textbf{v}_3$ such that:

$$\textbf{v}_1=α\textbf{v}_2+β\textbf{v}_3$$

Using the cross product of both base vectors we get two equations:

$$\textbf{v}_2×\textbf{v}_1=α\textbf{v}_2×\textbf{v}_2+β\textbf{v}_2×\textbf{v}_3$$ $$\textbf{v}_3×\textbf{v}_1=α\textbf{v}_3×\textbf{v}_2+β\textbf{v}_3×\textbf{v}_3 $$ As the cross product with itself maps to zero we can write: $$\textbf{v}_2×\textbf{v}_1=β\textbf{v}_2×\textbf{v}_3$$ $$\textbf{v}_3×\textbf{v}_1=α\textbf{v}_3×\textbf{v}_2 $$ Now the coefficients can be extracted for example using the fact the cross products are perpendicular to the plane defined by $\textbf{n}_z$. The dot product of the normal vector of the plane maps to the scalar: $$\textbf{n}_z\cdot (\textbf{v}_2×\textbf{v}_1 )=β\textbf{n}_z\cdot (\textbf{v}_2×\textbf{v}_3)$$ $$\textbf{n}_z\cdot (\textbf{v}_3×\textbf{v}_1)=α\textbf{n}_z\cdot (\textbf{v}_3×\textbf{v}_2) $$ Resulting in the following system: $$\alpha = \frac{\textbf{n}_z\cdot (\textbf{v}_3×\textbf{v}_1)}{\textbf{n}_z \cdot (\textbf{v}_3×\textbf{v}_2)}$$ $$ \beta = \frac{\textbf{n}_z\cdot (\textbf{v}_2×\textbf{v}_1)}{\textbf{n}_z \cdot (\textbf{v}_2×\textbf{v}_3)}$$

The nice thing is that this also work for non-unit length vectors and we need no matrix inversions.

Just directly building upon bubba's answer, in Matrix form this can be computed by doing $$Y = \left(M^TM\right)^{-1}M^TW$$ where $M=\left\{\begin{matrix} U \;\Big\vert\; V \end{matrix}\right\}$ is the matrix that contains the basis vectors and $Y=\left\{\begin{matrix} h\\ k \end{matrix}\right\}$ is the projection of $W$ (your original vector) on said basis.