The candidates include:
There are problems in using these as more than local charts, to do with their multiple-valued nature, and singularities. That is, one must be careful above all to work only with diffeomorphisms in the definition of chart. Problems of this sort are inevitable, since SO(3) is diffeomorphic to real projective space RP3, which is a quotient of S3 by identifying antipodal points, and charts try to model a manifold using R3.
This explains why, for example, the Euler angles appear to give a variable in the 3-torus, and the unit quaternions in a 3-sphere. The uniqueness of the representation by Euler angles breaks down at some points (cf. gimbal lock), while the quaternion representation is always a double cover, with q and −q giving the same rotation.
If we use a skew-symmetric matrix, every 3×3 skew-symmetric matrix is determined by 3 parameters, and so at first glance, the parameter space is R3. Exponentiating such a matrix results in an orthogonal 3×3 matrix of determinant 1--in other words, a rotation matrix, but this is a many-to-one map. It is possible to restrict these matrices to a ball around the origin in R3 so that rotations do not exceed 180 degrees, and this will be one-to-one, except for rotations by 180 degrees, which correspond to the boundary S2, and these identify antipodal points. The 3-ball with this identification of the boundary is RP3. A similar situation holds for applying a Cayley transform to the skew-symmetric matrix.
Axis angle gives parameters in S2×S1; if we replace the unit vector by the actual axis of rotation, so that n and −n give the same axis line, the set of axis becomes RP2, the real projective plane. But since rotations around n and −n are parameterized by opposite values of θ, the result is an S1 bundle over RP2, which turns out to be RP3.
Fractional linear transformations use four complex parameters, a, b, c, and d, with the condition that ad-bc is non-zero. Since multiplying all four parameters by the same complex number does not change the parameter, we can insist that ad-bc=1. This suggests writing (a,b,c,d) as a 2×2 complex matrix of determinant 1, that is, as an element of the special linear group SL(2,C). But not all such matrices produce rotations: conformal maps on S2 are also included. To only get rotations we insist that d is the complex conjugate of a, and c is the negative of the complex conjugate of b. Then we have two complex numbers, a and b, subject to |a|2+|b|2=1. If we write a+bj, this is a quaternion of unit length.
Ultimately, since R3 is not RP3, there will be a problem with each of these approaches. In some cases, we need to remember that certain parameter values result in the same rotation, and to remove this issue, boundaries must be set up, but then a path through this region in R3 must then suddenly jump to a different region when it crosses a boundary. Gimbal lock is a problem when the derivative of the map is not full rank, which occurs with Euler angles and Tait-Bryan angles, but not for the other choices. The quaternion representation has none of these problems (being a two-to-one mapping everywhere), but it has 4 parameters with a condition (unit length), which sometimes makes it harder to see the three degrees of freedom available.
One area in which these considerations, in some form, become inevitable, is the kinematics of a rigid body. One can take as definition the idea of a curve in the Euclidean group E(3) of three-dimensional Euclidean space, starting at the identity (initial position). The translation subgroup T of E(3) is a normal subgroup, with quotient SO(3) if we look at the subgroup E+(3) of direct isometries only (which is reasonable in kinematics). The translational part can be decoupled from the rotational part in standard Newtonian kinematics by considering the motion of the center of mass, and rotations of the rigid body about the center of mass. Therefore any rigid body movement leads directly to SO(3), when we factor out the translational part.