Procedural generation of global cloud cover

This article is about generating realistic planetary cloud cover. In an initial attempt I applied random rotation fields to a sphere with Worley noise and posted it on Reddit. However the results did not look convincing. A helpful comment by mr_bitshift pointed out the publication Bridson et al. “Curl noise for procedural fluid flow”. Also in a response to a later Reddit post smcameron shared an impressive result of a gas giant generated using curl noise.

In two dimensions curl noise is fairly easy to understand and implement. For a thorough description of 2D curl noise see Keith Peters’ article “Curl noise, demystified”. Basically one starts with a potential field such as multiple octaves of Worley noise. One then extracts the 2D gradient vectors and rotates them by 90°.

To generate curl vectors for a spherical surface one can use 3D Worley noise and sample the gradients on the surface of the sphere. The gradient vectors then need to be projected onto the sphere. This can be achieved by projecting the gradient vector onto the local normal vector of the sphere using vector projection. By subtracting the projected vector from the gradient vector one obtains the tangential component of the gradient vector.

The resulting vector p needs to be rotated around the normal n by 90°. This can be achieved by rotating the vector p into a TBN system, rotating by 90° around N and then transforming back. The GLSL functions for the rotation (without OpenGL tests) are shown below:

#version 410 core

vec3 orthogonal_vector(vec3 n)
{
  vec3 b;
  if (abs(n.x) <= abs(n.y)) {
    if (abs(n.x) <= abs(n.z))
      b = vec3(1, 0, 0);
    else
      b = vec3(0, 0, 1);
  } else {
    if (abs(n.y) <= abs(n.z))
      b = vec3(0, 1, 0);
    else
      b = vec3(0, 0, 1);
  };
  return normalize(cross(n, b));
}

mat3 oriented_matrix(vec3 n)
{
  vec3 o1 = orthogonal_vector(n);
  vec3 o2 = cross(n, o1);
  return transpose(mat3(n, o1, o2));
}

// note that axis needs to be a unit vector
vec3 rotate_vector(vec3 axis, vec3 v, float cos_angle, float sin_angle)
{
  mat3 orientation = oriented_matrix(axis);
  vec3 oriented = orientation * v;
  mat2 rotation = mat2(cos_angle, sin_angle, -sin_angle, cos_angle);
  vec3 rotated = vec3(oriented.x, rotation * oriented.yz);
  return transpose(orientation) * rotated;
}

In OpenGL one can create a cubemap where each pixel on each surface contains a 3D warp vector.

• Using a fragment shader the cubemap is initialised to be an identity transform for unit vectors.
• A second fragment shader is used to initialise a cubemap with the curl vectors which are tangential to the sphere.
• A third fragment shader is called multiple times to renormalize and increment the identity transform to become a warp field.
• A final fragment shader uses the cubemap warp field to perform lookups in a 3D Worley noise field to generate a cubemap of the global cloud cover.

If one uses octaves of Worley noise one obtains vortices rotating in one direction. To obtain prevailing winds and vortices with different direction of rotation depending on the latitude one can use the function (1+sin(2.5*latitude))/2 to mix positive and negative Worley noise.

Below is a result obtained using the method described in this article.

Also see here for a video.

See cover.clj for source code.

Enjoy!

Update

Another detail I forgot to mention is that the fragment shaders and the cubemap texture lookups use modified vectors to avoid performing lookups in the texture clamping regions which would lead to seams in the cloud cover. I.e. when converting fragment coordinates, one increases the range of the index by half a pixel on both ends:

vec2 x = (gl_FragCoord.xy - 0.5) / (size - 1);

Furthermore when performing lookups, two coordinates of the lookup vector are scaled down by half a pixel:

vec3 convert_cubemap_index(vec3 idx, int size)
{
  vec3 scale;
  if (abs(idx.x) >= abs(idx.y)) {
    if (abs(idx.x) >= abs(idx.z))
      scale = vec3(size, size - 1, size - 1);
    else
      scale = vec3(size - 1, size - 1, size);
  } else {
    if (abs(idx.y) >= abs(idx.z))
      scale = vec3(size - 1, size, size - 1);
    else
      scale = vec3(size - 1, size - 1, size);
  };
  return idx * scale / size;
}

The following picture illustrates the two related conversions.