Why do we choose the Dirac delta function as the eigenstate of position operator? When we try to find the eigenstates of the position operator, we get that the product of (x-y) and the eigenstate must be zero. It is obvious then that for x different than y, the eigenstate must be zero.
Now for x equal to y, how do we know that the eigenstate is infinite so that we get the Dirac delta function? What if we choose any other form of the eigenstate for x=y - it would still be zero and satisfy our math, right?
 A: The position space representation of the position operator $\hat{X}$ is just the variable $x$ so the eigenfunction equation looks like $$x\psi(x)=x_0 \psi(x)$$ where $x_0$ is just a constant. The only (generalized) function that satisfies this relation is $\delta(x-x_0)$.
A: If $\psi$ were $0$ except for some finite value at $x=x_0$, then the integral of $|\psi|^2$ — that is, the probability of finding the particle anywhere at all — would be $0$.  Yet experimenters keep finding particles in various places.
A: In elementary treatments, wavefunctions are sometimes described as $\mathbb C$-valued functions which are square-integrable, i.e.
$$\int_\mathbb R |\psi(x)|^2 \mathrm dx < \infty$$
This is almost true, but not quite.  The problem is that the mathematical structure of QM is built on the concept of a Hilbert space, which is (loosely) a vector space equipped with an inner product. One of the requirements of a well-defined inner product is that it must be positive-definite, which means that if a vector $\psi$ is nonzero, then we must have that $\langle \psi,\psi\rangle >0$.
Therein lies the problem with identifying the space of wavefunctions as the vector space of square-integrable functions. If $\psi(x)$ is nonzero on a countable set of points (or more generally, on a set of points with Lebesgue measure zero), then $\langle \psi,\psi\rangle = 0$ despite the fact that $\psi$ is not the zero function. As a result, $\langle\cdot,\cdot\rangle$ does not constitute an inner product.
The solution is to identify two functions $\psi$ and $\phi$ as equivalent if the square of their difference integrates to zero:
$$\psi \sim \phi \iff \int_\mathbb R |\psi(x)-\phi(x)|^2 =0$$
Therefore, the Hilbert space of wavefunctions is smaller than the vector space of all square-integrable functions, because we are to regard two functions as representing exactly the same vector if they are equal almost everywhere.
This is a mathematical requirement, but it is neither arbitrary nor unphysical.

*

*Physically, if $\psi$ and $\phi$ are two functions which are equivalent according to the definition above, then any physical quantity computed from them is will be precisely the same. They produce exactly the same predictions, and are in every way identical in terms of physical content. It is natural to consider them to be the same.

*Mathematically, if we don't consider $\psi$ and $\phi$ to be the same, then our space of wavefunctions loses some very important properties. For example, a convergent sequence or series no longer converges to a single limit, in the sense that if $\{\psi_n\}$ converges to $\psi$ and $\phi\sim \psi$, then $\{\psi_n\}$ converges to $\phi$ as well.

The choice to identify functions as described above is not carved onto a stone tablet and brought down off of a mountain by mathematicians; it is well-motivated both physically and mathematically, and makes our lives easier rather than harder.

Given that technical aside, the answer to your question is that eigenvectors of linear operators are by definition nonzero, but the function
$$\psi(x) = \begin{cases} 1 & x=x_0\\ 0 & \text{else}\end{cases}$$
is a representation of the zero vector, because $\psi(x)=0$ almost everywhere (as the singleton set $\{x_0\}$ has measure zero).
The next logical question is then in what sense is the "delta function" an eigenvector? After all, $\delta(x-x_0)$ is not even square-integrable, so the delta function isn't even a member of the Hilbert space.  And that's true - in the standard (rigorous) formalism of QM, the position operator (and any operator with a purely continuous spectrum) does not have eigenvectors.
However, as you are learning already, it is extremely useful to talk about objects like $|x\rangle$ and $|p\rangle$, as well as plane waves and many other "wavefunctions" which don't technically belong to the Hilbert space. In order to formalize these concepts, we need to expand our horizons and develop the rigged Hilbert space formalism, in which $|x\rangle$ is called a generalized eigenvector.
To do this rigorously requires substantial technical machinery, however. As a result, elementary treatments of QM (and indeed, a great many working physicists) are content with using $|x\rangle$ and $|p\rangle$ subject to a handful of seemingly ad-hoc rules such as

*

*$\langle x|x'\rangle = \delta(x-x')$ (never write $\langle x|x\rangle$!)

*The identity operator can be written $\mathbb I = \int \mathrm dx \ |x\rangle\langle x|$

*$\langle x|\psi\rangle = \psi(x)$
The rules are simple enough to follow, but if you want to put them on a truly rigorous footing (which most physicists are not particularly interested in doing, for perfectly good reason), then it will require quite a bit of work.

A brief addendum:
There are a number of comments which suggest that mathematical rigor is irrelevant to physics - a sentiment with which I wholeheartedly disagree.
I agree that most working physicists do not need to have a complete understanding of the technical underpinnings of the structures they use in their work, in the same sense that most plumbers do not need to know the chemical structure of PVC.  However, somebody needs to know it, and the study of such things benefits plumbing as a profession - sometimes by developing new tools, and sometimes by developing a new way to use or think about the tools they already have.
Physicists do not need to be mathematically rigorous all of the time. But if a particular area of physics cannot be made rigorous (or rather, if nobody knows how to do so), then that is a sure sign that there are fundamental questions that need to be answered. And historically there have been many occasions where searching for the answers to technical questions has yielded insights and opened doors to new physics.
A: This is because the position operator is $x\delta (x-x') $. This may sound like I'm just kicking the can down the road  but I'm not.
Naively, you would expect the position operator to be $x \delta _x ^{x'}$, where the $\delta$ is the Kronecker delta. But keep in mind that $x$ is a continuous index. So the action of the operator on a state cannot be defined to be  a sum:
$$\sum _ {x'}A(x, x') f(x') $$
It's instead defined to be an integral:
$$\int A(x, x') f(x') dx'$$
You can see that if we choose $A=x\delta (x-x') $, this yields the correct action of the position operator:
$$X (f(x))=xf(x) $$
The action of the operator, defined using the integral, will be zero on any vector if you chose the Kronecker delta as the position operator.
You can try using the eigenvector as $f(x)=1$ at $x=a$ and $=0$, everywhere else and the Kronecker delta as the position operator. The action of this operator will yield zero because the integral will be zero.
This still leaves the option of having the operator as $x\delta (x-x') $ but your proposed functions as eigenvectors. The problem with this is that you cannot write an arbitrary vector $f(x)$ as a linear combination of your proposed eigenvectors:
$$f(x)=\int f(x_0) E(x, x_0) dx_0$$
$f(x_0) $ is the component of the vector projected along the eigenvector $E(x, x_0) $. You can see that $E(x, x_0) $ must be the set of delta functions to satisfy the linear combination equation.
