Моделируем поверхность Земли в пару строчек

Вы когда-нибудь играли в 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 координаты это сферические углы

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

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

См. UPD в конце статьи

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

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

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

UPD 10.01.2025

Как упомянули в комментариях такой способ проекции облаков очень наивный с очевидной особенностью на полюсах.

Хотелось быстро решить проблему, но оказалось, что в ней целая наука.

В целом нужно решить две проблемы

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

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;

В повторяемости можно убедиться

Table[cloudTexture, {2}, {2}] // Grid 

• преобразовать либо сетку, либо исходный шум, чтобы избежать искажений на полюсах

Я сразу упомяну, что не являюсь экспертом в области геодезии, а лишь «слышал» о некоторых концептах из вычислительной геометрии и 3D графики 😉

В любом случае мы будет «натягивать» эту текстуру в том или ином виде (в виде полигонов) на сферу. Это означает мы преобразуем все ее точки следуя одной из картографических проекций. В 3D графике есть распространенный способ UV-mapping, что является аналогом равнопромежуточной картографической проекции с точностью до коэффициентов, отвечающих за сдвиги углов.

При условии, что наши исходные координаты uv нормированы в пределах от 0 до 1

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

Представим, если мы совершим некое обратное преобразование над нашей исходной текстуре, то при проецировании ее на сферу они должны компенсировать друг-друга и мы увидим неискаженное изображение облаков. Тогда вопросы:

• Чем является сейчас наша исходная текстура и что мы хотим увидеть в конце?

• Что является неискаженным изображением?

• Как мы в целом можем говорить о том, что мы хотим увидеть на сфере глядя на прямоугольное изображение шума Перлина в декартовой системе координат?

На эти вопросы можно было бы ответить в этой статье, если бы я был экспертом в поставленной теме. К сожаления, я не могу удовлетворять этому критерию, поэтому хотя бы поищем путь избавления от артефактов на полюсах. Глядя на текстуры из игр можно увидеть закономерность

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

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

Приведение проекции Мольвейде к равнопромежуточной картографической проекции

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]};  cloudTexture  newCloudTexture = ImageForwardTransformation[   cloudTexture,   mollweidetoequirect,   DataRange -> {{-2 Sqrt[2], 2 Sqrt[2]}, {-Sqrt[2], Sqrt[2]}},   PlotRange -> All ]

Остается повторить Marching Cubes на ней и преобразовать все вершины с помощью формулы выше

With[{plain = ImageData[newCloudTexture]}, Table[plain Exp[-( i)^2/200.], {i, -20,20}]];  {vertices, normals} = CMarchingCubes[%, 0.2, "CalculateNormals" -> False];  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 = Map[Function[v,     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]]}     ]   ] , vertices];  {   clouds = GraphicsComplex[0.017 pvertices, Polygon[1, Length[vertices]]] } // Graphics3D

Теперь собираем все вместе

Бонус: нормали облаков

Плоское затенение это может быть не так красиво. Давай-те включим расчёт нормалей и и спроецируем их также на сферу

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]];  {pverticesn, 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

Если в GraphicsComplex задано явно VertexNormals, тогда при затенении поверхность будет интерполироваться по всем трем нормалям треугольника, вместо одной (рассчитывается автоматически графической библиотекой если других данных нет, тогда полигон выглядит как плоскость aka эффект из демосцен 80х).

И теперь соберем финальное изображение