Use of the Dirac delta as a sampling operator
The OP's "equation"
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\tag{1}
$$
is incorrect; the sum on the right side of $(1)$ is the sampled signal $x_s(t)$, not $x(nT)$ at all.  That is,
$$x_s(t) = \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right).\tag{2}$$
This is not an ordinary function since it is a weighted sum of Dirac deltas, but if we pretend that it is indeed an ordinary function, then the  "value" of $x_s(t)$ at $t=mT$ is
$$x_s(mT) = \sum_{n=-\infty}^{\infty}x\left(nT\right)\delta\left(mT-nT\right) = x(mT)\delta(0),$$
while if $t$ is not an integer multiple of $T$, then the "value" of $x_s(t)$ is $0$. Note the quotation marks around the word value and avoid the usual nonsense of saying that $\delta(0) = \infty$ that is found in so many low-level DSP books.  $x_s(t)$ is not a function and it does not have any value at any real number $t$. What the sampled signal
$x_s(t)$ as defined in $(2)$ is is an impulse train; a collection of impulses spaced $T$ seconds apart with coefficient or amplitude of the impulse at time $t=mT$ being $x(mT)$, that is, the sampled value of $x(t)$ at time $mT$. As stated earlier, $x_s(t)$ is not an ordinary function of time; it is, as all the other answers have pointed out, a functional, that is, a function that operates on functions as a whole, not on individual specific values that the function being operated on takes on as ordinary functions do.   Now, it is a standard notion that
$$x(t)\delta(t-t_0) = x(t_0)\delta(t-t_0)$$ provided that $x(t)$ is continuous at $t=t_0$, and so if we re-write $(2)$  as
$$x_s(t) = \sum_{n=-\infty}^\infty x\left(nT\right)\delta\left(t-nT\right)
\tag{3}$$ we see that $(3)$ formally resembles the commonly used sum
$$x_m(t) = \sum_{n=-\infty}^{\infty}x\left(nT\right)g\left(t-nT\right)\tag{4}$$
which describes the modulation of a stream of pulses (of the form $g(t)$ spaced $T$ seconds apart) by the sample values of $x(t)$ spaced $T$ seconds apart. Small wonder then that textbook writers notice the analogy between $x_s(t)$ as given in $(3)$ and $x_m(t)$ in $(4)$ and impulsively describe the process that leads to $(2)$ or $(3)$ as a modulation process instead of a multiplication process. (At this point, it is worth reflecting that several modulation processes (e.g. AM-DSBSC which can be described as $x(t) \mapsto 
x(t)\cos(2\pi f_ct)$) are indeed just multiplication of a carrier signal by a baseband signal).

So, why is this rigamarole useful to anyone?  Well, let's look at the Fourier transform of $x_s(t)$:

\begin{align}
 \int_{-\infty}^\infty x_s(t) \exp(-j2\pi ft) \,\mathrm dt
&amp;= \int_{-\infty}^\infty \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt\\
&amp;= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt
&amp;{\scriptstyle{\text{sifting property of Dirac delta}}}\\
&amp;= \sum_{n=-\infty}^{\infty} x(nT)\exp(j2\pi f(nT)) \tag{5}
\end{align}
which is the DTFT of the discrete-time sequence $\big\{x(nT)\big\}_{n=-\infty}^\infty$ of sample values of $x(t)$ spaced $T$ seconds apart. More generally, if we have a (linear time-invariant BIBO) filter whose impulse response (dirty word!) is $g(t)$, then applying the input $x_s(t)$ to this filter, we get that the output is
\begin{align}
\int_{-\infty}^\infty \left( \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\right) g(\tau-t)
\,\mathrm dt
&amp;= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) g(\tau-t)
\,\mathrm dt\\
&amp;= \sum_{n=-\infty}^{\infty} x(nT) g(\tau-nT)\\
&amp;= x_m(\tau)&amp; \text{as defined in }(4)
\end{align}
So, Equation $(2)$ does have some utility whereas the jury is still out on the OP's $(1)$.

The OP asks:


  Why isn't the sampling process modeled as
  $$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau ??
$$


Well, quite apart from the royal confusion that the OP is creating by using $n$, a parameter on the left side of the OP's alleged definition, as an index of summation, note that the sifting integral property says that the value of that integral is $x(nT)$.  Consequently, the OP's purported definition reduces to
$$x(nT) = \sum_{n=-\infty}^\infty x(nT),$$
which is nonsensical, and matters don't improve very much if a different symbol is used for the parameter on the left side because all one gets is
$$x(mT) = \sum_{n=-\infty}^\infty x(nT),$$
which suggests that $x(nT) = 0$ for all integers $n$ that are not equal to the chosen $m$.  That's no way to run a railroad!



Turning to the frequency domain, the OP asserts that


  If 
  $$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
  is a sinusoid, then
  $$
X\left(f\right)=\delta\left(f-f_{0}\right) \tag{4}
$$
  and
  $$
Y\left(f\right)=H\left(f_{0}\right).\tag{5}
$$


While $(4)$ is correct, $(5)$ is not. Actually,
$$Y(f) = H(f)X(f) = H(f)\delta(f-f_0) = H(f_0)\delta(f-f_0)\tag{6}$$ which does not equal
$H(f_0)$ as the OP claims.  Indeed, in $(5)$ the left side is a function of $f$ while the right side is a constant so that the standard interpretation of $(5)$ is that $Y(f)$ is a constant, which in turn implies that
$y(t)$ is an impulse or Dirac delta!  Since $H(f)$ is arbitrarily chosen, what $(5)$ seems to be implying is that the response of any filter to a complex sinusoid is an impulse, which is manifestly nonsense.  What is true is that if we assume the system with transfer function $H(f)$ (a.k.a frequency response) has input $\exp(j2\pi f_0 t)$, then the output is 
$$y(t) = \mathcal F^{-1}(Y(f)) = \int_{\infty}^\infty H(f)\delta(f-f_0)\exp(j2\pi tf)\, \mathrm df = H(f_0)\exp(j2\pi f_0t)$$
which is the epitome of the definition of the frequency response function: 


  If for each real number $f_0$, the complex sinusoidal input $\exp(j2\pi f_0t)$ at frequency $f_0$ produces complex sinusoidal output $H(f_0)\exp(j2\pi f_0t) = |H(f_0)|\exp\big(j(2\pi f_0t+\angle H(f_0))\big)$ (which is also at frequency $f_0$ but has amplitude
  $|H(f)|$ and phase $\angle H(f_0)$), then $H(f)$ is called the frequency response function or transfer function of the linear time-invariant system.


In short, the OP's whole question is based on massaging "alternative facts" (that do not withstand scrutiny) to arrive at incorrect conclusions.