Say I am given a point in an x1,y1,z1 coordinate system. I have a different coordinate system, x2,y2,z2 that has the same origin as the x1,y1,z1 system, but the axis are not aligned. I have roll, pitch, and yaw sensors fitted on the x2,y2,z2 coordinate system.
This picture may help:
How do I go about transforming the point in x1,y1,z1 coordinates to x2,y2,z2 coordinates? My original thought was that I simply just rotate by "head" degrees about y1 by the heading first, then rotate by "roll" degrees about z, then rotate by "pitch" degrees about x. This seems too easy and direct. Is this the correct approach?

Choose a global reference frame X, Y, Z. Then choose a local reference frame x, y, z that is attached to the rotating body (e.g. the person's head). Align you frames of reference, i.e. align x with X, y with Y, z with Z with the same origin.
One can represent the xyz coordinate system with unit vectors $\hat{x}, \hat{y}, \hat{z}$.
Now the task in hand is to find the rotation operation $R$ that correctly transforms $\hat{x}_1, \hat{y}_1, \hat{z}_1$ to $\hat{x}_2, \hat{y}_2, \hat{z}_2$, i.e. from the initial state to the final state of the local reference frame. This rotation operator has two purposes; $R$ represents the state of the local coordinate system with respect to the global one and by doing so it permits you to transforms all vectors between the final local reference frame and the initial/global reference frames (not just the unit vectors).
If done by hand, you can combine three Euler rotations to obtain $R$. But try using this Python library (
from scipy.spatial.transform import Rotation) which allows you to use and to convert between different rotation representations such as quaternions (to avoid coordinate singularities), Euler angles, rotation matrixes, and much more.