There exists a remarkably simple explanation for how electromagnetic phenomena must be transformed in a moving frame if: (1) the physics of that frame is to remain invariant; and (2) the speed of light c must remain invariant for both the at-rest observer and within the moving frame. These are Einstein's original two SR postulates, of course. (There is actually a third assumption needed about dimensional scaling orthogonal to the direction of motion. Lorentz discussed that issue explicitly, but Einstein seems to assumed as a given.)
The technique is to use light pendulums to define proper time and length in both frames, then show geometrically how this affects electromagnetics phenomena. I've only used it myself to bring out the energy relations, which is part of Einstein's subsequent very short paper that led to $E=mc^2$ (what a delightfully strange technique that paper uses!). All of the Maxwell transformations are necessarily implicit those energy transformations of moving light pendulums. They have to be, else SR would not work.
I have some light pendulum diagrams on hand that I had intended for use in time dilation discussions. I'll append them to this answer shortly, and explain at least briefly how they can be used to re-interpret electromagnetics and Maxwell's equations. I don't have time for a full treatment (flying tomorrow), but the light pendulum perspective is both helpful and deeply linked to Einstein's very first (and also much later) explanations of why special relativity is an unavoidable consequence of holding c invariant in all frames of motion.
So, the promised Graphical Overkill ensues below. I don't even have time tonight to explain the graphics, but I to promise I'll get back to them.
The initial foray into EM implications is, alas, the very last bullet in the very last slide. But it does show how transformations of both wavelength and energy are inherent and unavoidable within the curious constraints of having to preserve the physics and speed of light in both the frame doing the observing, and in the frame being observed.
This final graphic is at the very heart of how Einstein first came up with the theory of special relativity. He postulated two things: (1) that the speed of light remains invariant in all unaccelerated frames of reference, and (2) that all physics, including not just mechanics (Galilean relativity) but also electromagnetics, must remain unchanged regardless of frame.
Oh my... what those two postulates do in combination!
Light pendulums give a nicely succinct way of exploring all of those implications using a single device. The rho clocks I use here are simply light pendulums with enough design constraints and frame-specific labels attached to allow those implications to be explored.
So, take a look at this final figure:
Notice that the sides of the gray rectangles represent light paths drawn out in spacetime. So, if the speed of light is invariant, guess what? An object that is moving cannot have the shortest possible inbound an outbound light paths, because the light has to constantly pace the object to keep up with it. If the object moves very, very fast, the time it takes for the outbound light to catch up with its leading edge can become very long indeed, as viewed by an observer the "rest" frame labeled 0 ($\phi_0$).
How long? Well, look at the figure: It's just the height of the scaled gray box, which is what the rho clock ($\rho_1$) for frame 1 ($\phi_1$) looks like when viewed from the rest frame $\phi_0$. For shorthand, I call that situation $\rho_{1:>0}$, where "1:>0" just means that frame 1 is being observed by frame 0 (the ":>" is two eyes looking left). And yes, the label on that height is the traditional $\gamma$ of special relativity. The idea of $\gamma$ just emerges a lot more naturally (and a lot more geometrically) in rho diagrams.
As for electromagnetics, recall that both frames must see their physics unchanged. But notice how the pulses of blue light on the left side of the gray rectangle get crowded together to make that true. That doesn't happen just for pulse spacing, it happens for all of electromagnetic theory. So, for example, light traveling along the left branch of the frame 1 rho rectangle ("rho" actually stands for rectangle, not relativity) must increase in frequency to keep the internal physics of frame 1 rho clock unchanged (shorthand: $\rho_{1:>1}$ must remain invariant).
That has huge energy implications over in frame 0, however, since the faster the object travels, the more energetic that side of the rho clock rectangle becomes. If you combine that with some very unusual arguments from Einstein in his second special relativity paper (I like to call it his "asymptotic tautology argument"), you wind up eventually with the famous equation $E=mc^2$. (Side comment: That equation is more famous, but $E^2 = (pc)^2 + (mc^2)^2$ is a lot more useful; just ask any particle physicist.)
The nice thing about Maxwell's equations is that they already accommodated and allowed such unusual forms of scaling long before Einstein and special relativity came around. That is one of the reason why you will see almost breathless praise for Maxwell from Einstein and other physicists involved in the early days of special relativity. The new theory helped them appreciate in new ways just how deep Maxwell's insights had been.
Best Answer
I'm no expert on the historical development of the subject, however I will offer a derivation.
Consider two frames of reference $S$ and $S'$, and suppose that $S'$ moves with speed $\textbf v$ with respect to $S$. Coordinates in $S$ and $S'$ are related by a Galileian transformation: $$\begin{cases} t' = t \\ \textbf x' = \textbf x-\textbf vt\end{cases}$$ To find how the fields transform, we note that a Lorentz transformation reduces to a Galileian transformation in the limit $c \to \infty$. In fact, under a Lorentz transformation the fields transform like: $$ \begin{cases} \textbf E' = \gamma (\textbf E + \textbf v \times \textbf B) - (\gamma-1) (\textbf E \cdot \hat{\textbf{v}}) \hat{\textbf{v}}\\ \textbf B' = \gamma \left(\textbf B - \frac{1}{c^2}\textbf v \times \textbf E \right) - (\gamma-1) (\textbf B \cdot \hat{\textbf{v}}) \hat{\textbf{v}}\\ \end{cases}$$ Taking the limit $c\to \infty$ so that $\gamma\to 1$, we obtain the Galileian transformations of the fields: $$ \begin{cases} \textbf E' = \textbf E + \textbf v \times \textbf B\\ \textbf B' = \textbf B\\ \end{cases}$$ We can then invert the transformation by sending $\textbf v \to -\textbf v$: $$ \begin{cases} \textbf E = \textbf E' - \textbf v \times \textbf B'\\ \textbf B = \textbf B'\\ \end{cases}$$ By the same reasoning, can obtain the Galileian transformation of the sources: $$ \begin{cases} \textbf J = \textbf J' + \rho' \textbf v\\ \rho = \rho'\\ \end{cases}$$ We know that the fields and sources satisfy Maxwell's equations in $S$: $$ \begin{cases} \nabla \cdot \textbf E = \rho/\epsilon_0\\ \nabla \cdot \textbf B = 0\\ \nabla \times \textbf E = -\frac{\partial \textbf B}{\partial t}\\ \nabla \times \textbf B = \mu_0 \left(\textbf J +\epsilon_0 \frac{\partial \textbf E}{\partial t} \right)\\ \end{cases}$$ Replacing the fields and sources in $S$ with those in $S'$ we obtain: $$ \begin{cases} \nabla \cdot \textbf (\textbf E' - \textbf v \times \textbf B') = \rho'/\epsilon_0\\ \nabla \cdot \textbf B' = 0\\ \nabla \times \textbf (\textbf E' - \textbf v \times \textbf B') = -\frac{\partial \textbf B'}{\partial t}\\ \nabla \times \textbf B' = \mu_0 \left(\textbf J' + \rho' \textbf v +\epsilon_0 \frac{\partial (\textbf E' - \textbf v \times \textbf B')}{\partial t} \right)\\ \end{cases}$$ As a last step, we need to replace derivatives in $S$ with derivatives in $S'$. We have: $$\begin{cases} \nabla = \nabla' \\ \frac{\partial }{\partial t} = \frac{\partial }{\partial t'} - \textbf v \cdot \nabla\end{cases}$$ Substituting and removing the primes and using vector calculus, we obtain: $$ \begin{cases} \nabla \cdot \textbf E + \textbf v \cdot (\nabla \times \textbf B) = \rho/\epsilon_0\\ \nabla \cdot \textbf B = 0\\ \nabla \times \textbf E = -\frac{\partial \textbf B}{\partial t}\\ \nabla \times \textbf B = \mu_0 \left(\textbf J + \rho \textbf v +\epsilon_0 \frac{\partial}{\partial t}( \textbf E - \textbf v \times \textbf B) - \epsilon_0 \textbf v \cdot \nabla (\textbf E - \textbf v \times \textbf B) \right)\\ \end{cases}$$
In a vacuum, we can take the curl of the fourth equation to obtain: $$c^2\nabla^2 \textbf B = \frac{\partial^2 \textbf B}{\partial t^2} + (\textbf v \cdot \nabla)^2 \textbf B - 2 \textbf v \cdot \nabla \left(\frac{\partial \textbf B}{\partial t}\right)$$ Substituting a wave solution of the form $\textbf B \sim \exp{i(\textbf k \cdot \textbf x -\omega t)}$ We obtain an equation for $\omega$, which we can solve to obtain: $$\omega = -\textbf v \cdot \textbf k \pm c |\textbf k|$$ Therefore the speed of propagation is the group velocity: $$\frac{\partial \omega}{\partial \textbf k} = -\textbf v \pm c \hat{\textbf{ k}}$$ which gives you the expected $c\pm v$ with an appropriate choice of $\textbf v$ and $\textbf k$.