Modelling Earth

*using GeoData, surface plot, and a bunch of other beautiful stuff*

Have you ever played Outer Wilds? The planets there are incredibly beautiful. This actually became the main motivation to create my own simple model of a planet, using real geographical elevation data and a bit of Wolfram Language magic.

Preparations

Get elevation data

points = SpherePoints[5000];
latLon = (*FB[*)((180)(*,*)/(*,*)(Pi))(*]FB*) {90Degree - ArcCos[#[[3]]], ArcTan[#[[1]], #[[2]]]} &/@ points;
elevation = GeoElevationData[GeoPosition[latLon]];
elevation = (elevation + Min[elevation])/Max[elevation];

Make surface

surface = MapThread[(#1 (0.8 + 0.1 #2))&, {points, elevation}];

Generate cloud texture

n = 256;
ratio = 1/2;

k2 = Outer[Plus, #, #] &[RotateRight[N@Range[-n, n - 1, 2]/n, n/2]^2];

k2 = k2[[;; ;; 1/ratio, All]];

spectrum = With[{d := RandomReal[NormalDistribution[], {n ratio, n}]},
   (1/n) (d)/(0.001 + k2)]; 
spectrum[[1, 1]] *= 0;

im[p_] := Clip[Re[InverseFourier[spectrum Exp[I p]]], {0, ∞}]^0.5

p0 = p = Sqrt[k2];

cloudTexture = im[p0 += p] // Image

Transform for UV mapping

lat[y_, rad_:1] := ArcSin[(2 theta[y, rad] + Sin[2 theta[y, rad]])/Pi];
lon[x_, y_, rad_:1, lonzero_: 0] := lonzero + Pi x/(2 rad Sqrt[2] Cos[theta[y, rad]]);
theta[y_, rad_:1] := ArcSin[y/(rad Sqrt[2])];
mollweidetoequirect[{x_, y_}] := {lon[x, y, 1], lat[y]};

newCloudTexture = ImageForwardTransformation[
  cloudTexture,
  mollweidetoequirect,
  DataRange -> {{-2 Sqrt[2], 2 Sqrt[2]}, {-Sqrt[2], Sqrt[2]}},
  PlotRange -> All
];

Image[%, Magnification->1]

Generate polygons

PacletRepositories[{
    Github -> "https://github.com/JerryI/wl-marching-cubes" -> "master"
}]

<<JerryI`MarchingCubes`

With[{plain = ImageData[newCloudTexture]}, Table[plain Exp[-( i)^2/200.], {i, -20,20}]];

{vertices, normals} = CMarchingCubes[%, 0.2];

vertices = With[{
  offset = {Min[vertices[[All,1]]], Min[vertices[[All,2]]], 0},
  maxX = Max[vertices[[All,1]]] - Min[vertices[[All,1]]],
  maxY = Max[vertices[[All,2]]] - Min[vertices[[All,2]]]
}, Map[Function[v, With[{p = v - offset}, {p[[1]]/maxX, p[[2]]/maxY, p[[3]]}]], vertices]];

{pvertices, pnormals} = MapThread[Function[{v,n},
    With[{\[Rho] = 50.0 + 0.25 (v[[3]] - 10), \[Phi] =  2 Pi Clip[v[[1]], {0,1}], \[Theta] =  Pi Clip[v[[2]], {0,1}]},
      {
        {\[Rho] Cos[\[Phi]] Sin[\[Theta]], \[Rho] Sin[\[Phi]] Sin[\[Theta]], \[Rho] Cos[\[Theta]]},

        n[[3]] {Cos[\[Phi]] Sin[\[Theta]], Sin[\[Phi]] Sin[\[Theta]], Cos[\[Theta]]} + 
        n[[2]] {Cos[\[Theta]] Cos[\[Phi]],Cos[\[Theta]] Sin[\[Phi]],-Sin[\[Theta]]} + 
        n[[1]] {-Sin[\[Theta]] Sin[\[Phi]],Cos[\[Phi]] Sin[\[Theta]],0}
      }
    ]
  ]
, {vertices, -normals}] // Transpose;

{
  clouds = GraphicsComplex[0.017 pvertices, Polygon[1, Length[vertices]], VertexNormals->pnormals]
} // Graphics3D

Final composition

Put all together and animate

rainbow = ColorData["DarkRainbow"];
lightPos = {-2.4909, 4.069, 3.024};

rotationMatrix = RotationMatrix[0., {0,0,1}];

EventHandler[InputRange[0,1, 0.01], Function[value,
  lightPos = RotationMatrix[2Pi value, {1,1,1}].{-2.4909, 4.069, 3.024} // N;
  rotationMatrix = RotationMatrix[2Pi value, {0,0,1}] // N;
]]

ListSurfacePlot3D[
  MapThread[(#1 (0.8 + 0.1 #2))&, {points, elevation}], 
  Mesh->None, MaxPlotPoints->100, 
  ColorFunction -> Function[{x,y,z}, rainbow[1.5(2 Norm[{x,y,z}]-1)]], 
  ColorFunctionScaling -> False, 
  Lighting->"Default",
  PlotStyle->Directive["Shadows"->True, "CastShadow"->True],

  Prolog -> {
    Directive["Shadows"->True, "CastShadow"->True],
    GeometricTransformation[clouds, rotationMatrix // Offload],
    HemisphereLight[LightBlue, Orange // Darker // Darker],
    SpotLight[Orange, lightPos // Offload]
  },
  
  Background->Black,
  ImageSize->800
]