Partial trace over continuous degrees of freedom Let us suppose we have some quantum system whose Hilbert space admits a bipartition $\mathscr{H}\simeq \mathscr{H}_A\otimes \mathscr{H}_B$. Let $|n\rangle_A$ be a basis of $\mathscr{H}_A$ and $|m\rangle_B$ be a basis of $\mathscr{H}_B$. Then $|n,m\rangle_{A,B}$ is a basis of the composite system.
A general mixed state $\rho$ can be written always as $$\rho=\sum \rho_{nmn'm'}|n,m\rangle_{AB}\langle n',m'|$$
And we can define the partial trace over $\mathscr{H}_A$ to be $$\operatorname{Tr}_A\rho=\sum \rho_{mm'}|m\rangle_B\langle m'|,\quad \rho_{mm'}=\sum_n \rho_{nmnm'}$$
This seems to work well if $|n\rangle_A,|m\rangle_B$ are discrete basis.
Now, in Quantum Mechanics we often work with "continuous bases". These are not really rigorous since the "position eigenstates" $|x\rangle$ are not well defined. Still, this is very useful.
So now suppose that for $\mathscr{H}_A$ we take such a $|x\rangle_A$ basis. Then a basis of $\mathscr{H}$ is $|x,m\rangle_{A,B}$ and the state is $$\rho=\sum \int \rho_{mm'}(x,x') |x,m\rangle_{A,B}\langle x',m'|$$
And we can in analogy take a partial trace $$\operatorname{Tr}_A\rho=\sum \rho_{mm'}|m\rangle_B\langle m'|,\quad \rho_{mm'}=\int \rho_{mm'}(x,x)dx$$
It seems everything is fine, but now take the pure state $\rho = |z,\sigma\rangle\langle z,\sigma|$ (for instance, a particle located at $z$ with spin $S_z$ being $\sigma$). If we evaluate $$\rho_{mm'}(x,x')=\delta_{m\sigma}\delta_{m'\sigma}\delta(x-z)\delta(x'-z)$$
and the partial trace involves $(\delta(x-z))^2$ which is ill-defined.
Still, there should certainly be a way of talking about partial traces when continuous degrees of freedom like position and momentum are involved. 
How is that done correctly? How to avoid this issue of getting a delta function squared when taking a partial trace over continuous degrees of freedom?
 A: I wouldn't say that the square of the delta function is ill-defined, there are ways to define functions of $\delta$, see e.g. this question. Rather, the problem is that it's the wrong result.
You are essentially looking for a formal way to write a $\rho(x,y)$ such that:
$$\int dxdy \,\rho(x,y) |x\rangle\!\langle y|=|z\rangle\!\langle z|.\tag A$$
I guess you got the expression with the square of the Dirac delta by using
$\langle x|z\rangle=\delta(z-x)$.
This expression is, however, incorrect: integrating $\langle x|z\rangle$ over $x$ (or $z$) does not give you the identity. It is therefore more correct to use something like
$$\langle x|z\rangle=\sqrt{\delta(x-z)},\tag B$$
so that $\int dx |\langle x|z\rangle|^2=1$, as it should be. There is no problem in defining this object if you think of distributions as "underdetermined functions" (that is, functions only defined via their integral properties): $\sqrt{\delta(x)}$ is a function whose square integrates to the identity and is zero outside of $x=0$. On a more rigorous level, you can make sense of these objects e.g. via holomorphic calculus, see this question.
Using (B), you get $\rho(x,y)=\delta(z-x)^{1/2}\delta(z-y)^{1/2}$.
Integrating, you get
$$\int dxdy\,\delta(z-x)^{1/2}\delta(z-y)^{1/2}|x\rangle\!\langle y|=|x\rangle\!\langle x|.$$
You can probably prove that this holds using specific properties of delta functions. A simple, if handwavy, argument is is based on understanding the delta functions as the limit of functions concentrated at $x=0$.
Using rectangle functions for ease of calculation, you then have
$$\int dxdy \, \delta(z-x)^{1/2}\delta(z-y)^{1/2}f(x,y)
\simeq f(z,z) \Delta \frac{1}{\sqrt\Delta}\frac{1}{\sqrt\Delta}=f(z,z).$$
What I'm doing here is using the fact that $\delta(x)=0$ for $x\neq0$ to restrict the domain of integration to a "very small" region $\Delta$ around, in this case, $z$. If $f$ is continuous, then for $\Delta$ small enough I can always approximate $f(x,y)\sim f(z,z)$. Then, using rectangular approximations for the deltas, we have that both $\delta$ equal $1/\sqrt\Delta$ in this infinitesimal region. Because there are two such $\delta$, we get a factor $(1/\sqrt\Delta)^2$, which cancels out with the $\Delta$ coming out of the integration region.
