Hilbert space of harmonic oscillator: Countable vs uncountable?
This question was first posed to me by a friend of mine; for the subtleties involved, I love this question. :-)
The "flaw" is that you're not counting the dimension carefully. As other answers have pointed out, $\delta$-functions are not valid $\mathcal{L}^2(\mathbb{R})$ functions, so we need to define a kosher function which gives the $\delta$-function as a limiting case. This is essentially done by considering a UV regulator for your wavefunctions in space. Let's solve the simpler "particle in a box" problem, on a lattice. The answer for the harmonic oscillator will conceptually be the same. Also note that solving the problem on a lattice of size $a$ is akin to considering rectangular functions of width $a$ and unit area, as regulated versions of $\delta$-functions.
The UV-cutoff (smallest position resolution) becomes the maximum momentum possible for the particle's wavefunction and the IR-cutoff (roughly max width of wavefunction which will correspond to the size of the box) gives the minimum momentum quantum and hence the difference between levels. Now you can see that the number of states (finite) is the same in position basis and momentum basis. The subtlety is when you take the limit of small lattice spacing. Then the max momentum goes to "infinity" while the position resolution goes to zero -- but the position basis states are still countable!
In the harmonic oscillator case, the spread of the ground state (maximum spread) should correspond to the momentum quantum i.e. the lattice size in momentum space.
The physical intuition
When we consider the set of possible wavefunctions, we need them to be reasonably behaved i.e. only a countable number of discontinuities. In effect, such functions have only a countable number of degrees of freedom (unlike functions which can be very badly behaved). IIRC, this is one of the necessary conditions for a function to be fourier transformable.
ADDENDUM: See @tparker's answer for a nice explanation with a slightly more rigorous treatment justifying why wavefunctions have only countable degrees of freedom.
The Hilbert space ${\cal H}$ of the one-dimensional harmonic oscillator in the position representation is the set $L^2(\mathbb{R})={\cal L}^2(\mathbb{R})/{\cal N}$ (of equivalence classes) of square integrable functions $\psi:\mathbb{R}\to\mathbb{C}$ on the real line. The equivalence relation is modulo measurable functions that vanish a.e.
The Dirac delta distribution $\delta(x-x_{0})$ is not a function. It is a distribution. In particular, it is not square integrable, cf. this Phys.SE post.
One may prove that all infinite-dimensional separable complex Hilbert spaces are isomorphic to the set $${\ell}^{2}(\mathbb{N})~:=~\left\{(x_n)_{n\in\mathbb{N}}\mid\sum_{n\in\mathbb{N}} |x_n|^2 <\infty\right\}$$ of square integrable complex sequences.
The previous answers are all correct, but I thought I'd give a more conceptual explanation for why the delta-function basis is the "wrong" basis in which to expand when counting degrees of freedom. Since the situation is much, much more complicated in QFT, for simplicity I'll only consider first-quantized wavefunctions for a system with a fixed, finite number of particles, so that the configuration space is just $\mathbb{R}^n$ for some finite $n$. (If you don't know what "configuration space" is, all that really matters for this question is that for a single-particle system, it's the same as real space.)
Physicists often say that for these systems, "the Hilbert space $L^2(\mathbb{R}^n)$ is the space of square-integrable functions on $\mathbb{R}^n$, with inner product $\langle f | g \rangle := \int_{\mathbb{R}^n} d^nx\ f^*({\bf x})\, g({\bf x}). $" But this definition is wrong, because that isn't actually a valid inner product on that space! The problem is that it violates the positive-definiteness requirement for the inner product that $||\psi|| = 0 \implies | \psi \rangle = 0$: if a function $f$ is supported on a nonempty set of Lebesgue measure zero, then the "norm" $\int_{\mathbb{R}^n} d^nx\ |f({\bf x})|^2 = 0$. Since this "norm" is zero for some nonzero vectors, it is more properly only a seminorm on the space of square-integrable functions on $\mathbb{R}^n$. This function space is denoted $\mathcal{L}^2(\mathbb{R}^n)$ (note the different script for the "$\mathcal{L}$") and is therefore only a seminormed vector space.
To convert $\mathcal{L}^2(\mathbb{R}^n)$ into a true Hilbert space, we need to mod it out by the vector space of functions whose support has Lebesgue measure zero. In other words, we define an equivalence relation $f \sim g$ between functions $f({\bf x})$ and $g({\bf x})$ that agree almost everywhere, and then define the Hilbert space $L^2(\mathbb{R}^n)$ to be the space of equivalence classes under this equivalence relation. So two square-integrable functions $\psi({\bf x})$ and $\phi({\bf x})$ which are equal almost everywhere, but differ on a set of Lebesgue measure zero, actually correspond to the exact same state $|\psi\rangle$ in the Hilbert space. This fixes the problem, because now all those problematic functions whose support has Lebesgue measure zero correspond to the zero vector of the Hilbert space, so it's fine for them to have norm zero.
This is more than just a technical trick only performed in order to satisfy the mathematical definition of an inner product - it's actually the right thing to do physically. Remember that the value of $|\psi({\bf x})|^2$ at a particular point ${\bf x}$ isn't actually a probability - it's a probability density, which is not a directly physical quantity. You can't directly measure the probability density at a single point; you can only measure the probability $P(V) = \int_V d^nx\, |\psi({\bf x})|^2$ for a particle to be found in a (potentially very small) region $V$. But if two wavefunctions $\psi, \phi \in \mathcal{L}^2(\mathbb{R}^n)$ only differ on a set of Lebesgue measure zero, then $P(V) = \int_V d^nx\, |\psi({\bf x})|^2 = \int_V d^nx\, |\phi({\bf x})|^2$ will be the same for any region $V$. Therefore all physically measurable quantities will be the same for these two wavefunctions, and so they correspond to the same physical state $| \psi \rangle \in L^2(\mathbb{R}^n)$.
The point of all this is that any wavefunction $\psi({\bf x})$ carries a whole lot of extra, unphysical information (beyond just the overall phase factor, which you're probably used to). Changing its value at any set of points of Lebesgue measure zero doesn't actually change the state. The (uncountable) delta-function basis is too "fine" and picks out all these irrelevant unphysical degrees of freedom. The (countable) oscillator-eigenstate basis, on the other hand, is much less sensitive to the details of the wavefunction: changing $\psi({\bf x})$ on any set of Lebesgue measure zero doesn't change any of the expansion coefficients $\langle \psi_n | \psi \rangle$. These coefficients therefore only record information about the physical degrees of freedom, of which there are only countably many.
By the way, the Hilbert space $L^2 \left( \mathbb{R}^d \right)$ is the same for the free particle as for the harmonic oscillator, so everything in this answer carries over directly to the companion question about the free-particle Hilbert space.