If I understand correctly, your question basically comes down to identifying a basis for the space of square-integrable functions, $L^2(\mathbb{R})$, since any physical state $|\Psi\rangle$ can be constructed by performing the integral you listed in your question with a function $\Psi_x(x)\in L^2(\mathbb{R})$. $L^2$ is known to be a vector space, so a basis has to exist. Off the top of my head, I think an example would be
$$f_k(x) = e^{-x^2/a^2 - ikx}$$
which is just the normal plane wave basis $e^{-ikx}$ multiplied by a Gaussian envelope $e^{-x^2/a^2}$ where $a$ is some constant. Multiplying by this Gaussian envelope ensures that the functions will be square-integrable, but since you use the same envelope for every element of the basis, you can factor it out of the Fourier transform, so it doesn't change any of the essential properties of momentum-space decomposition.
P.S. I found a question on math.SE that seems related, and which motivated this answer.
There are many approaches. But first I want to make sure you know that when you have $n$ spin zero particles in a 3d space and have a wavefunction that the function is from $\mathbb R^{3n},$ i.e. from configuration space. But also I want you to know when someone says infinite dimensional Hilbert Space, they mean the size of a set of mutually orthonormal vectors.
Now when you have some particles and know their spin you could start with a specific concrete set of wavefunctions from configuration space into a tensor product of the spin states. You can then use the inner product on the spin states to get an inner product of the wavefunctions. And this is basically like picking a basis when do vector analysis.
Given the wavefunctions you could then talk about the operators. The operators form a vector space, then have a norm defined on them, they have an adjoint operation. And some times $A-\lambda I$ is not invertible as an operator (for instance when $\lambda$ is an eigenvalue of $A$). And it is possible to reverse the process.
So if you start with a bunch of things that act like operators, i.e. that form a vector space, that have a multiplication, have a norm, have a star operator and they satisfy the rules you'd expect. Then you can call it a $C^*$ algebra. Then it turns out the standard assumptions in mathematics assert that there is a Hilbert space and a subset of the operators on the Hilbert space that as a vector space with a multiplication, a norm, and the adjoint (as the star operator) will act just like your $C^*$ algebra.
So you can start with a Hilbert Space and then define operators. Or you can start with a $C^*$ algebra and then define your Hilbert Space as one of the Hilbert Spaces that has a subalgebra of operators just like your $C^*$ algebra. The biggest deal is knowing how to connect to experimental results since either way you have an algebra, some operators and a Hilbert Space and any triple you have is going to look like the other triple you made starting from the other end.
As for a spectrum. Since the algebra has a multiplication you can discuss whether $A-\lambda I$ is invertible without needing to have vectors that are eigen to any operators, the object just have an inverse or it doesn't, from the set of things you can multiply.
You could imagine messing with matrices, adding, scaling, multiplying, taking adjoints. And discussing whether they have inverses, without every mentioning that there are vectors the matrices could act on. And then you can be more abstract and say the fact they are matrices wasn't important, just have things with a scaling operation, a norm operation, and adjoint operation, and a multiplication operation, that do the right things.
How does one choose a Hilbert space to associate with a system?
You could start with it, or start with your algebra. Or you can focus on connecting a Hilbert Space with the outcomes of certain experimental results. For instance the completion of the set of eigenvectors for a maximal set of compatible measurements.
How does one deal with this ambiguity given by that theorem?
Having different vector spaces that correspond to your physical setup is no worse than being able to choose your x y and z axis to point any direction you want in the lab. Truly no different. They are all isomorphic, and 3d vector space. And it isn't a problem.
If the observables determine the space, how does this happen if one can only define the operators after knowing the space?
The algebra doesn't need the space. You can have a bunch of square matrices without having column or row vectors. You can have the algebra with the scaling, adding, multiplying, morning, and adjointing without having the matrices. And you still can have a spectrum.
Best Answer
Your idea of what $\psi(\vec x,t)$ is supposed to be is essentially correct. Given a space of states $\mathcal{H}$, the "Schrödinger state" is a map $$ \psi : \mathbb{R}\to\mathcal{H}, t\mapsto\lvert\psi(t)\rangle$$ where $\lvert \psi(t)\rangle\in\mathcal{H}$ for every instant $t$. If $\mathcal{H}$ is a space of functions in a variable $\vec x$, then $\lvert \psi(t)\rangle$ is often written $\psi(\vec x,t)$.
The space of wavefunctions in usual quantum mechanics is crucially not the space of smooth functions $C^\infty(\mathbb{R}^3,\mathbb{C})$, but the space of equivalence classes of square integrable functions $L^2(\mathbb{R^3},\mathbb{C})$(link to Wikipedia article on $L^p$-spaces). This space contains the smooth compactly supported functions $C_c^\infty(\mathbb{R}^3,\mathbb{C})$ (non-compactly supported functions are not necessarily square-integrable) and also all smooth square-integrable functions, but is larger.
The reason for this is twofold: For one, there is no reason to demand a generic wavefunction be continuous - the whole physical content of the wavefunction is encapsulated in the probability density $$ \rho(\vec x) = \lvert\psi(x)\rvert^2$$ and this needs to have $\int\rho(\vec x)\mathrm{d}^3x = 1$ to be probability density, hence $\psi$ must be square integrable, but nothing else.
Another reason is that the smooth compactly supported functions do not form a Hilbert space under the scalar product $$ (f,g) = \int \overline{f(x)}g(x)\mathrm{d}^3 x$$ since they are not complete - there are sequences which are Cauchy with respect to the $L^2$-norm induced by this inner product, but which do not converge. The space of square-integrable functions is precisely the completion of the smooth compactly supported functions under this norm, and hence a Hilbert space. In other words, every wavefunction may be arbitrarily accurately be approximated by a smooth compactly supported function, but is itself not guaranteed to be a smooth function.
Among other difficulties commonly overlooked, this means that, strictly speaking, one cannot evaluate wavefunctions at points, since the $L^2$-space elements are only defined up to a zero measure set, and points are of zero measure. This is again meaningful when considering $\rho(\vec x)$: The value of a probability density on a zero measure set is physically meaningless, since only the integration of it over a set of non-zero measure gives a physically relevant probability.