Inverting the equation for $T_{\mu\nu}$ in terms of $F_{\mu\nu}$ The Stress-Energy Tensor for electromagnetism is given by: 
$$ T_{\mu \nu} = F_{\mu}\,^{\alpha}F_{\nu\alpha}-\frac{1}{4}g_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta} $$
How can I find $F_{\mu\nu}$ in terms of $T_{\mu\nu}$?
Rewriting the above equation using:
$$ T_{\mu\nu}=- F_{\mu \alpha} g^{\alpha\beta} F_{\beta\nu} + \frac{1}{4} g_{\mu \nu}g^{\alpha\beta}F_{\beta\delta}g^{\delta\gamma} F_{\gamma\alpha}$$
from which we can write the following $4\times4$ matrix equation for the three matrices $T,\,F,\,g$, where $T$ is symmetric, $F$ is anti-symmetric and $g$ is symmetric and invertible:
$$ T = -F g^{-1} F+\frac{1}{4}\left(\mathrm{Tr}\, \left[g^{-1}Fg^{-1}F\right]\right)\,g$$
The only way I can think of is writing down 10 equations (as there are free components in $T^{\mu\nu}$) and then trying to find the 6 unknowns (as there are free components of $F^{\mu\nu}$).
Is there a better way to do this? 
 A: See Edit below, the original answer is not completely correct.
There is no gauge freedom in $F$. $F$ is gauge invariant.
In fact, $F$ is completely measurable. It's components are the Electric and Magnetic fields, so you just go out with a set of test charges and measure $E$ and $B$ and you've got $F$.
One hint that $T$ and $F$ do not contain the same amount of information is that they have different numbers of independent components. $F$ has 6 independent components as an antisymmetric tensor, while $T$ has 10 as a symmetric one. This isn't a proof of anything, but a hint that they are capturing different things. 
If you are working locally (ie, at a point), the simple way to see this explicitly is to use Lorentz transformations. The stress energy tensor has $10$ independent components since it is a symmetric tensor, we can use the $6$ Lorentz transformations to diagonalize $T$. Then we have 4 equations
\begin{eqnarray}
T_{00} &=& \frac{1}{2}\left(E^2 + B^2\right) \\
T_{ii} &=& (E_i^2 - \frac{1}{2}E^2) + (B_i^2 - \frac{1}{2}B^2)
\end{eqnarray}
There is no sum over $i$ implied in the second equation, it's just a quick way of writing the 3 spatial equations.
You can see that there is no way to solve these. For one thing, there are more components in $E$ and $B$ than there are in $T$ in this frame. For another, since the fields appear squared, there is no way to determine the sign of any of the components of $E$ or $B$. 
Additionally you can't tell the difference between $E$ and $B$ (ie, given $T_{00}$, who is to say whether you had $E^2=0$ or $B^2=0$ or neither)? This last point is a consequence of the electromagnetic duality: in the absence of matter, the physics of E/M is invariant under $E\rightarrow B$, $B\rightarrow -E$.
EDIT: 
The above is not quite correct in detail (though I think the conclusion is correct). For whatever reason I neglected the fact that there are always 10 components of $T_{\mu\nu}$, so there are always 10 equations, even in the frame in which $T$ is diagonal. In particular, there are also conditions like
\begin{eqnarray}
0 &=& E_x E_y + B_x B_y \\
0 &=& E_x B_y - E_y B_x
\end{eqnarray}
So my counting argument, "There are more variables than equations," was incorrect. This fits with the idea that $T$ has more components than $E$--if anything based on counting you would think that computing $T$ from $E$ was the harder thing to do. (In fact this is generically true--the stress energy tensors you get from field theory are not the most general stress energy tensors you can write down. There are plenty of stress energy tensors you can write down that won't come from a lagrangian).
The real reason this won't work, as far as I can tell, is the electromagnetic duality as well as the fact that everything is squared. There just isn't a way to distinguish $E$ from $B$ if you write out all the components. In other words, the duality means that the equations have a degeneracy, so there are fewer equations than it naively appears, so you can't solve for all the components.
On the other hand, if you know $T$ everywhere, not just locally, that is a totally different story. That's because (1) if you know $T$ everywhere you can differentiate it, and (2) $\partial_\mu T^{\mu\nu}=0$ is just maxwells equations $\partial_\mu T^{\mu\nu}=\partial_\mu F^{\mu\nu}$, possibly up to an overall factor. So then, up to the usual caveats about needing to know the boundary conditions, if you know $T$ everywhere you can solve maxwell's equations to obtain $F$.
Moral: don't believe everything you read on the internet. 
A: The easiest way I can think of in Minkowski space, short of doing the algebra in terms of matrices, is to use $$\begin{split} f^a &= \rho E^a + \epsilon^{abc} J_b B_c = \partial_b T^{ab} - \epsilon_0 \mu_0 \partial_t S^a\\ \frac{\partial T^{00}}{\partial t} &= - \vec{J}\cdot \vec{E}- \vec{\nabla} \cdot \vec{S} \end{split},$$ with $S^a \equiv \frac{1}{\mu_0} \epsilon^{abc} E_b B_c = T^{0a},$ $a,b \in \{1,2,3\}$, and hope that the fields will be easy to discern.
Perhaps the more symmetric form $$\frac{1}{2} (F_\mu{}^\rho F_{\nu\rho} + \star F_\mu{}^\rho \star\! F_{\nu\rho}),$$ with the star denoting the Hodge dual, will prove easier to handle in curved spacetime, if you manage to break $T_{\mu\nu}$ into a sum of matrices of similar structure.
