plotting - How to make a 3D topographic globe?

I've made some attempt with this dataset of elevation.

But I have some trouble with ListSurfacePlot3D, which can not show the globe correctly:

And I've checked that the data is not the problem, since the same data drawn by ListPointPlot3D shows a globe:

By the way, my goal is to make a topographic globe which shows 3D mountains and sea basins like this:

Here is my code:

(*elev data input*)
elev1d = BinaryReadList["D:\\topo\\ETOPO5.DAT", {"Integer16"}, ByteOrdering -> +1];
elev2d = ArrayReshape[elev1d, {2160, 4320}];
lati = Flatten @
   Transpose @ Table[Rest @ Table[i, {i, 90, -90, -1/12}], {4320}];

long = Flatten @ Table[Rest @ Table[i, {i, 0, 360, 1/12}], {2160}];
(* make a {lat, lon, altitude} matrix*)
elevlatlon = Transpose @ {lati, long, Flatten @ elev1d};
(*select part of the huge amount of data, add mean earth radius to altitude*)
elevlatlonInUse = (elevlatlon[[;; ;; 12, All]] /.
{m_, n_, o_} -> {m, n, o/200 + 6721}) /. {x_, y /; y > 180, z_} -> {x, y - 360, z};
coordsToXYZ[list_] := Transpose[{Cos[#[[1]]*Pi/180.]*Cos[#[[2]]*Pi/180.]*#[[3]], 
                                 Cos[#[[1]]*Pi/180.]*Sin[#[[2]]*Pi/180.]*#[[3]], 
                                 Sin[#[[1]]*Pi/180.]*#[[3]]} & @ Transpose[list]]
xyz = First[coordsToXYZ /@ {elevlatlonInUse}];

ListPointPlot3D[xyz, BoxRatios -> {1, 1, 1}]
ListSurfacePlot3D[xyz, BoxRatios -> {1, 1, 1}]

It's a little different from How to make a 3D globe?. That's a globe with a 2D texture covering it, but this is a real 3D globe with elevations shown in 3D as well.

P.S. Someone reminded me that, compared with the radius of the Earth (about $6371 \text{ km}$), even Mt. Everest ($8.8\text{ km}$) and the Marianas Trench ($-11\text{ km}$) can be ignored. That's true, I know, but to draw a globe with bumps, we can just scale the elevation. A visualized topographic globe is just for presentation, and not for calculation.

Answer

This answer is intended to demonstrate a neat method I'd recently learned for constructing interpolating functions over the sphere.

A persistent problem dogging a lot of interpolation methods on the sphere has been the subject of what to do at the poles. A recently studied method, dubbed the "double Fourier sphere method" in this paper (based on earlier work by Merilees) copes remarkably well. This is based on constructing a periodic extension/reflection of the data over at the poles, and then subjecting the resulting matrix to a low-rank approximation. The first reference gives a sophisticated method based on structured Gaussian elimination; in this answer, to keep things simple (at the expense of some slowness), I will use SVD instead.

As I noted in this Wolfram Community post, one can conveniently obtain elevation data for the Earth through GeoElevationData[]. Here is some elevation data with modest resolution (those with sufficient computing power might consider increasing the GeoZoomLevel setting):

gdm = Reverse[QuantityMagnitude[GeoElevationData["World", "Geodetic",

              GeoZoomLevel -> 2, UnitSystem -> "Metric"]]];

The DFS trick is remarkably simple:

gdmdfst = Join[gdm, Reverse[RotateLeft[gdm, {0, Length[gdm]}]]];

This yields a $1024\times 1024$ matrix. We now take its SVD:

{uv, s, vv} = SingularValueDecomposition[gdmdfst];

To construct the required low-rank approximations, we treat the left and right singular vectors (uv and vv) as interpolation data. Here is a routine for trigonometric fitting (code originally from here, but made slightly more convenient):

trigFit[data_?VectorQ, n : (_Integer?Positive | Automatic) : Automatic,

        {x_, x0_: 0, x1_}] :=
    Module[{c0, clist, cof, k, l, m, t},
           l = Quotient[Length[data] - 1, 2]; m = If[n === Automatic, l, Min[n, l]];
           cof = If[! VectorQ[data, InexactNumberQ], N[data], data];
           clist = Rest[cof]/2;
           cof = Prepend[{1, I}.{{1, 1}, {1, -1}}.{clist, Reverse[clist]}, First[cof]];
           cof = Fourier[cof, FourierParameters -> {-1, 1}];
           c0 = Chop[First[cof]]; clist = Rest[cof];
           cof = Chop[Take[{{1, 1}, {-1, 1}}.{clist, Reverse[clist]}, 2, m]];
           t = Rescale[x, {x0, x1}, {0, 2 π}];

           c0 + Total[MapThread[Dot, {cof, Transpose[Table[{Cos[k t], Sin[k t]},
                                                           {k, m}]]}]]]

Now, convert the singular vectors into trigonometric interpolants (and extract the singular values as well):

vals = Diagonal[s];
usc = trigFit[#, {φ, 2 π}] & /@ Transpose[uv];
vsc = trigFit[#, {θ, 2 π}] & /@ Transpose[vv];

Now, build the spherical interpolant, taking as many singular values and vectors as seen fit (I arbitrarily chose $\ell=768$, corresponding to $3/4$ of the singular values), and construct it as a compiled function for added efficiency:

l = 768; (* increase or decrease as needed *)

earthFun = With[{fun = Total[Take[vals, l] Take[usc, l] Take[vsc, l]]}, 
                Compile[{{θ, _Real}, {φ, _Real}}, fun, 
                        Parallelization -> True, RuntimeAttributes -> {Listable}, 
                        RuntimeOptions -> "Speed"]];

Now, for the plots. Here is an appropriate color gradient:

myGradient1 = Blend[{{-8000, RGBColor["#000000"]}, {-7000, RGBColor["#141E35"]},
                     {-6000, RGBColor["#263C6A"]}, {-5000, RGBColor["#2E5085"]},
                     {-4000, RGBColor["#3563A0"]}, {-3000, RGBColor["#4897D3"]},
                     {-2000, RGBColor["#5AB9E9"]}, {-1000, RGBColor["#8DD2EF"]},

                     {0, RGBColor["#F5FFFF"]}, {0, RGBColor["#699885"]},
                     {50, RGBColor["#76A992"]}, {200, RGBColor["#83B59B"]},
                     {600, RGBColor["#A5C0A7"]}, {1000, RGBColor["#D3C9B3"]},
                     {2000, RGBColor["#D4B8A4"]}, {3000, RGBColor["#DCDCDC"]},
                     {5000, RGBColor["#EEEEEE"]}, {6000, RGBColor["#F6F7F6"]},
                     {7000, RGBColor["#FAFAFA"]}, {8000, RGBColor["#FFFFFF"]}}, #] &;

Let's start with a density plot:

DensityPlot[earthFun[θ, φ], {θ, 0, 2 π}, {φ, 0, π},
            AspectRatio -> Automatic, ColorFunction -> myGradient1,

            ColorFunctionScaling -> False, Frame -> False, PlotPoints -> 185, 
            PlotRange -> All]

Due to the large amount of terms, the plotting is a bit slow, even with the compilation. One might consider using e.g. the Goertzel-Reinsch algorithm for added efficiency, which I leave to the interested reader to try out.

For comparison, here are plots constructed from approximations of even lower rank ($\ell=128,256,512$), compared with a ListDensityPlot[] of the raw data (bottom right):

Finally, we can look at an actual globe:

With[{s = 2*^5},
     ParametricPlot3D[(1 + earthFun[θ, φ]/s)

                      {Sin[φ] Cos[θ], Sin[φ] Sin[θ], -Cos[φ]} // Evaluate,
                      {θ, 0, 2 π}, {φ, 0, π}, Axes -> None, Boxed -> False, 
                      ColorFunction -> (With[{r = Norm[{#1, #2, #3}]}, 
                                              myGradient1[s r - s]] &),
                      ColorFunctionScaling -> False, MaxRecursion -> 1,
                      Mesh -> False, PlotPoints -> {500, 250}]] // Quiet

(I had chosen the scaling factor s to make the depressions and elevations slightly more prominent, just like in my Community post.)

Of course, using all the singular values and vectors will result in an interpolation of the data (tho it is even more expensive to evaluate). It is remarkable, however, that even the low-rank DFS approximations already do pretty well.