Is there a Hamiltonian for the (classical) electromagnetic field? If so, how can it be derived from the Lagrangian?

Yes. There is a standard way to generalize to field theory.

Let a theory of $n\geq 1$ fields $\phi^i$ with a Lagrangian density $\mathcal L = \mathcal L(\phi^i, \partial_\mu\phi^i)$ be given. Here we use that standard abuse of notation in which $\phi^i$ denotes the vector whose components are the fields; $\phi^i = (\phi^1, \dots, \phi^n)$.

To obtain the corresponding hamiltonian density, one first defines the following canonical momentum corresponding to the field $\phi^i$: \begin{align} \pi_i(x) = \frac{\partial \mathcal L}{\partial\dot \phi^i}(\phi^i(x), \partial_\mu\phi^i(x)), \qquad \dot\phi^i := \partial_t\phi^i \tag{1} \end{align} Then, the Hamiltian density is \begin{align} \mathcal H = \pi_i\dot\phi^i - \mathcal L \end{align} where a sum over $i$ is implied. Note that as in classical mechanics, on the right hand side of this expression, $\dot \phi^i$ should be replaced with its expression in terms of $\pi_i,\phi^i$ so that the Hamiltonian is a function of $(\pi_i, \phi^i)$ only, namely \begin{align} \mathcal H(\pi_i, \phi^i) = \pi_j \,\dot\phi^j(\pi_i,\phi^i) - \mathcal L(\phi^i, \dot\phi(\pi_i,\phi^i)). \end{align} Again we have abused notation slightly here in writing $\dot\phi^i$ as a function of $\pi_i$ and $\phi^i$. What we mean is the expression for $\dot\phi^i$ is obtained by solving the definition $(1)$ of the canonical momentum for $\dot\phi^i$ in terms of $\pi_i$ and $\phi^i$.

In your case, the fields are $A^\mu$ with corresponding momenta $\pi_\mu$.

Yes. Not only is it possible to find a Hamiltonian density but it is even possible to find a plain Hamiltonian.

I gave the construction in this EPL ArXiv1303.6143: Euro Physics Letters, 103 (2013) 28004

This Letter is short but dense, I can only outline the method here.

The reason why you can't write an hamiltonian for the electromagnetic field is that you represent it with a continuum of degrees of freedom. So you can only find a hamiltonian density. This is in constrast to "regular" mechanical systems which are represented by a discrete set of degrees of freedom. In this case you find a hamiltonian function not a density.

So the main step is therefore to discretize the field, but without committing any kind of approximation. This sounds contradictory, but actually it is not. Take the example of a periodic continuous function. It has an infinity of degrees of freedom. But since it is periodic you can represent it with a Fourier series. The Fourier coefficients are an infinite set of discrete parameters representing exactly your function. So you succeeded to discretize your function without committing any approximation. These coefficients will be the canonical variables of the hamiltonian.

In the case of the electromagnetic field we need a mathematical tool found by the russian mathematician Israel Gel'fand. This Gel'fand transformation is very similar to the Fourier series in many aspects. Thanks to it (and to modal decomposition) it is possible to discretize the field and to find its hamiltonian. Note that no approximation are committed and minimal hypothesis are needed. The only requirement is that the boundary conditions of the domain are spacially periodic (period $d$).

The result (equation 22 of the letter) is that the electromagnetic field is represented by a chain of coupled harmonic oscillators. The canonical variables are the couples of conjugate variables $V_n$ and $I_n$ of each individual oscillator: $$\mathcal{H}_\mathrm{em}(V_n,I_n, n \in \mathbb{Z}) = \frac{1}{2} \sum_{n \ge m} \Omega_{n-m}(V_{n}V_{m}+I_{n}I_{m})$$ The coupling coefficients are given by the Fourier series of the propagation dispersion curve $\omega(\beta d)$ of the structure: $$\omega(\beta)=\sum_{m=-\infty}^{+\infty}\Omega_m e^{m \beta d} $$ The electric (resp. magnetic) field is given by the $V_n$ (resp. $I_n$) and by a function $E_0$ (resp. $H_0$) of space which depends only on the geometry of the domain: $$E(x,t)=\sum_n V_n(t) E_0(x-nd)$$ This function of space $E_0$ is the inverse Gel'fand transformation of the modes propagating in the domain.