After revisiting my question, I think I have managed to find a proof that there are no other examples of lattices with the requested property. I'm posting it as a self-answer in case someone is interested. Here is the outline:
The case of even dimension $\mathrm{dim}\: \Lambda \ge 4$ was already described$^\dagger$ in the "Attempt" section of my question; the most important part (bounding the generalized Bernoulli numbers) is done in Section 3 of this 2011 paper by M. Johnson which I cited in the question.
We can also adapt the strategy in that paper to the odd-dimensional case. As in the even case, we start by observing that the theta series must be a certain Eisenstein series of half-integer weight $k$ and prescribed level $4N$ and character $\chi$, normalized such that the zeroth coefficient is $1$. One needs an expression for the Fourier coefficients of these series, which I found for $\mathrm{dim}\: \Lambda = 2k \ge 5$ in page 17 of this 2015 thesis by M. Owen. We are only interested in the coefficient of $q^1$, which simplifies to
$$c_1 = \left(-2\pi i/N\right)^k \:\Gamma(k)^{-1} A(1) X(1),$$
where
$$A(1) = \sum_{r=1}^{4N} \epsilon_r^{2k} \left(\frac{4N}{r}\right) \chi(r) e^{2\pi ir/N}, \quad X(1)= \prod_{l\nmid 4N}(1+\chi(l) \epsilon_l^{2k-1} l^{1/2-k})$$
(as in Owen's paper, here $l$ is prime, $\left(\frac{4N}{r}\right)$ is a Kronecker symbol, and $\epsilon_n$ denotes the principal branch of $\sqrt{\left(\frac{-1}{n}\right)}$. Note also that I use $4N$ where the author uses $N$, since for half-integer weight forms the level is always a multiple of 4). This $c_1$ is an analogue of the generalized Bernoulli numbers for half-integral weight.
The $A(1)$ factor is a sum of roots of unity whose modulus is obviously bounded by $4N$, and by the same argument as in page 9 of Johnson's paper, $|X(1)|$ is bounded by $\prod_{l}(1+l^{1/2-k})=\frac{\zeta(k-1/2)}{\zeta(2k-1)}$. Putting everything together, we see that the bound on $|c_1|$ will be lower than $2$ unless $N$ and $k$ are low enough (concretely we have either $N=1$ and $k \le 9/2$ or $1 < N \le 3$ and $k=5/2$), so recalling that $c_1$ must be a nonzero even integer, we can restrict our analysis to these cases only. All of them are then ruled out by manual computation of the coefficients.
The case $\mathrm{dim}\: \Lambda = 2$ is easily dealt with by using the relationship between 2D integer lattices and rings of integers of imaginary quadratic number fields, and applying the class number formula.
This only leaves the case $\mathrm{dim}\: \Lambda = 3$. It is known that there are no isospectral lattices in dimension lower than $4$ (see e.g. the end of Section 2 here), so in this case it suffices to appeal to the Siegel-Weil formula and check the single-class examples in dimension $3$ with the LMFDB online tool, which I already did (this could also be an alternative way to treat the case $\mathrm{dim}\: \Lambda = 2$). This completes the proof.
Remark: I admit that the proof above is somewhat "ugly", since it does not give any insight into the form of the classification. A more satisfying proof would explain why all of these lattices have the multiplicative structure of an order in a division algebra (and thus occur in dimensions 1, 2, 4, 8) except for the sporadic case $(E_8 \times E_8, D_{16}^+)$.
There is an observation to be made regarding that last case. For any lattice $\Lambda$ it is possible to define Siegel theta series of genus $g$ generalizing the usual theta series, which essentially count the number of $g$-dimensional sublattices of $\Lambda$. Since there is a bijective correspondence between lattices and their theta series in genus $g\ge \mathrm{dim} \: \Lambda$, the issue cannot arise that two lattices have the same theta series; in fact, the sporadic example of $E_8 \times E_8 \leftrightarrow D_{16}^+$ disappears already in genus 4 because of the Schottky form, and the corresponding Siegel-Weil formula for $g\ge 4$ takes the form $\mathrm{E}_8^{(g)} = \frac{405}{691} \Theta_{E_8 \times E_8}^{(g)} + \frac{286}{691} \Theta_{D_{16}^+}^{(g)}$ with $\Theta_{E_8 \times E_8}^{(g)} \neq \Theta_{D_{16}^+}^{(g)}$. In contrast, the octonionic "$\mathrm{E} = \Theta$" identity between the Siegel theta series of the $E_8$ lattice and the Siegel Eisenstein series of weight 4 does hold for all genera, and I believe the same happens for the other examples. This could perhaps be a starting point to explain why the sporadic case doesn't have a multiplicative structure.
$\dagger$- The extra possible example of dimension $4$ and level $42$ that I mentioned does not correspond to any lattice; perhaps the easiest way to see this is to inspect the coefficients of the candidate theta series, which starts as
$$1+2q+2q^2+2q^3+2q^4+12q^5+2q^6+2q^7+2q^8+2q^9+12q^{10}+\ldots$$
(the general coefficient is twice the sum of the divisors of $n$ that don't divide $42$).
Since the coefficient of $q^n$ counts the number of vectors of norm $2n$, such a lattice would have a vector $\mathbf{a}$ of squared norm $2$ and another, linearly independent vector $\mathbf{b}$ of squared norm $10$. We can also see that all vectors of squared norm lower than $10$, as well as those with squared norms strictly between $10$ and $20$, must be integer multiples of $\mathbf{a}$, since there are only two vectors of each norm in these intervals.
The squared norm of the sum $\mathbf{a}+\mathbf{b}$ is a positive even integer bounded by $(\sqrt{2}+\sqrt{10})^2 \approx 20.944$; by the above observations, it can only be $20$. But then by the parallelogram law, $|\mathbf{a}-\mathbf{b}|^2 = 2|\mathbf{a}|^2+2|\mathbf{b}|^2-|\mathbf{a}+\mathbf{b}|^2 = 4+20-20=4$, which is a contradiction since $\mathbf{a}-\mathbf{b} \not\propto \mathbf{a}$.