The answer is yes, since the exponential map $\exp: \mathbf{so}(3) \rightarrow \mathbf{SO}(3)$ is surjective (=onto).
Long answer:
Axis-angle can be represented using a $3$-vector $\omega$ while the magnitude $\theta=|\omega|$ is the rotation angle and $\mathbf{u}=^\omega/_\theta$ is the rotation axis. 3-Vectors are closed under the cross product:
$$\omega_1\in \mathbb{R}^3, \omega_2\in \mathbb{R}^3\Rightarrow (\omega_1\times \omega_2)\in\mathbb{R}^3.$$
Each such vector $\omega$ has an equivalent $3\times 3$ matrix representation
$\hat{\omega}$ (which is uniquely defined by $\hat{\omega}\cdot \mathbf{a} := \omega\times \mathbf{a}$ for $\mathbf{a}$ being a general 3-vector).
The space of matrices of the form $\hat{\omega}$ is called the Lie algebra $\mathbf{so}(3)$. Thus, one can show that matrices of the form $\hat{\omega}$ are closed under the Lie bracket $[A,B]=AB-BA$:
$$\hat{\omega_1}\in \mathbf{so}(3), \hat{\omega_2}\in \mathbf{so}(3)\Rightarrow [\hat{\omega}_1, \hat{\omega}_2]\in\mathbf{so}(3).$$
Now, let us consider the matrix exponential: $\exp(\mathtt{A}) = \sum_{i=0}^\infty \frac{\mathtt{A}^i}{i!} $. Two poperties can be shown:
(1) If $\hat{\omega}\in\mathbf{so}(3)$, then $\exp(\hat{\omega})\in\mathbf{SO}(3)$.
$\mathbf{SO}(3)$ is the special orthogonal group in three dimensions. Thus, it consists of matrices which are orthogonal ($\mathtt{R}\cdot \mathtt{R}^\top=\mathtt{I}$) and the determinant is 1. In other word, it is the group of pure rotations.
(2) The exponential map $\exp: \mathbf{so}(3) \rightarrow \mathbf{SO}(3)$ is surjective.
So, (1) says that every $\exp(\hat{\omega})$ is a rotation matrix. And, (2) says that for each rotation matrix $\mathtt{R}$, there is at least one axis-angle
representation $\omega$ so that $\exp(\hat{\omega})=\mathtt{R}$
Proofs of (1) and (2) are in corresponding text books, e.g. [Gallier, page 24].
The simplest way to model 3D rotations with quaternions is to view 3D space as the "imaginary quaternions", that is, use $i,j,k$ vectors like you do in physics to model 3D space. (But no real part.)
Every quaternion with norm 1 produces a rotation of this 3D space by conjugation $x\mapsto qxq^{-1}$, where $q$ is any length 1 quaternion.
The formula is surprisingly simple: pick a vector that points along the axis you want to rotate around, and scale it down to a unit vector, call it $h$. Then, take your angle of rotation $\theta$ and form $q=\cos(\theta/2)+h\sin(\theta/2)$. The result is a unit length quaternion that produces the rotation you wanted. ($\theta$ may be positive or negative, depending on the direction of rotation.)
I'm not sure I understand what you mean by "nullify rotation in the x and y axes", but I was guessing you meant you would like to be able to pick which axis you are rotating around.
Best Answer
Once again, the easiest way of doing this (IMHO, at least) turns out to be quaternions: a rotation of $\theta$ about the (normalized) vector $\hat{v} = (v_0, v_1, v_2)$ is represented by the quaternion ${\bf q}=\mathrm{cos}(\theta/2) + \mathrm{sin}(\theta/2)*(v_0{\bf i}+v_1{\bf j}+v_2{\bf k})$; and quaternion multiplication is easy to compute. All you need are the base axioms ${\bf i}^2 = {\bf j}^2 = {\bf k}^2 = {\bf i}{\bf j}{\bf k} = -1$; note that the latter implies for instance that ${\bf i}{\bf j} = {\bf k}$ by simply multiplying both sides on the right by ${\bf k}$ — on the right, because multiplication of quaternions, like composition of rotations, isn't commutative: $\bf{i}\bf{j}=\bf{k}$ but $\bf{j}\bf{i}=-\bf{k}$!. Once you have your resultant quaternion, the axis can be extracted as the (normalized) imaginary part of the result and the angle can be found by taking the arc-cosine of the scalar part. See the Wikipedia page on quaternions and rotations for more details on just how this works.