How can I make Mathematica solve this set of four equations?

The matrix logarithm in Alex's answer will give one out of many possible (complex!) solutions, in complete analogy with the scalar case.

One way to go about this is to simultaneously reduce cmat and the matrix within the exponential to the Jordan form:

{sm, jm} = JordanDecomposition[2 π {{a, b}, {c, d}}]
{sr, jr} = JordanDecomposition[{{Sin[1], Cos[1]}, {-Cos[1], Sin[1]}}] // FullSimplify

Conveniently, 1. both jm and jr are diagonal matrices; and 2. both sm and sr are normalized such that their second row is one. We then recall that the Jordan vectors of $\mathbf A$ and $\exp(\mathbf A)$ should be the same, so:

GroebnerBasis[Thread[First[sm] == First[sr]], {a, b, c, d}]
   {b + c, a - d, -2 I c + Sqrt[a^2 + 4 b c - 2 a d + d^2]}

Immediately, we find that $c=-b$ and $a=d$. We can use this to simplify the next set of equations:

eq = Simplify[TrigToExp[Thread[Diagonal[jr] == Exp[Diagonal[jm]]] /.
                        {c -> -b, d -> a}], b < 0]
   {-I + E^(I + 2 a π + 2 I b π) == 0, 
    E^(2 (a - I b) π) == -I E^I}

Feeding this to Solve[] (Solve[%, {a, b}] // FullSimplify) and then plugging the results into the original matrix will generate a set of parametrized solutions:

{{I u, 1/4 - 1/(2 π) + v},
 {-((-2 + π + 4 π v)/(4 π)), I u}}

and

{{1/2 I (1 + 2 u), -((2 + π)/(4 π)) + v},
 {(2 + π - 4 π v)/(4 π), 1/2 I (1 + 2 u)}}

where I have replaced the C[k] with simpler parameters for clarity. Here, u and v are integers. In particular, Alex's solution corresponds to the first set, with u = 0 and v = 0.


There is a solution with using MatrixLog

ClearAll[a, b, c, d];
T = 2 Pi; m = T {{a, b}, {c, d}}; q = 1;
cmat = {{Sin[q], Cos[q]}, {-Cos[q], Sin[q]}};


 NSolve[m == MatrixLog[cmat], {a, b, c, d}]

Out[]= {{a -> -1.76697*10^-17 + 0. I, b -> 0.0908451 + 0. I, 
  c -> -0.0908451 + 0. I, d -> 1.76697*10^-17 + 0. I}}

I have no idea why, but this seems to work. Starting with your code

T = 2 Pi;
bt = MatrixExp[{{a, b}, {c, d}}*T];
cmat = {{Sin[1], Cos[1]}, {-Cos[1], Sin[1]}};
eqs = {cmat[[1, 1]] == bt[[1, 1]], cmat[[1, 2]] == bt[[1, 2]], 
   cmat[[2, 1]] == bt[[2, 1]], cmat[[2, 2]] == bt[[2, 2]]};

and "eliminating" one of the variables

elim = FullSimplify @ Eliminate[eqs, a]

gives a new set of equations. I put "eliminate" in quotes because it doesn't actually do that, elim still has the variable a in it. However, plugging this to Solve

Solve[elim, {a, b, c, d}]

gives a bunch of solutions, one of which is the numerical solution you posted.

Solutions

Using Resolve instead, like

FullSimplify @ Reduce[elim, {a, b, c, d}]

gives a parametrized solution with $a=d$ and $b=-c$, like in J. M.'s elegant answer.

Parametrized solution set