$$\newcommand{\a}{\mathbf{\ddot{\vec{r}}}} \newcommand{\v}{\mathbf{\dot{\vec{r}}}} \newcommand{\x}{\times} \newcommand{\r}{\mathbf{\vec{r}}} \newcommand{\h}{\mathbf{\vec{h}}} \newcommand{\e}{\mathbf{\vec{e}}} \newcommand{\L}{\mathbf{\vec{L}}} \newcommand{\p}{\mathbf{\vec{p}}}\newcommand{\A}{\mathbf{\vec{A}}}\newcommand{\o}{\mathbf{\cdot}}$$
The eccentricity vector ($\e$) is a conserved quantity (somewhat long and meandering proof at end), so you can compute its magnitude anywhere in the orbit you like, but there are certain places where the answer is easier to figure out than others. Hunting for the easiest of them all brings you to periapsis, so that's what we use to define the standard formula.
$\h=\r\x\v$, so when you write its scalar magnitude, it's easiest to do that at apoapsis and periapsis, because then the answer is just $r\dot{r}$, the magnitude of the radius times the magnitude the velocity, with no trig function of angle term, because only along the line of apsides are $\r$ and $\v$ perpendicular. Then when you try to compute the magnitude of $\r\x\h = \v\x(\r\x\v)$, you again find it easiest to do at either periapsis or apoapsis, as $r\dot{r}^2$. Which of the two it points toward seems a matter of taste, but that it should point to one of them seems natural to me.
The proof that David Hammen's comment suggested I include looks like the result of what one might call a brute force search for conserved quantities, in which the sort of vectors one usually finds lying about in mechanics problems (position, velocity, momentum...) are combined in various ways, and we take time derivatives of them all, in search of something that equals zero. Why we take the derivatives of these particular combinations, rather than many others we might consider, is because these are the ones that led physicists such as Laplace (in 1799) to find interesting zeros and publish them. It looks like sleight of hand --- making seemingly random combinations of the position and velocity vectors just to see what happens, and then acting surprised to have stumbled on to a formula that gives the equation of the orbit --- but it's really just knowing the answer and working backwards.
Most textbooks follow the same basic approach; for example, the way Howard Curtis handles it in Orbital Mechanics for Engineering Students (4th edition, pages 67-74, 2020) closely follows the method of Herbert Goldstein's Classical Mechanics (2nd edition, pages 102-105, 1980), but adds a bunch of pictures and shows a lot more intermediate steps in the algebra. I learned this from Goldstein, but he starts from the Euler-Lagrange equation, which I won't assume everyone knows. I will, however, keep things as general as I can for as long as I can.
Physicists like momentum, both linear momentum ($\p$) and angular momentum ($\L=\r\x\p$), partly because in many situations of practical interest, one or both are conserved, meaning their values do not change with time. Expressions for $\p$ can get a little complicated, such as when electromagnetism is introduced, but for now we consider the simplest case, $\p=m\v$. Then $\L=m\r\x\v$, and $\h=\L/m=\r\x\v$. Consider now the derivative of $\h$ with respect to time, by means of the product rule:
$$ \frac{d\h}{dt} = \frac{d}{dt}\left(\r\x\v\right) = \frac{d\r}{dt}\x\v + \r\x\frac{d\v}{dt} = \v\x\v + \r\x\a$$
We know $\v\x\v$ is always zero by definition of the cross product. For any force which acts parallel to the radial position vector, $m\a=\mathbf{\vec{F}}=f(r)\r/r=f(r)\mathbf{\hat{r}}$, so $\r\x\a$ is proportional to $\r\x\r$, which is also always zero by definition, so for any such force (which we call central), angular momentum is conserved ($\mathbf{\dot{\vec{L}}}=m\mathbf{\dot{\vec{h}}}=0$). Since that worked well for us, let's do it again, but considering two other combinations.
First, take the time derivative of the squared magnitude of the position vector, $r^2=\|\r\|^2=\r\o\r.$ We have $$\frac{d}{dt}r^2 = 2r\dot{r}\ \ \mathrm{and} \ \ \frac{d}{dt}\r\o\r = \v\mathbf{\cdot}\r + \r\o\v = 2\r\o\v\ \ ,$$ which combine to tell us that $\r\o\v=r\dot{r}$, in general. This is another way of saying motion where $\v$ is always perpendicular to $\r$ is motion where $\dot{r}=0$, meaning the path is a circle. Be careful here: $\dot{r}$ is the derivative of the magnitude of the position vector, which does not generally equal the magnitude of the derivative of the position vector (which is the speed, and in circular motion, can be anything at all). Note that if the coordinates are defined with respect to an arbitrary observer, $r$ is the range of the moving object to the observer, and $\dot{r}$ is the range rate at which the object is getting closer (negative $\dot{r}$), getting farther away (positive $\dot{r}$) or staying the same distance away (because either $\v=0$ or $\r\o\v=0$), so $\dot{r}=\r\o\v/r$ gives the Doppler shift times the speed of the signal.
Second, take the time derivative of $\p\x\L = m\v\x m\h$, which equals $m^2\a\x\h$ since $\mathbf{\dot{\vec{h}}}=0$. Now the fun begins. If we don't know anything about the force law, we'd have to stop. But, for a central force, we know $m\a=f(r)\r/r$, so $$m^2\a\x\h = \frac{mf(r)}{r} \r \x (\r\x\v) = \frac{mf(r)}{r} \left[ \r (\r\mathbf{\cdot}\v) - \v(\r\o\r)\right] = \frac{mf(r)}{r} \left( r\dot{r}\r - r^2\v\right)$$ by means of a vector identity often called the "BAC-CAB" rule.
Taking yet another time derivative seemingly at random, $$\frac{d}{dt}\left(\frac{\r}{r}\right)=\frac{r\v-\dot{r}\r}{r^2} = \frac{\a\x\h}{r^3} \frac{r}{mf(r)}$$ This means that if $f(r)=-k/r^2$ for some constant $k$, the inverse square force we all know and love (including Gauss's law in electricity as well as Newton's law of gravitation), then $$\frac{d}{dt}\left(\p\x\L\right)=\frac{d}{dt}\left(\frac{mk\r}{r}\right)$$ so $\A=\p\x\L-mk\r/r$ has time derivative zero, and thus is the conserved quantity we were seeking.
Now, what do we know about $\A$? Well, $\A\o\L=0$, since $\L$ is perpendicular to $\p\x\L$ and $\r$ is perpendicular to $\L=\r\x\p$, so $\A$ must lie in the same plane as $\r$ and $\v$. Now, let's write $\A\o\r = Ar\cos\theta$ in the usual way, with $\theta$ as the angle between the vectors $\A$ and $\r$. Since the triple scalar product is cyclic, we can say $(\p\x\L)\o\r = (\r\x\p)\o\L = \L\o\L$, so $Ar\cos\theta=L^2-mkr$. Rewriting some more brings us to $$\frac{1}{r}=\frac{mk}{L^2}\left(1+\frac{A}{mk}\cos\theta\right)$$
The last step is some fiddling with energy balance equations to discover that the eccentricity of the orbit ellipse, in terms of its semimajor axis $a$, angular momentum $L$, mass $m$ and force constant $k$ is $e=\sqrt{1-L^2/mka}$, so $a(1-e^2)=L^2/mk$, which means our equation for $\A\o\r$ gives exactly the elliptical orbit formula $r=a(1-e^2)/(1+e\cos\theta)$, so long as $A=mke$ and $\theta$ is the true anomaly, which means $\A$ points to periapsis.
This proof got harder to follow as I went along, and the last part is perilously close to "and then a miracle happens", but I have written too much for today already. If someone else sees a good way to clean this up, please go ahead.
