Function calculating the curvature of earth

Question

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.

Answer 1

You don't state it very clearly, but I guess that you want to mesh a sphere with a grid, in such a way that the tiles are approximately of the same area (correct me if I am wrong).

From differential geometry, the element of area in spherical coordinates is

$$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.

Lookup "A general rule for disk and hemisphere partition into equal-area cells".