Is the seven-dimensional cross product unique?

My question has gathered several useful comments, which together have given me enough pieces of information that I think I now have the answer. I think this is correct, but please don't hesitate to correct me if I'm making any mistakes:

Throughout this answer, the letter $R$ will denote an element of (the fundamental representation of) $SO(7)$, and $\times_1$ and $\times_2$ will denote 7D cross products. We'll use square brackets to denote the rotated vector $R[a]$ and subscripts to denote any additional quantities that the rotation $R$ might depend on.

Since the definition of the cross product fixes its magnitude, we obviously have that $|a \times_1 b| \equiv |a \times_2 b|$ and so $$ a \times_2 b \equiv R_{\times_1, \times_2, a, b} [a \times_1 b]. $$ (I.e. the map is norm-preserving but not a priori linear in $a$ or $b$.)

However, it turns out that any two 7D cross products are related by a single rigid rotation of the entire vector space: that is, for any two cross products $\times_1$ and $\times_2$, we have $$ a \times_2 b \equiv \pm R_{\times_1, \times_2} \left[ R_{\times_1, \times_2}^{-1}[a] \times_1 R_{\times_1, \times_2}^{-1}[b] \right]. $$ (We can drop the "$\pm$" if we allow $R$ to be a general element of $O(7)$ instead of $SO(7)$.) This nontrivial result is what Massey means when he says that the 7D cross product is "unique up to isomorphism"; by "isomorphism" he means a single rigid rotation of the entire vector space.

Now let's think about this from a different perspective. Above, we started with two arbitrary fixed cross products $\times_1$ and $\times_2$ and considered a rotation $R_{\times_1, \times_2}$ that related them. Now let's instead start with a particular cross product $\times_1$ and imagine a generic "rotated version" $$ f_{\times_1,R}(a, b) := R\left[ R^{-1}[a] \times_1 R^{-1}[b] \right], \qquad R \in SO(7). $$ We can show that $f_{\times_1, R}$ is also a cross product $\times_2$ (with the same orientation as $\times_1$ if $R \in SO(7)$). This is in some sense the converse of the previous result; the previous result states that any two cross products are rotations of each other, and this result states that the rotation of any cross product is another cross product.

In 3D, it turns out that there is a single cross product $\times$ (up to orientation reversal) and $f_{R,\times} \equiv \times$ for all $R \in SO(3)$. But this is not true in 7D. In 7D it turns out that for any fixed cross product $\times_1$, the claim $f_{\times_1, R} \equiv \times_1$ is only true for a set of rotations $R$ that forms a proper subgroup of $SO(7)$ that is isomorphic to the exceptional Lie group $G_2$. (As Jason DeVito explains in the comments, the exact subgroup of $SO(7)$ depends on $\times_1$, but it is always isomorphic to $G_2$.)

The set of distinct 7D cross products is therefore $X := SO(7)/G_2$. (More precisely, this is the set of transformations that maps an arbitrary initial cross product into any other cross product.) Since $SO(n)$ is always a simple Lie group, $G_2$ is not a normal subgroup of $SO(7)$. So $X$ does not have the full structure of a Lie group, but only of a homogeneous space. In fact, it turns out that this homogeneous space $SO(7)/G_2$ is isometric to the real projective space $\mathbb{R}P^7$.

Now let's fix an orthonormal basis of $\mathbb{R}^7$ and consider the set of cross products that map basis vectors to basis vectors. One way to do this is to choose an arbitrary such cross product $\times_1$ and then consider the cosets $R\, G_2 \in X$ that transform $\times_1$ to another cross product with the same property of mapping basis vectors to basis vectors. We find that this is a finite subset of $X \cong \mathbb{R}P^7$ with 480 elements. (But the other elements of $X$ are still perfectly valid cross products.)

TLDR: There are an uncountably infinite number of distinct 7D cross products, and they form a Reimannian manifold isometric to the real projective space $\mathbb{R}P^7$.