# Function calculating the curvature of earth

I am interested in calculating a box, based on a center point and its ∆ values. The box is an excerpt of the earth' surface. Note that I am not working with a projection of the earth (2D), but with values accounting for its original form - being a sphere (3D).

Given is the center (c) point - in geographical coordinates (latitude (x) / longitude (y)) -> c(x/y). Another information I am having is ∆ latitude (∆y) and ∆ longitude (∆x).

Here an abstraction of the earth and how the points are located.

       _______90°______
/                \
∆(y)               \
/                    \
c(x/y) --- ∆(x)--------\
/                        \ Longitude
/                          \
|--------0° Equator --------|
\                         /
\                       /
\----- latitude ------/
\                   /
\                 /
------ -90°----


Given this information, I can calculate x/y from p1..p4, i.e. p2 = x + ∆x or p4 = x - ∆x. Once I have p1...p4, I have the coordinates for x1...x4 which are the corners.

x1------p1------x2
|       |       |
|       ∆y      |
|       |       |
p4--∆x--c(x/y)--p2
|       |       |
|       |       |
|       |       |
x4------p3------x3


When I draw the polygon created from x1...x4, it covers more than it should. This gets more extreme the closer I am to north - or south pole. I assume this is because of the curvature of the earth. The more I am moving north or south, the more distance is between two latitude lines

How would I correct the calculated polygon regarding the curvature of the earth? I tried a linear function f(x) = 1/90 * x but that only approximates it roughly.

• The "geographical coordinates" of course behave very badly near the poles. Why do you want to use them? From the mathematical standpoint, in order to get reasonable results you should use several different charts to cover the surface of the earth. If you like "latitudes and longitudes" consider a "multipolar world" where in addition to the north and south poles, you have few more (at least two more). – Moishe Kohan Feb 16 '18 at 13:49
• @MoisheCohen I do not have choice. The only information I am getting are geographical coordinates... I know its not good - but that is the only information I am having. Could you explain the chart approach? – Stophface Feb 16 '18 at 13:50
• You mean that for the center point $c$ you have the "standard" coordinates, but that's OK. My suggestion is, as I said, to use several different charts. You would have to compute the "transition maps" near your center point $c$. You thus, will have, say, two different "boxes" centered at $c$: One is standard and the other is obtained from the first one by a rotation. In a band which avoids the standard poles you will use one box (the standard one) while near the "standard" poles you will use another one. The more precision you need, the more boxes/charts you would have to use. – Moishe Kohan Feb 16 '18 at 13:56

$$dS=r^2\sin\phi\, dr\,d\phi\,d\theta$$ so that with a constant $\theta$ step, the $\phi$ step should grow like $\csc\phi$. The function that achieves such a behavior is the antiderivative,
$$\int\dfrac{d\phi}{\sin\phi}=-\log(\cot(\phi)+\csc(\phi)).$$
This function has vertical asymptotes at $\phi=0$ and $\phi=\pi$, which was to be expected as it is impossible to achieve constant area at the poles. A possible solution is to use plain disks around the poles, locally dropping the mesh topology.