A video recording of the final shader.

Protoplanet

This story begins four and a half billion years ago, with a lump of molten rock…

The early earth was a protoplanet, red hot and heavily cratered by asteroid impacts. As my earth simulation is entirely procedurally generated, with no pre-rendered textures, the first task is to generate a map of this terrain. To calculate the height of the terrain at a given latitude and longitude, first translate to 3D cartesian coordinates:

vec3 p = 1.5 * vec3(
    sin(lon*PI/180.) * cos(lat*PI/180.),
    sin(lat*PI/180.),
    cos(lon*PI/180.) * cos(lat*PI/180.));

Now, as asteroids come in a variety of sizes, so do the resulting craters. To accommodate this, the shader iterates over five levels of detail, layering craters of decreasing size over each other. To make the craters have a realistic rugged appearance, this is mixed with some fractional Brownian motion noise, and scaled so that the largest craters have the most impact on the terrain.

float height = 0.;
for (float i = 0.; i < 5.; i++) {
    float c = craters(0.4 * pow(2.2, i) * p);
    float noise = 0.4 * exp(-3. * c) * FBM(10. * p);
    float w = clamp(3. * pow(0.4, i), 0., 1.);
    height += w * (c + noise);
}
height = pow(height, 3.);

The craters themselves are generated on a 3D grid, from which a sphere is carved out for the surface terrain. To avoid visible regularity, the crater centres are given a pseudo-random offset from the grid points, using a hash function. To calculate influence of a crater at a given location, take a weighted average of the craters belonging to the nearby grid points, with weights exponentially decreasing with distance from the centre. The crater rims are generated by a simple sine curve.

float craters(vec3 x) {
    vec3 p = floor(x);
    vec3 f = fract(x);
    float va = 0.;
    float wt = 0.;
    for (int i = -2; i <= 2; i++)
     for (int j = -2; j <= 2; j++)
      for (int k = -2; k <= 2; k++) {
        vec3 g = vec3(i,j,k);
        vec3 o = 0.8 * hash33(p + g);
        float d = distance(f - g, o);
        float w = exp(-4. * d);
        va += w * sin(2.*PI * sqrt(d));
        wt += w;
    }
    return abs(va / wt);
}

The final procedurally generated heightmap looks like this:

Although relatively simple, after filling the low-lying regions with water, this procedural terrain resembles what scientists believe the early earth actually looked like:

Artistic impression of the early earth, by NASA.

Water contained within was vaporised by the heat, which escaped and began circulating through the early atmosphere forming around the planet. As time progressed and the rock cooled, the water vapour began to condense into oceans. The flow of liquid water across the surface carved valleys in the terrain, leaving an accumulation of sediment in its wake.

Tectonic plates

The formation of mountains, ocean trenches, and familiar continental landforms requires a model of tectonic movement. The simulation randomly generates seed locations for plates, with an initial velocity. These plates grow in size over time with a simple aggregation model, which randomly selects neighbouring points and adds them to a plate if they have not already been assigned to another plate. All of the pixels within a plate store the velocity of the plate’s movement. The aggregation model is similar to that of a diffusion-limited aggregation (but without the diffusion):

Continuous movement of the plates is difficult, as it would require plate boundaries to account for movements measured in fractions of a pixel. To avoid this, the plates are instead moved at discrete time-steps, by a whole pixel either horizontally or vertically. These times are randomised for each plate such that the average velocity is maintained at the set speed and direction, and also so that it is unlikely that neighbouring plates will move simultaneously.

Plate collisions occur when some boundary pixels of one plate move onto a location previously occupied by pixels belonging to another plate. This causes subduction, which is modelled by simply slightly increasing the elevation of the terrain at the locations of the collision. Although this only occurs at the pixels along the boundary of a plate, the impact is gradually spread to neighbouring pixels through a simple thermal erosion model, which pushes the elevation of a pixel in the direction of the average of its neighbours.

Altogether this provides a decent simulation of the formation of continents with mountain ranges (which will be further improved with the introduction of hydraulic erosion in the next section):