Вы когда-нибудь играли в Outer Wilds? Планеты там невероятно красивы. Это собственно стало основной мотивацией создать свою простую модель планеты, используя реальные географические данные о высотах и немножко магии Wolfram Language

Как разбросать точки на сфере

Очевидно, даже если мы стащим карту высот земли возникнет проблема проекции ее на сферу. Это очень сложная тема с триангуляцией (ведь потом надо сделать из этого полигоны) и прочим, про нее довольно много написано. К примеру

• https://www.redblobgames.com/x/1842-delaunay-voronoi-sphere/

• https://medium.com/@oscarsc/four-ways-to-create-a-mesh-for-a-sphere-d7956b825db4

Как один из вариантов - мы пойдем от обратного - сначала создадим пробные точки на сфере, а затем уже запросим по ним из гео-данных значения высот.

Итак, нам нужно равномерно распределить N точек по сфере. Не нужно изобретать велосипед, если есть богатая стандартная библиотека

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

Географические данные

К счастью или сожалению, в стандартной библиотеке WL уже лежит грубая карта всей Земли (на всякий случай). Если нужно поточнее, она пойдет в интернет и добудет больше. Собственно, к делу

GeoElevationData[GeoPosition[Here]]

Quantity[490, "Meters"]

здесь Here возвращает текущую широту и долготу в градусах. В функцию можно передать не одно значение, но и целый список. Собственно этот список мы возьмем из точек на сфере

points = SpherePoints[5000];
latLon = (180/Pi) {90Degree - ArcCos[#[[3]]], ArcTan[#[[1]], #[[2]]]} &/@ points;

elevation = GeoElevationData[GeoPosition[latLon]];
elevation = (elevation + Min[elevation])/Max[elevation];

Здесь мы переходим из декартовой системы координат в сферическую (геодезическую точнее)

получаем высоты и нормируем их.

Остается связать их с исходными точками на сфере, используя нормированную высоту как расстояние по нормали

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

здесь мы масштабируем их на глаз, чтобы высота над уровнем моря лишь слегка "модулировала" поверхность сферы, таким образом мы получим видимый рельеф Земли

rainbow = ColorData["DarkRainbow"];

ListSurfacePlot3D[
  surface, 
  Mesh->None, MaxPlotPoints->100, 
  ColorFunction -> Function[{x,y,z}, rainbow[1.5(2 Norm[{x,y,z}]-1)]], 
  ColorFunctionScaling -> False
]

Генерируем облака

Какая же Земля без облаков? Какие же облака без Шума Перлина? Следующий кусок я честно украл на одном из форумов (никаких GPT вам!)

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

Остается другая проблема - как сделать их трехмерными? Я не придумал ничего лучше, как использовать технику Marching Cubes и сгенерить low-poly облака. Однако для начала "растянем" двумерное изображение в трехмерное с затуханием по краям

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

Для того, чтобы превратить это скалярное поле в полигоны можно воспользоваться функцией ImageMesh, однако с ней довольно тяжело работать, если требуется применить нелинейные преобразование к этой сетке. А нам придется их сделать, иначе, как потом натянуть этот квадрат на сферу?

Воспользуемся внешней библиотекой wl-marchingcubes

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

<<JerryI`MarchingCubes`

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

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

Результат записывается в vertices. Библиотека генерирует по умолчанию неиндексированный набор вершин треугольников. Это означает, что можно напрямую интерпретировать их как

Polygon[{
  {1,2,3}, //треугольник 1
  {4,5,6}, //треугольник 2
  ...
}]

где ни одна из вершин не переиспользуется другим треугольником. Такой формат особенно прост для GPU, так как нужно отправить лишь один такой список в один из буферов WebGL. К счастью примитив Polygon поддерживает такой вариант

GraphicsComplex[vertices, Polygon[1, Length[vertices]]] // Graphics3D

Здесь я намеренно отключил генерацию нормалей, так как мне не удалось привести их к сферической системе координат без видимых искажений.

Натягиваем на сферу

Здесь мы сталкиваемся с такой же проблемой, что и ранее. Как натянуть это на сферу? Для иллюстративных целей попробуем Такое преобразование в лоб, как будто наши x,y координаты это сферические углы

где . Вопрос зачем здесь некий угол сдвига - так как такое преобразование создает серьезные артефакты на полюсах сферы, можно попытаться косметически поправить их сдвинув полюс в область, где меньше геометрии.

Если кто-то из хабровчан знает способ получше - прошу в комментарии ?

Итак, попробуем этот вариант

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

Собираем все вместе

Можно расширить исходный SurfacePlot опциями и добавить через них дополнительные графические примитивы. Ах да, еще отключить стандартное освещение. Мы же в космосе

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, 
  Lighting->None,
  PlotStyle->Directive["Shadows"->True, "CastShadow"->True],

  Prolog -> {
    Directive["Shadows"->True, "CastShadow"->True],
    clouds,
    HemisphereLight[LightBlue, Orange // Darker // Darker],
    SpotLight[Orange, {-2.4909, 4.069, 3.024}]
  },
  
  Background->Black
, ImageSize->800]

Можно поиграться с освещением в реальном времени, добавив обработчик на свойство transform мелкой сферы на месте источника света

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


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->None,
  PlotStyle->Directive["Shadows"->True, "CastShadow"->True],

  Prolog -> {
    Directive["Shadows"->True, "CastShadow"->True],
    clouds,
    HemisphereLight[LightBlue, Orange // Darker // Darker],
    SpotLight[Orange, lightPos // Offload],
    EventHandler[Sphere[lightPos, 0.001], {"transform" -> Function[t,
      lightPos = t["position"]
    ]}]
  },
  
  Background->Black
, ImageSize->800]

Анимация

Добавим простенькую анимацию движения источника света, а также вращения сферы облаков, чтобы разнообразить сцену

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
]

Чтобы запустить движение нужно "пнуть" обработчик, тогда обновление позиции источника света вызовет прерывание, которое снова обновит позицию светового тела и так далее синхронно с циклом отрисовки браузера

Посмотреть в живую

Среда разработки

Все примеры были показаны в открытой среде-разработки-блокноте WLJS Notebook

Спасибо за внимание ??‍♂️