I started writing this before the accepted answer was posted, and am posting it anyway for the sake of completeness.
I will do my best to provide a comprehensive answer for the two questions which were asked here. I will not be addressing the shape of the earth, the utility (or lack thereof) of the proposed model, or any broader implications.
Question 1: Is the "proof" provided in the linked blog post valid?
Difficult to say. It's unclear exactly what the author is assuming and what they intend to "prove." If I were really invested in it, I would contact the author and ask for clarification. I am not that invested, however, so I'm going to provide my best guess as to what I think the author is trying to say using the interpretation which makes the post "as right as possible."
Premises:
All mass in the universe is distributed uniformly through a region of space bounded by two parallel planes separated by a distance $d$. One of these planes is the surface of the earth, which is assumed flat and infinite in area.
Gravitation is dictated by Gauss's law for gravity. The law is given by the equation:
$${\subset\!\supset}\llap\iint_{\partial\Omega}\mathbf g\cdot\mathbf n\ dA=-4\pi Gm\qquad m=\iiint_\Omega\rho\ dV,$$
where $\Omega$ is a region bounded by a closed surface $\partial\Omega$, $\mathbf g$ is the gravitational acceleration at a point, $\mathbf n$ is the [outward-facing] unit normal vector at a point on $\partial\Omega$, $m$ is the mass within $\Omega$, and $\rho$ is the density at a point. The second equation states that the mass within a three dimensional region is the integral of the density over that region. When this density is the same everywhere, this is equalent to $m=\rho V$, where $V$ is the volume of $\Omega$.
[Note: The author makes a minor mistake in confusing volume density (the standard notion of "density" in units of mass per unit volume) with area density, but this is easily corrected.]
Neither premise is in agreement with observation.
The first implies, among other things, that the sun, moon, and stars all have negligible masses and so do not exhibit any significant gravitational interaction. The second implies, among other things, the reality of instantaneous communication. However, the validity of the premises do not affect the validity of the conditional assertion which I have decided the author is arguing.
Conclusion:
If the premises hold, then the gravitational acceleration at any point above the surface of the earth is $\mathbf g=(0,0,-g)\ \text m\cdot \text s^{-2}$ - that is, $g\ \text m\cdot \text s^{-2}$ perpendicular to and in the direction of the surface of the earth - where $g$ is a constant, finite value.
Argument:
For the sake of computation, let's pick some point on the surface of the infinite flat earth (henceforth "earth") to call the origin, and establish Euclidean coordinates $(x,y,z)$, with $x$ and $y$ spanning the surface of the earth and $z$ perpendicular to it. The density distribution $\rho$, is then given by
$$\rho(x,y,z)=\begin{cases}\bar\rho&-d\le z\le0\\0&\text{otherwise}\end{cases}$$
where $\bar\rho\approx 5.51\ \text g\cdot\text{cm}^{-3}$ is the average density of the earth and $d$ (value unknown) is its "depth." It should be noted that the assumption of uniform density on its own is fairly reasonable - it's a good enough approximation for many practical calculations in real life space science.
Now we have to correct a mistake. The author writes that the mass contained within a given region is equal to the density (units of mass per unit volume) multiplied by the area of that region. This is false. Since we are working in three dimensions, the product of density with area is a linear density, not mass. Instead, we have the mass given by the second equation in our second assumption, which leads to
$$m=\iiint_\Omega\rho\ dV=\iiint_{\Omega_0} 0\ dV+\iiint_{\Omega_{\bar\rho}}\bar \rho\ dV=\bar\rho V_{\bar\rho},$$
where $\Omega_0$ is the portion of $\Omega$ where the density is $0$, $\Omega_{\bar\rho}$ is the portion $\Omega$ where the density is $\bar\rho$, and $V_{\bar\rho}$ is the volume of that portion.
However, because the density within the earth is is assumed constant in every direction, we can replace the density with area density $\rho_A(x,y)=\bar\rho$ in our calculations. This kind of approximation - where we exploit some symmetry to "compress" everything into a two dimensional surface - is commonly used in describing the mechanics of sufficiently thin objects. Since the thickness of the earth relative to its presumed surface area is negligible, this is a reasonable simplification. Following this,
$$m=\iiint_\Omega \rho\ dV=\iint_{\pi_{xy}\Omega}\rho_A\ dA=\rho_A A$$
where $\pi_{xy}\Omega$ is the projection of $\Omega$ to the $xy$-plane. Note that in order for this to work, we have to choose $\Omega$ so that the cross sections parallel to the $xy$-plane are identical at every height.
Now consider the case where $\Omega$ is a [right circular] cylinder of radius $r$, whose axis passes through the origin at a right angle, of sufficient height and positioned so that the two circular faces are located above and below the upper and lower planes, respectively. Returning to Gauss's law, we have
$${\subset\!\supset}\llap\iint_{\partial\Omega}\mathbf g\cdot\mathbf n\ dA=-4\pi G\rho_AA,$$
where $A$ is the area of the projection of $\Omega$ to the $xy$-plane.
To compute the surface integral, we'll partition $\Omega$ into three parts - the upper circular face $U$, the lower circular face $L$, and the central "tube" which we'll call $C$.
$${\subset\!\supset}\llap\iint_{\partial\Omega}\mathbf g\cdot\mathbf n\ dA=\iint_U\mathbf g\cdot\mathbf n\ dA+\iint_L\mathbf g\cdot\mathbf n\ dA+\iint_C\mathbf g\cdot\mathbf n\ dA$$
While we don't know the value of $\mathbf g$ just yet, we know that the acceleration in any direction parallel to the plane is countebalanced by that in the opposite direction. This tells us 1. that $\mathbf g=(0,0,g^*)$ for some (possibly variable) $g^*$, at every point, 2. that the value of the third integral in the above equation is zero. Since the mass is distributed symmetrically to either side of the plane $z=-d/2$, we also know that the value of $\mathbf g$ at a given distance above the upper surface of the earth is the reflection of its value at the same distance below the lower surface along the $z$ direction. Since the unit normals of $U$ and $L$ are constant and the mass contained within $\Omega$ is independent of the height of $\Omega$, we can infer that $\lvert g^*\rvert=g$, for some constant value $g$, at all points outside of the earth. Thus,
$$\mathbf g(x,y,z)=\begin{cases}(0,0,\pm g)&z>0\\?&-d\le z\le 0\\(0,0,\mp g)&z<-d\end{cases}.$$
The $?$ is there because we neither know nor care what happens inside the earth. The $\pm$ is because we still don't know the direction of $\mathbf g$, only that it is parallel to the $z$-axis and flips when we cross the midplane. This doesn't actually matter, since the sign of $g$ is a matter of convention and we can always change it later. In any case,
$$\iint_U\mathbf g\cdot\mathbf n\ dA=\iint_L\mathbf g\cdot\mathbf n\ dA=\pm g A,$$
so
$${\subset\!\supset}\llap\iint_{\partial\Omega}\mathbf g\cdot\mathbf n\ dA=\pm2gA,$$
where $A$ is again the area of $\Omega$ projected to the $xy$-plane (equivalently, it is the area of $U$ or $L$.) Altogether this gets us
$$\pm 2gA=-4\pi G\rho_AA.$$
Solving for $\pm g$, we have $\pm g=-2\pi G\rho_A$. The gravitational constant and area density are conventionally positive, making the right-hand side of the equation negative. Adopting the sign convention of the author, we then have
$$g=2\pi G\rho_A.$$
This suffices to show that, given our premises, the gravitational acceleration at any point above or below the upper or lower surface of the earth, respectively, is indeed finite.
Note: For subsequent calculations, the author uses $\rho_A=\rho d$. This is a necessary correction if we want to pretend that the model is even plausible - any other means of determining the area density of the presumed infinite flat earth will reveal either the finite extent, significantly variable gravitation, asymmetric density, or nonplanar surface of the earth if attempted.
Answer 1: Ignoring the confusion of volume and area density, yes, the "proof" is valid - the premises entail the conclusion. However, this is predicated on the interpretation of the argument as a hypothetical, we cannot conclude that the premises or the conclusion are true.
Answer 1.5: It is unclear whether the author intends to evidence the possibility - in the absence of evidence to the contrary - of an infinite flat earth, or the agreement of the proposed model with the stated numerical values. Independently of its validity, the proof cannot serve as evidence for the latter, because our calculations only predict the relation between $g$, $d$, and $\bar\rho$, not their values.
The author uses the experimental values of $g$ and $\bar\rho$ to calculate the value of $d$. In order to check that this value agrees with observation, we need to measure $d$ as well. This is something that we can actually do, right now (and without directly contradicting the premises), using seismology. I desperately hope that someone does, and that, under the given assumptions, it comes out to $d\approx4195\ \text{km}$, because that would would be an absolutely hilarious coincidence.
Question 2: How is the [area] density of an infinite sheet defined?
The definition of area density (henceforth "density") at a point is given by its relation to mass. Specifically, given a region $\Omega$ bounded by a closed curve $\partial\Omega$,
$$m=\iint_\Omega\rho\ dA,$$
where $m$ is the mass contained in $\Omega$. If $A$ is the area of $\Omega$, then $\bar\rho=m/A$ gives the average density within $\Omega$. If the density is the same at every point in $\Omega$, then $\rho=\bar\rho$ at every point. You are correct that we cannot determine the average density of an infinite surface from its mass and area directly. However, if we assume that the density is constant, then the mass contained within any region divided by the area of that region will be the same as the density at every point on the surface. In the linked post, the author assumes that density is constant, so the value of $\rho$ could, in principle, be determined by measuring the mass contained within any finite chunk.
If the density is not constant at every point, then we would find it by looking at smaller and smaller regions. Let $\Omega_1,\Omega_2,\Omega_3,\ldots$ be a sequence of regions which all share a point $p$, such that $\Omega_n$ is contained within $\Omega_m$ wherever $m<n$, and $\lim_{n\to\infty} A_n=0$, where $A_n$ is the area of $\Omega_n$. Then the density at $p$ is $\rho=\lim_{n\to\infty}m_n/A_n$, where $m_n$ is the mass inside of $\Omega_n$. Physically, you can think of this as taking a big sheet of aluminum, weighing it and measuring its area, then removing a section from the center and repeating the process with that, until we're left with something so small we can't measure it. We extrapolate the density at the center of the original sheet from the sequence of densities we get by doing this.
As for the average density of an infinite surface, we get this by going in the opposite direction, looking at bigger and bigger regions. Let $\Omega_n$ be a sequence of regions such that $\Omega_m$ is contained inside of $\Omega_n$ whenever $m<n$, and the union of all regions covers the entire surface. Then $\bar\rho=\lim_{n\to\infty} m_n/A_n$, where $m_n$ and $A_n$ are again the mass and area, respectively, contained in $\Omega_n$.
