Skip to main content

Modeling Earth

· 2 min read
Kirill Vasin

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

Have you ever played Outer Wilds ? Planets there are so beautiful. That was the main motivation to rebuild one using real elevation data and a bunch of WL's magic

Download original notebook

Spreading points on a sphere​

A small task to do

Graphics3D[{
  Sphere[],
  Red, PointSize[0.01], SpherePoints[1000]//Point
}]

(*VB[*)(FrontEndRef["6133ae90-57e8-4be5-8ed4-54d14e466715"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKmxkaGyemWhrompqnWuiaJKWa6lqkppjompqkGJqkmpiZmRuaAgB9OhUh"*)(*]VB*)

Convert to Geo coordinates​

Find latitude and langitude

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];

In the case if Geo did not work. Here is a dump

{points, latLon, elevation} = (*VB[*)(Get[FileNameJoin[{".iconized", "iconized-337e.wl"}]])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeH5DwTM7Py6xKdcvMSXXKr8j0YOVmAADrbgrN"*)(*]VB*);
rainbow = ColorData["DarkRainbow"];

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]
(*VB[*)(FrontEndRef["ef4101c8-35c1-415b-a881-83c0e3a8f8f1"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKp6aZGBoYJlvoGpsmG+qaGJom6SZaWBjqWhgnG6QaJ1qkWaQZAgB+RxVt"*)(*]VB*)

That looks like Earth!

Generating clouds​

One can simply use Perlin noise and perform marching cubes

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

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

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

p0 = p = Sqrt[k2];

Image[im[p0 += p]]
(*VB[*)(FrontEndRef["f77ace65-5c0e-4968-abdf-b65fdb70367f"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKp5mbJyanmpnqmiYbpOqaWJpZ6CYmpaTpJpmZpqUkmRsYm5mnAQCRdxZG"*)(*]VB*)

It is hard to wrap on the sphere RegionMesh function from the standard library. Therefore we will go with an external library

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

<<JerryI`MarchingCubes`

Now generate vertices and map them to sphere

With[{plain = im[p0+=p]}, Table[plain Exp[-( i)^2/200.], {i, -20,20}]];

{vertices, normals} = CMarchingCubes[%, 0.2, "CalculateNormals" -> False];

vertices = Map[Function[v,
    With[{\[Rho] = 50.0 + 0.25 (v[[3]] - 10), \[Phi] = 2.0 Pi v[[1]]/127.0, \[Theta] =  Pi/2 + Pi v[[2]]/127.0},
      {\[Rho] Cos[\[Phi]] Cos[\[Theta]], \[Rho] Sin[\[Phi]] Cos[\[Theta]], \[Rho] Sin[\[Theta]]}
    ]
  ]
, vertices];

{
  clouds = GraphicsComplex[0.017 vertices, Polygon[1, Length[vertices]]]
} // Graphics3D;

We shifted θ\theta angle for purpose, to avoid visibles artifacts on poles.

Combine​

Plot together

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

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

animation = CreateUUID[];
EventHandler[animation, Function[Null,
  lightPos = RotationMatrix[1 Degree, {1,1,1}].lightPos;
  rotationMatrix = RotationMatrix[angle, {0,0,1}];
  angle += 0.5 Degree;
]];

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]
  },

  Epilog -> AnimationFrameListener[lightPos // Offload, "Event"->animation],
  Background->Black
]
(*VB[*)(FrontEndRef["c8e48c06-39cf-41b3-b9e8-9d005bb356f5"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJ1ukmlgkG5jpGlsmp+maGCYZ6yZZplroWqYYGJgmJRmbmqWZAgCHGhWu"*)(*]VB*)

To start animation

EventFire[animation, True];