Perhaps surprisingly, the random variable $\ell(n)$ (with $n$ drawn uniformly from $[0,N)$ concentrates too much around $\log_2\log_2 N$ (where $\log_2$ denotes the logarithm to base $2$) to have a limiting parity probability - the variance stays bounded as $N \to \infty$, as opposed to growing to infinity. One only recovers a limiting law when the fractional part $\{\frac{1}{2} \log_2 \log_2 N \}$ of half the double logarithm of $N$ converges to a limit, and when one does so the parity probability will usually converge to a limit that deviates slightly from $1/2$.
To simplify the calculations a little let us assume that $N$ is of the form $N = 2^{2^{2k+1}}$ (I'll leave it as an exercise to the reader to handle the general case) in the asymptotic regime $k \to \infty$. Then the binary expansion of a randomly chosen element $n$ of $[0,N)$ consists of $2^{2k+1}$ independent Bernoulli variables (each taking $0$ and $1$ with values $1/2$). We think of this as the initial segment of an infinite sequence of Bernoulli variables. Now we perform the standard trick of viewing this sequence as a renewal process. After each $0$, the number of $1$s one encounters before one reaches the next $0$ is $a-1$ where $a$ is distributed according to a geometric distribution of expectation $2$. One can thus interpret this sequence as $a_1-1$ zeroes followed by a one, then $a_2-1$ zeroes followed by a one, and so forth ad infinitum, where $a_1,a_2,\dots$ are iid geometric distributions of expectation $2$. By the law of large numbers, we see with probability $1-o(1)$ that the first $t$ for which $a_1+\dots+a_t$ exceeds $2^{2k+1}$ will lie in the range $[2^{2k}-2^{4k/3},2^{2k}+2^{4k/3}]$ (say). Also, by symmetry we see that with probability $1-o(1)$, the maximum value of the $a_i$ for $i \leq 2^{2k}+2^{4k/3}$ will already be attained for $i \leq 2^{2k}-2^{4k/3}$. Putting these two together, we see that with probability $1-o(1)$, $\ell(n)$ will equal $\sup_{1 \leq i \leq 2^{2k}} a_i-1$. So asymptotically we just need to understand the distribution of $\sup_{1 \leq i \leq 2^{2k}} a_i-1$. We have the exact formula $$ {\bf P}( \sup_{1 \leq i \leq 2^{2k}} a_i-1 < t ) = \prod_{i=1}^{2^{2k}} {\bf P}(a_i-1 < t)$$ $$ = (1-2^{-t})^{2^{2k}}$$ for any positive integer $t$, so in particular $$ {\bf P}( \sup_{1 \leq i \leq 2^{2k}} a_i-1 - 2k < s ) = \exp( - 2^{-s} ) + o(1)$$ for any fixed $s$. Thus in the limit $k \to \infty$, $\ell(n) - 2k$ converges in distribution to a discrete random variable $X$ with distribution function $$ {\bf P}( X < s ) = \exp( - 2^{-s} ).$$ (Is there a name for this sort of random variable? EDIT: it is a discrete Gumbel distribution, see update below.) The quantity $\frac{|E \cap [0,N)|}{N}$ then converges to the probability that $X$ is even, which is $$ \sum_{j \in {\bf Z}} \exp(-2^{-2j-1}) - \exp(-2^{-2j}) = 0.4998402\dots$$ which is very slightly less than $1/2$. (If one picked a different subsequence of $N$ one would obtain a different limit; for instance if $N = 2^{2^{2k}}$ then the same analysis would ultimately give the complementary limiting probability of $0.500158\dots$.)
UPDATE: after a tip in the comments, I'll remark that the above analysis will eventually show (after some additional effort) that the distribution of $\ell(n)$ is asymptotic to the integer part $\lfloor \mathrm{Gumbel}(\log_2 \log_2 N, \log_2 e)\rfloor$ of a Gumbel distribution, in the sense that the Levy metric between the two metrics goes to zero as $N \to \infty$ (without any further restriction on the natural number $N$). In retrospect this was a natural guess, given the usual role of the Gumbel distribution in extreme value theory.