How would you go about finding the real area of a polygon which is defined by latitude/longitude points? Remember that the real area is different as a map is distorted towards the poles.
[Math] Area of a geopolygon or polygon defined by longitude/latitude points
geometryspherical-geometry
Related Solutions
Since you only want to calculate the number roughly and the distance of a mile is small compared to the Earth's radius, it should be sufficient to use a linear approximation with the surface element of the transformation to Cartesian coordinates. For a rough calculation, we can take the Earth as spherical and assume that the latitude and longitude are spherical coordinates $\theta$ and $\phi$, respectively. The surface element for spherical coordinates is $r^2\sin\theta\mathrm d\theta\mathrm d\phi$. Loosely speaking, that means that an infinitesimal area $\mathrm d\theta\mathrm d\phi$ in the coordinate space corresponds to an infinitesimal surface area $r^2\sin\theta\mathrm d\theta\mathrm d\phi$ on the actual surface of the sphere. For your purposes, we can take $\sin\theta$ as constant, in which case the areas are simply proportional to each other, with proportionality constant $r^2\sin\theta$. For your purposes, we can calculate the area of the circle as if it were flat, so a circle with a radius of one mile has an area $A$ of $\pi$ square miles. The radius $r$ of the Earth is roughly $4000$ miles, and at latitude $\theta=41.174104^\circ$ we have $\sin\theta\approx0.66$, so the corresponding area in coordinate space is $A/(r^2\sin\theta)\approx\pi/(4000^2\cdot0.66)\approx3\cdot10^{-7}$. The area per grid point in a grid of coordinates with $d$ decimals is $(10^{-d})^2=10^{-2d}$, so for instance if you use $d=6$ digits you'd have roughly $3\cdot10^{-7}/10^{-12}=3\cdot10^5$ different coordinates within a one-mile radius, whereas if you use $d=3$ digits the result is $3\cdot10^{-7}/10^{-6}=0.3$, so you wouldn't expect any other coordinates within that radius.
Given the latitudes of two points, $\beta_1$ and $\beta_2$, and their difference in longitude, $\Delta\lambda=\lambda_2-\lambda_1$, compute $\delta\in[0,\pi]$ with $$ \cos(\delta)=\sin(\beta_1)\sin(\beta_2)+\cos(\beta_1)\cos(\beta_2)\cos(\Delta\lambda)\tag{1} $$ Then compute $\alpha$ using $$ \tan(\alpha/2)=\frac{\sin(\Delta\lambda)\cos(\beta_1)\cos(\beta_2)}{\sin{\beta_2}-\sin(\delta-\beta_1)}\tag{2} $$
After computing $\delta$ and $\alpha$, compute $\beta'\in[-\frac{\pi}{2},\frac{\pi}{2}]$ with $$ \sin(\beta')=\sin(\beta_1)\cos(k\delta/5)+\cos(\beta_1)\sin(k\delta/5)\cos(\alpha)\tag{3} $$ Then compute $\Delta\lambda'=\lambda'-\lambda_1$ using $$ \tan(\Delta\lambda'/2)=\frac{\sin(k\delta/5)\sin(\alpha)\cos(\beta_1)}{\cos(k\delta/5)+\cos(\beta_1+\beta')}\tag{4} $$ Where $k$ in $(3)$ and $(4)$ ranges in $\{1,2,3,4\}$.
Best Answer
The area of a spherical polygon is proportional to its excess angle. In the plane, a $n$-gon would have an interior angle sum of $(n-2)\pi$. On the sphere, the angle sum $\alpha$ is larger than that. Compute the difference, multiply with a suitable constant (which is $r^2$ as I'll show below) and you have the area.
$$A = \bigl(\alpha-(n-2)\pi\bigr)r^2$$
The fact that the proportionality factor has to be $r^2$ can be deduced like this: On the unit sphere, the hemisphere can be interpreted as a $2$-gon with two angles of $\pi$, stretching from pole to pole. Therefore you get $A=2\pi$, which is exactly half of the surface area of the complete sphere, as expected. Scale your sphere by a factor of $r$, and all areas will scale by $r^2$.
To compute these interior angles from latitude and longitude, you can have a look at how to compute the initial course for great-circle navigation. The difference between bearings towards adjacent corners will give you the angle at the current corner.