Why intuitively do the quaternions satisfy the mixture of geometric and algebraic properties that they do? [I completely rewrote the question to see if I could make it clearer. The comments below won't make any sense. In fact, my original question has been answered by Eric Wolfsey, so I may restore it.]
When you read about the quaternions on Wikipedia and on many other sources, they are defined with the relation $$i^2 = j^2 = k^2 = ijk = -1$$
This comes across as completely random. There is no explanation of where this relation comes from.
These sources then go on to prove that the quaternions satisfy various algebraic and geometric properties. Among these properties, they act on 3D vectors as rotations (but as a double-covering), they act on 4D space as rotations, they're associative, they're distributive. Etc. All of this comes across as a random coincidence when you compare it to the defining relation.
Some of the time quaternions are thought of as 3D rotations. But this is a lossy interpretation, as a quaternion $q$ and its negation $-q$ represent the same rotation. For similar reasons, addition of quaternions starts to seem more bizarre.
When defined as 4D rotations, they happen to be an "isoclinic rotation". I haven't thought enough about this concept...
When defined as a complex matrix that satisfies a linear algebra relationship (sourced from here: https://qchu.wordpress.com/2011/02/12/su2-and-the-quaternions/), closure under multiplication is baked in, but closure under addition comes across as a coincidence. This relation is: $$M\text{ is a $2\times2$ matrix over $\mathbb C$}, M^\dagger M \in \mathbb R, \det(M) \geq 0 $$
Geometric Algebra textbooks show that they're an instance of a larger family of algebras with similar interpretations, called Geometric Algebra (or Clifford Algebras). But then why do Clifford Algebras exist? Why do they have their mixture of algebraic and geometric properties? It feels like you're replacing a small mystery with an even bigger mystery.
 A: The idea is that rotations have a certain rigidity in two dimensions: given a basis $\{e_1,e_2\}$, once you know where a rotation sends $e_1$, the image of $e_2$ is uniquely determined.  Moreover, the conditions that uniquely define the image of $e_2$ can be written linearly in the image of $e_1$ (if $\{e_1,e_2\}$ is an orthogonal basis, then the conditions are that the image of $e_2$ must be orthogonal to the image of $e_1$ and must make the matrix have determinant $1$).  As a result, scalar multiples of rotations form a linear subspace of all matrices.
Here's how it works in detail.  Let $V=\mathbb{R}^2$. Consider the bilinear map $\det:V\times V\to\mathbb{R}$ which takes a pair of vectors $(v,w)$ to the determinant of the matrix with columns $v$ and $w$.  We can think of this bilinear instead as a linear map $V\to V^*$, or, using the canonical isomorphism $V^*\cong V$ given by the inner product, as a linear map $J:V\to V$.  Explicitly, $J(v)$ is the unique vector which satisfies $\langle w,J(v)\rangle=\det(v,w)$ for all $w\in V$.  Even more explicitly, if $\{e_1,e_2\}$ are the standard basis for $V$, then $J(e_1)=e_2$ and $J(e_2)=-e_1$ (so $J$ is in fact just multiplication by $i$ when we identify $ae_1+be_2$ with $a+bi$).
The key observation is now that a matrix $T$ with columns $v$ and $w$ is a scalar multiple of an element of $SO(2)$ iff $w=J(v)$.  To prove this, we may scale $T$ to assume $v$ is a unit vector (if $v=0$ it is easy to see both conditions are equivalent to $w=0$).  Then $T$ is a scalar multiple of an element of $SO(2)$ iff it is an element of $SO(2)$ iff $w$ is orthogonal to $v$ and $\det(v,w)=1$.  Note that since the orthogonal complement of $v$ is $1$-dimensional, there is a unique $w$ orthogonal to $v$ such that $\det(v,w)=1$.  Now notice that $J(v)$ is orthogonal to $v$ since $\langle v,J(v)\rangle=\det(v,v)=0$ and that $\det(v,J(v))=\langle v,v\rangle=1$.  So in fact $J(v)$ is the unique $w$ that will make $T$ a scalar multiple of an element of $SO(2)$.  Thus $T$ is a scalar multiple of an element of $SO(2)$ iff $w=J(v)$.
Finally, the set of matrices which satisfy $w=J(v)$ is obviously closed under addition, and so we conclude that so is the set of scalar multiples of elements of $SO(2)$.  Moreover, we see that the map taking such a matrix to its first column is a linear isomorphism.  In this way we cover the standard identification between $\mathbb{C}$ (considered as the set of scalar multiples of elements of $SO(2)$) and $\mathbb{R}^2$.
[This discussion can be made more basis-free in various ways.  In particular, instead of explicitly talking about the columns of matrices, you can fix a unit vector in $\bigwedge^2 V$, use that unit vector to get an isomorphism $\bigwedge^2 V\cong\mathbb{R}$ and replace $\det$ by the wedge product map $V\times V\to \bigwedge^2 V\cong\mathbb{R}$ when defining $J$.  Then, one direction of the second part of the argument can be recast as a proof that if $T$ is a scalar multiple of an element of $SO(2)$ then $T$ commutes with $J$.  For the converse, I don't see a way to completely avoid using a basis.]
For quaternions and $SU(2)$, the story is the same, except that you take $V=\mathbb{C}^2$ and $J$ is conjugate-linear instead of linear (since the inner product is sesquilinear).
A: I will try to address this part of your question, as edited. (If you restore the original version, which I  have not seen, this may no longer be an answer, though it should stand as a piece of history.)

This comes across as completely random. There is no explanation of
  where this relation comes from.
These sources then go on to prove that the quaternions satisfy various
  algebraic and geometric properties. Among these properties, they act
  on 3D vectors as rotations (but as a double-covering), they act on 4D
  space as rotations, they're associative, they're distributive. Etc.
  All of this comes across as a random coincidence when you compare it
  to the defining relation.

Mathematicians are naturally curious. Each question they answer suggests others. Sometimes answers are useful, sometimes satisfying for their own sake.
Complex numbers were invented (or discovered, depending on your philosophical bent) in the 16th century to solve polynomial equations - mathematicians were surprised and puzzled to see that they were necessary in the formula for the roots of cubics even when all three roots were real. The wikipedia page is a good reference for this and the ensuing history. Extending calculus to complex functions of a complex variable  led to extraordinarily useful and interesting new mathematics. Moreover, complex numbers turned out to be a good way to study  rotations in the plane. 
Once mathematicians  defined multiplication of vectors in the plane (addition was always there) it was natural to wonder about how to define multiplication for vectors in space. The hope was to find an algebraic structure that could model rotations in space analogous to how the complex numbers model rotations in the plane. 
Many tried, and failed. Hamilton succeeded - in a way - by using a fourth dimension and inventing quaternions. The seemingly arbitrary $i^2 = j^2 = k^2 = -1$ and the anticommutation rules were forced on him so he could capture how rotations compose. That motivation is lost when you just start with a definition. In fact, matrix algebra (invented at about the same time) provided a better way to work with rotations algebraically. Quaternions languished until a recent resurgence as a reasonable data structure for manipulations in 3d computer graphics.
In the meanwhile mathematicians proved that  only in dimensions a power of $2$ was there even a hope for a reasonable vector multiplication -  and then you have to relax the meaning of "reasonable". Quaternion multiplication is noncommutative.  Multiplication in the Cayley algebra in dimension $8$ isn't associative.
This from the end of your question is a feature of mathematics, not a problem:

It feels like you're replacing a small mystery with an even bigger
  mystery.

