Difficulty interpolating near poles using lat/lons

I've been writing a 2D interpolation scheme for use on the globe which utilizes lat/lon values. Once I try to use values closer to the poles, latitude distortion becomes an issue. For example:

270,-86 ---------- 0, -90          This looks like an ordinary square on the globe
|        |                  but when these values are projected to a plane,
|        |                  they create a shape that's unworkable.
|        |                  Note that opposite corners have the same latitude.
225,-84 ---------- 180,-86

I've been converting from spherical to cartesian coordinates in an attempt to create a 2D plane but there's still the issue of these planes existing in three dimensions and I've had difficulty "flattening" them. Surely there must be a common method for dealing with the latitude distortion near the poles but I don't know what it is. Any help would be much appreciated.

Thanks.

• Which projection are you trying to use? Also, are you working with the entire Earth or just the section in your question? Aug 14 '20 at 19:02
• Cross-site duplicate: gis.stackexchange.com/q/429/136572
– Rob
Aug 14 '20 at 19:46
• Welcome. Yep, ask the GIS guys as well. It would be good to know the tools you're using and what you're up to. How about a cube map projection that projects each ray from the center of an ellipsoid onto one of the six faces (there's the 2D) of a cube ? No singularities at the poles and easy to make handy portions for big planets. marlam.de/ecm with reference implementation. Btw., interpolation is done automatically when skilfully using the rendering and rasterization pipeline of a graphics API ...
– user20217
Aug 15 '20 at 10:05
• Is this assuming a spherical Earth? It looks like that assumption is made, but when I answered accordingly I got critical comments because the spherical Earth is only an approximation. Aug 16 '20 at 14:33
• Positve criticism, of course :-) Converting cartesian to geodetic is an iterative algorithm, needs surface normal calculation and height offset and some trigonometry. A sphere is muuuch simpler, but not the shape of the earth, depending on what's the goal. Georeferenced data won't fit a sphere. I think we should be able to judge if the problem is data driven or methodological.
– user20217
Aug 16 '20 at 14:39

$$\color{blue}{\text{The answer here presumes that Earth is a sphere.}}$$

$$\color{blue}{\text{This is not fully true, but the question seems to build in that assumption.}}$$

$$\color{blue}{\text{That should be clarified to determine if this answer is relevant.}}$$

Do not use a 2D-plane which does not generally exist for four points. Use a tetrahedron. You can interpolate inside the tetrahedron and then project the interpolated point (which will be inside the globe) radially onto the globe's surface.

Say you pick the South Pole, S, and three points A, B, C at 60° south latitude and longitude 0°, 90° east, 90° west. Converting these coordinates to Cartesian form with the South Pole at $$(0,0,-1)$$ prime meridian passing through $$(1,0,0)$$, we get the coordinates

South Pole = $$(0,0,-1)$$

A = $$(1/2,0,-\sqrt3/2)$$

B = $$(0,1/2,-\sqrt3/2)$$

C = $$(0,-1/2,-\sqrt3/2)$$

Now suppose you want to interpolate with a weight of 1/2 for S and 1/6 for each of the other three points. Take the linear combinations of the above coordinates with those coefficients to get a point P inside the tetrahedron:

$$P=(1/2)(0,0,-1)+(1/6)(1/2,0,-\sqrt3/2)+(1/6)(0,1/2,-\sqrt3/2)+(1/6)(0,-1/2,-\sqrt3/2)=(1/12,0,-(2+\sqrt3)/4)$$

Point P has the following physical significance: if you draw the smaller tetrahedron PABC, opposite S, then this tetrahedron has 1/2 the volume of the big tetrahedron ABCS, where 1/2 is the coefficient you put in for S. For this choice of P the corresponding tetrahedra opposite A,B,C will each have 1/6 the volume of the big one.

Now we need to project this point onto the surface ofthe globe. Work out its distance $$d$$ to the origin:

$$d^2=(1/12)^2+0^2+((2+\sqrt3)/4)^2=(16+9\sqrt3)/36\approx 0.8775, d\approx 0.9367$$

Note that $$d<1$$. Divide this $$d$$ into the coordinates of P getting the approximate result

$$P'\approx (0.0890, 0, -0.9960)$$

Converting back to spherical coordinates then gives a latitude of $$\approx 84.90°=84°54'$$ south and longitude zero, the latter due to equally weighting the eastern and western hemispheres in this example. The interpolated point is more than halfway from the 60° south latitude circle to the South Pole because we interpolated from a distribution of points along this circle rather than a single point; the centroid of of this distribution ( $$\triangle ABC$$) is at more than 60° south latitide.

• Then I need some clarification to the question. The OP seemed to me to be eendering the Earth as a sphere. I will ask for clarification in a comment to the question. Aug 16 '20 at 14:26