Thanks for all the comments. Using the paper linked by the helpful commenter (to which I found a pdf version, see here) - I could see how Planck did it (reminding me of all my old statistical mechanics in the process). Thought I'd sketch it out here to answer the OP.
We want to identify the energy distribution of a blackbody. Planck conceptualised the situation as a system of oscillators capable of absorbing and emitting quanta of radiation. In thermal equilibrium, he used some thermodynamic arguments to arrive at the distribution of energy of such a system.
First, define the overall energy of the system of $N$ oscillators as $E_N = NE$, where E is the average energy of one oscillator. Similarly the entropy is extensive: total entropy of the system is $S_N = NS$.
The quantisation step is to realise that to arrive at the known right answer (Planck found the energy distribution essentially empirically, before working out how to derive it) you must assume the energy stored by any oscillator to be quantised (not continuous) such that the total energy $E_N$ can be decomposed into $P$ units of a quanta of energy $\epsilon$: $E_N = P\epsilon$.
The number of ways of arranging the $P$ units across all $N$ oscillators is (as written above) given by the weight function:
$$
W = \frac{(N+P-1)!}{P!(N-1)!} \simeq \frac{(N+P)!}{P!N!} \simeq \frac{(N+P)^{N+P}}{P^PN^N}
$$
where we've assumed number of oscillators and number of energy units are large and then applied Stirling's approximation.
What I'd forgotten was the clever steps to utilise this within Boltzmann's equation (which was actually written by Planck for the first time) to get the entropy of the set of N resonators:
\begin{align}
S_N = k \ln W &= k \big[ (N+P)\ln (N+P) - N\ln N - P\ln P \big] \\
&= kN \big[ \left(1+\frac{P}{N}\right) \left( \ln N + \ln\left(1+\frac{P}{N}\right) \right) - \ln N - \frac{P}{N}\ln P \big] \\
&= kN \big[ \left(1+\frac{P}{N}\right) \ln\left(1+\frac{P}{N}\right) - \frac{P}{N}\ln \frac{P}{N} \big].
\end{align}
Using $P = E_N/\epsilon = NE/\epsilon$ and dividing out the number of oscillators to get the average entropy per oscillator we get:
$$
S = k \left[ \left(1+\frac{E}{\epsilon}\right) \ln\left(1+\frac{E}{\epsilon}\right) - \frac{E}{\epsilon}\ln \frac{E}{\epsilon} \right]
$$
This is all framed in terms of the unknown energy quanta Planck introduced to be able to generate the form of the weight function he knew would lead to the correct result. To go beyond this he uses Wien's displacement law, which states that the peak blackbody intensity occurs at a frequency proportional to temperature of the blackbody - i.e. $\nu_{peak} \propto T$.
Planck goes through a series of steps (which I have to admit I haven't looked at in detail) to show that this implies that the entropy of a resonator sitting in a diathermic medium (which I think means basically in thermodynamic equilibrium) must be a function of $E/\nu$, i.e. $S = f(E/\nu)$.
Our expression for $S$ is framed in terms of $E/\epsilon$ which means that (up to a constant of proportionality) $\epsilon \propto \nu$, or $\epsilon = h\nu$ introducing Planck's constant. This 1901 paper (and a 1900 paper that seems to cover mostly the same ground) I think are some of the first places where the Planck constant is introduced and quantified - which I thought was quite cool!
Anyway, substituting this into the entropy equation we can find an expression for the energy of an average oscillator related to the temperature $T$:
$$
\frac{1}{T} = \frac{dS}{dE} = \frac{k}{h\nu} \ln\left( \frac{1 + E/h\nu}{E/h\nu} \right) = \frac{k}{\epsilon} \ln\left( 1 + \frac{h\nu}{E} \right) \Rightarrow E = \frac{h\nu}{e^{h\nu/kT}-1}.
$$
To go all the way to the energy radiance, it is possible to re-use the Rayleigh-Jeans law, namely that the radiance of a blackbody is $B_\nu = \frac{2\nu^2}{c^2}kT = \frac{2\nu^2}{c^2} E$, where $E$ is the average energy of a mode of radiation in the blackbody. This expression is incorrect since it assumes that each mode of radiation has (according to equipartition of energy) average energy $kT$. The reason I'm borrowing it here is because it contains the results of doing various integrals to deal with the 3D geometry so we don't have to (Planck refers to a previous paper of his to sidestep doing this in his paper).
As a result we can convert from energy of a mode of radiation to spectral energy radiance - i.e. the result quoted above in the OP - using Planck's expression for the average energy (as derived by considering the distribution of energy quanta across oscillators in the cavity wall), giving (finally) the result:
$$
B_\nu = \frac{2\nu^2}{c^2} E = \frac{2h\nu^3}{c^2}\frac{1}{e^{h\nu/kT}-1}.
$$
It is worth noting in the paper linked above that Planck writes out the spectral energy density not the energy radiance as I've done - the two can be switched between according to the relation: $u_\nu = \frac{4\pi}{c}B_\nu$.
