I'm working on a model of landscape evolution and would like to make a nice 3D movie of how the landscape over time. How do I do this?

My landscape evolution model, which is based on this paper and its accompanying code, is as follows:

#include <array>
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <fenv.h> //Used to catch floating point NaN issues
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <random>
#include <omp.h>  //Used for OpenMP run-time functions
#include <vector>

//Used to handle situations in which OpenMP is not available
//(This scenario has not been extensively tested)
#ifndef _OPENMP
  #define omp_get_thread_num()  0
  #define omp_get_num_threads() 1
  #define omp_get_max_threads() 1

///This is a quick-and-dirty, zero-dependency function for saving the outputs of
///the model in matrix data format.
void PrintXYZ_CSV(
  const std::string filename,
  const std::vector<double>& h,
  const int width,
  const int height
  std::ofstream fout(filename.c_str());
  //Since the outer ring of the dataset is a halo used for simplifying
  //neighbour-finding logic, we do not save it to the output here.
  for(int y=1;y<height-1;y++){
    for(int x=1;x<width-1;x++)

///The entire model is contained in a handy class, which makes it easy to set up
///and solve many such models.
class FastScape_RBPQ {
  //Value used to indicate that a cell had no downhill neighbour and, thus, does
  //not flow anywhere.
  const int    NO_FLOW = -1;
  const double SQRT2   = 1.414213562373095048801688724209698078569671875376948; //Yup, this is overkill.

  //NOTE: Having these constants specified in the class rather than globally
  //results in a significant speed loss. However, it is better to have them here
  //under the assumption that they'd be dynamic in a real implementation.
  const double keq       = 2e-6;   //Stream power equation constant (coefficient)
  const double neq       = 2;      //Stream power equation constant (slope modifier)
  const double meq       = 0.8;    //Stream power equation constant (area modifier)
  const double ueq       = 2e-3;   //Rate of uplift
  const double dt        = 1000.;  //Timestep interval
  const double dr[8]     = {1,SQRT2,1,SQRT2,1,SQRT2,1,SQRT2}; //Distance between adjacent cell centers on a rectangular grid arbitrarily scale to cell edge lengths of 1
  const double tol       = 1e-3;   //Tolerance for Newton-Rhapson convergence while solving implicit Euler
  const double cell_area = 40000;  //Area of a single cell

  int width;                 //Width of DEM
  int height;                //Height of DEM
  int size;                  //Size of DEM (width*height)
  unsigned long rand_seed;   //Random seed used to generate landscape
  std::string output_prefix; //Start of output file: PREFIX-#.csv

  //The dataset is set up so that the outermost edge is never actually used:
  //it's a halo which allows every cell that is actually processed to consider
  //its neighbours in all 8 directions without having to check first to see if
  //it is an edge cell. The ring of second-most outer cells is set to a fixed
  //value to which everything erodes (in this model).

  //Rec directions (also used for nshift offsets) - see below for details
  //1 2 3
  //0   4
  //7 6 5

  std::vector<double> h;        //Digital elevation model (height)
  std::vector<double> accum;    //Flow accumulation at each point
  std::vector<int>    rec;      //Direction of receiving cell
  std::vector<int>    donor;    //Indices of a cell's donor cells
  std::vector<int>    ndon;     //How many donors a cell has
  std::array<int,8>   nshift;   //Offset from a focal cell's index to its neighbours in terms of flat indexing

  int stack_width;  //Number of cells allowed in the stack
  int level_width;  //Number of cells allowed in a level

  void GenerateRandomTerrain(){
    std::mt19937 eng;

    std::uniform_real_distribution<> d;
    using parm_t = std::uniform_real_distribution<>::param_type;

    for(int y=0;y<height;y++)
    for(int x=0;x<width;x++){
      const int c = y*width+x;
      h[c] = d(eng, parm_t{0, 1});

      //Outer edge is set to 0 and never touched again. It is used only as a
      //convenience so we don't have to worry when a focal cell looks at its
      if(x == 0 || y==0 || x==width-1 || y==height-1)
        h[c] = 0;

      //Second outer-most edge is set to 0 and never touched again. This is the
      //baseline to which all cells would erode where it not for uplift. You can
      //think of this as being "sea level".
      if(x == 1 || y==1 || x==width-2 || y==height-2)
        h[c] = 0;

  ///Initializing code
    const int width0,
    const int height0,
    const std::string &output_prefix,
    const unsigned long rand_seed
    //Initialize code for finding neighbours of a cell
    : nshift{-1,-width0-1,-width0,-width0+1,1,width0+1,width0,width0-1},
    width  = width0;
    height = height0;
    size   = width*height;

    h.resize(size);   //Memory for terrain height

    GenerateRandomTerrain();     //Could replace this with custom initializer

  ///The receiver of a focal cell is the cell which receives the focal cells'
  ///flow. Here, we model the receiving cell as being the one connected to the
  ///focal cell by the steppest gradient. If there is no local gradient, than
  ///the special value NO_FLOW is assigned.
  void ComputeReceivers(){
    //Edge cells do not have receivers because they do not distribute their flow
    //to anywhere.

    #pragma omp for collapse(2)
    for(int y=2;y<height-2;y++)
    for(int x=2;x<width-2;x++){
      const int c      = y*width+x;

      //The slope must be greater than zero for there to be downhill flow;
      //otherwise, the cell is marked NO_FLOW.
      double max_slope = 0;        //Maximum slope seen so far amongst neighbours
      int    max_n     = NO_FLOW;  //Direction of neighbour which had maximum slope to focal cell

      //Loop over neighbours
      for(int n=0;n<8;n++){
        const double slope = (h[c] - h[c+nshift[n]])/dr[n]; //Slope to neighbour n
        if(slope>max_slope){    //Is this the steepest slope we've seen?
          max_slope = slope;    //If so, make a note of the slope
          max_n     = n;        //And which cell it came from
      rec[c] = max_n;           //Having considered all neighbours, this is the steepest

  ///The donors of a focal cell are the neighbours from which it receives flow.
  ///Here, we identify those neighbours by inverting the Receivers array.
  void ComputeDonors(){
    //The B&W method of developing the donor array has each focal cell F inform
    //its receiving cell R that F is a donor of R. Unfortunately, parallelizing
    //this is difficult because more than one cell might be informing R at any
    //given time. Atomics are a solution, but they impose a performance cost
    //(though using the latest and greatest hardware decreases this penalty).

    //Instead, we invert the operation. Each focal cell now examines its
    //neighbours to see if it receives from them. Each focal cell is then
    //guaranteed to have sole write-access to its location in the donor array.

    //Remember, the outermost ring of cells is a convenience halo, so we don't
    //calculate donors for it.

    #pragma omp for collapse(2) schedule(static) nowait
    for(int y=1;y<height-1;y++)
    for(int x=1;x<width-1;x++){
      const int c = y*width+x;
      ndon[c] = 0; //Cell has no donor neighbours we know about
      for(int ni=0;ni<8;ni++){
        const int n = c+nshift[ni];
        //If the neighbour has a receiving cell and that receiving cell is
        //the current focal cell c
        if(rec[n]!=NO_FLOW && n+nshift[rec[n]]==c){
          donor[8*c+ndon[c]] = n;

  ///Cells must be ordered so that they can be traversed such that higher cells
  ///are processed before their lower neighbouring cells. This method creates
  ///such an order. It also produces a list of "levels": cells which are,
  ///topologically, neither higher nor lower than each other. Cells in the same
  ///level can all be processed simultaneously without having to worry about
  ///race conditions.
  void GenerateOrder(
    std::vector<int>& stack,
    std::vector<int>& levels,
    int &nlevel
    int nstack = 0;

    levels[0]  = 0;
    nlevel     = 1;

    //Outer edge
    #pragma omp for schedule(static) nowait
    for(int y=1;y<height-1;y++){
      stack[nstack++] = y*width+1;          assert(nstack<stack_width);
      stack[nstack++] = y*width+(width-2);  assert(nstack<stack_width);

    #pragma omp for schedule(static) nowait
    for(int x=2;x<width-2;x++){
      stack[nstack++] =          1*width+x; assert(nstack<stack_width);
      stack[nstack++] = (height-2)*width+x; assert(nstack<stack_width);

    //End of outer edge
    levels[nlevel++] = nstack; //Last cell of this level

    //Interior cells
    //TODO: Outside edge is always NO_FLOW. Maybe this can get loaded once?
    //Load cells without dependencies into the queue
    //TODO: Why can't I use nowait here?
    #pragma omp for collapse(2) schedule(static)
    for(int y=2;y<height-2;y++)
    for(int x=2;x<width -2;x++){
      const int c = y*width+x;
        stack[nstack++] = c;
    levels[nlevel++] = nstack; //Last cell of this level

    //Start with level_bottom=-1 so we get into the loop, it is immediately
    //replaced by level_top.
    int level_bottom = -1;         //First cell of the current level
    int level_top    =  0;         //Last cell of the current level

    #pragma omp barrier //Must ensure ComputeDonors is done before we start the following

    while(level_bottom<level_top){ //Ensure we parse all the cells
      level_bottom = level_top;    //The top of the previous level we considered is the bottom of the current level
      level_top    = nstack;       //The new top is the end of the stack (last cell added from the previous level)
      for(int si=level_bottom;si<level_top;si++){
        const auto c = stack[si];
        //Load donating neighbours of focal cell into the stack
        for(int k=0;k<ndon[c];k++){
          const auto n = donor[8*c+k];
          stack[nstack++] = n;

      levels[nlevel++] = nstack; //Start a new level

    //End condition for the loop places two identical entries
    //at the end of the stack. Remove one.


  ///Compute the flow accumulation for each cell: the number of cells whose flow
  ///ultimately passes through the focal cell multiplied by the area of each
  ///cell. Each cell could also have its own weighting based on, say, average
  void ComputeFlowAcc(
    const std::vector<int>& stack,
    const std::vector<int>& levels,
    const int &nlevel
    for(int i=levels[0];i<levels[nlevel-1];i++){
      const int c = stack[i];
      accum[c] = cell_area;

    //Highly-elevated cells pass their flow to less elevated neighbour cells.
    //The queue is ordered so that higher cells are keyed to higher indices in
    //the queue; therefore, parsing the queue in reverse ensures that fluid
    //flows downhill.

    //We can process the cells in each level in parallel. To prevent race
    //conditions, each focal cell figures out what contirbutions it receives
    //from its neighbours.

    //`nlevel-1` is the upper bound of the stack.
    //`nlevel-2` through `nlevel-1` are the cells which have no higher neighbours (top of the watershed)
    //`nlevel-3` through `nlevel-2` are the first set of cells with higher neighbours, so this is where we start
    for(int li=nlevel-3;li>=1;li--){
      const int lvlstart = levels[li];      //Starting index of level in stack
      const int lvlend   = levels[li+1];    //Ending index of level in stack
      for(int si=lvlstart;si<lvlend;si++){
        const int c = stack[si];
        for(int k=0;k<ndon[c];k++){
          const auto n = donor[8*c+k];
          accum[c]    += accum[n];

  ///Raise each cell in the landscape by some amount, otherwise it wil get worn
  ///flat (in this model, with these settings)
  void AddUplift(
    const std::vector<int>& stack,
    const std::vector<int>& levels,
    const int &nlevel
    //We exclude two exterior rings of cells in this example. The outermost ring
    //(the edges of the dataset) allows us to ignore the edges of the dataset,
    //the second-most outer ring (the cells bordering the edge cells of the
    //dataset) are fixed to a specified height in this model. All other cells
    //have heights which actively change and they are altered here.

    //Start at levels[1] so we don't elevate the outer edge
    #pragma omp simd
    for(int i=levels[1];i<levels[nlevel-1];i++){
      const int c = stack[i];
      h[c]       += ueq*dt;

  ///Decrease he height of cells according to the stream power equation; that
  ///is, based on a constant K, flow accumulation A, the local slope between
  ///the cell and its receiving neighbour, and some judiciously-chosen constants
  ///m and n.
  ///    h_next = h_current - K*dt*(A^m)*(Slope)^n
  ///We solve this equation implicitly to preserve accuracy
  void Erode(
    const std::vector<int>& stack,
    const std::vector<int>& levels,
    const int nlevel
    //The cells in each level can be processed in parallel, so we loop over
    //levels starting from the lower-most (the one closest to the NO_FLOW cells)

    //#pragma omp parallel default(none)
    for(int li=2;li<nlevel-1;li++){
      const int lvlstart = levels[li];
      const int lvlend   = levels[li+1];
      #pragma omp simd
      for(int si=lvlstart;si<lvlend;si++){
        const int c = stack[si];         //Cell from which flow originates
        const int n = c+nshift[rec[c]];  //Cell receiving the flow

        const double length = dr[rec[c]];
        //`fact` contains a set of values which are constant throughout the integration
        const double fact   = keq*dt*std::pow(accum[c],meq)/std::pow(length,neq);
        const double h0     = h[c];      //Elevation of focal cell
        const double hn     = h[n];      //Elevation of neighbouring (receiving, lower) cell
        double hnew         = h0;        //Current updated value of focal cell
        double hp           = h0;        //Previous updated value of focal cell
        double diff         = 2*tol;     //Difference between current and previous updated values
        while(std::abs(diff)>tol){       //Newton-Rhapson method (run until subsequent values differ by less than a tolerance, which can be set to any desired precision)
          hnew -= (hnew-h0+fact*std::pow(hnew-hn,neq))/(1.+fact*neq*std::pow(hnew-hn,neq-1));
          diff  = hnew - hp;             //Difference between previous and current value of the iteration
          hp    = hnew;                  //Update previous value to new value
        h[c] = hnew;                     //Update value in array


  ///Run the model forward for a specified number of timesteps. No new
  ///initialization is done. This allows the model to be stopped, the terrain
  ///altered, and the model continued. For space-efficiency, a number of
  ///temporary arrays are created each time this is run, so repeatedly running
  ///this function for the same model will likely not be performant due to
  ///reallocations. If that is your use case, you'll want to modify your code
  void run(const int nstep){
    stack_width = std::max(300000,5*size/omp_get_max_threads()); //Number of stack entries available to each thread
    level_width = std::max(1000,size/omp_get_max_threads());     //Number of level entries available to each thread

    accum.resize(  size);  //Stores flow accumulation
    rec.resize  (  size);  //Array of Receiver directions
    ndon.resize (  size);  //Number of donors each cell has
    donor.resize(8*size);  //Array listing the donors of each cell (up to 8 for a rectangular grid)

    //! initializing rec
    #pragma omp parallel for
    for(int i=0;i<size;i++)
      rec[i] = NO_FLOW;

    #pragma omp parallel for
    for(int i=0;i<size;i++)
      ndon[i] = 0;

    #pragma omp parallel
      std::vector<int> stack(stack_width); //Indices of cells in the order they should be processed

      //A level is a set of cells which can all be processed simultaneously.
      //Topologically, cells within a level are neither descendents or ancestors
      //of each other in a topological sorting, but are the same number of steps
      //from the edge of the dataset.

      //It's difficult to know how much memory should be allocated for levels. For
      //a square DEM with isotropic dispersion this is approximately sqrt(E/2). A
      //diagonally tilted surface with isotropic dispersion may have sqrt(E)
      //levels. A tortorously sinuous river may have up to E*E levels. We
      //compromise and choose a number of levels equal to the perimiter because
      //why not?
      std::vector<int> levels(level_width);
      int  nlevel = 0;

      for(int step=0;step<=nstep;step++){
        ComputeReceivers  ();
        ComputeDonors     ();
        GenerateOrder     (stack,levels,nlevel);
        ComputeFlowAcc    (stack,levels,nlevel);
        AddUplift         (stack,levels,nlevel);
        Erode             (stack,levels,nlevel);
        #pragma omp barrier //Ensure threads synchronize after erosion so we calculate receivers correctly

        PrintXYZ_CSV(output_prefix + "-" + std::to_string(step) + ".csv", h, width, height);

        #pragma omp master
        if( step%20==0 ) //Show progress
          std::cout<<"p Step = "<<step<<std::endl;

    //Free up memory, except for the resulting landscape height field prior to
    //exiting so that unnecessary space is not used when the model is not being
    accum .clear();   accum .shrink_to_fit();
    rec   .clear();   rec   .shrink_to_fit();
    ndon  .clear();   ndon  .shrink_to_fit();
    donor .clear();   donor .shrink_to_fit();

  ///Returns a pointer to the data so that it can be copied, printed, &c.
  std::vector<double>& getH() {
    return h;

int main(int argc, char **argv){
  //Enable this to stop the program if a floating-point exception happens

    std::cerr<<"Syntax: "<<argv[0]<<" <Dimension> <Steps> <Output Prefix> <Seed>"<<std::endl;
    return -1;

  const int         width         = std::stoi (argv[1]);
  const int         height        = std::stoi (argv[1]);
  const int         nstep         = std::stoi (argv[2]);
  const std::string output_prefix = argv[3] ;
  const auto        rand_seed     = std::stoul(argv[4]);

  //Uses the RichDEM machine-readable line prefixes
  //Name of algorithm
  std::cout<<"A FastScape RB+PQ"<<std::endl;
  //Citation for algorithm
  std::cout<<"C Richard Barnes TODO"<<std::endl;
  //Random seed used to produce outputs
  std::cout<<"m Random seed = "<<rand_seed<<std::endl;

  FastScape_RBPQ tm(width, height, output_prefix, rand_seed);

  return 0;

It can be compiled with:

g++ -O3 fastscape_RB+PQ.cpp -fopenmp

and run with:

./a.out 100 120 landscape 1234456

to generate a test dataset.

What is this model doing?

For those wondering, what's happening here is that a landscape is initialized to have some random noise.

Then, at each timestep, the landscape experiences uplift, which prevents it from being eroded away entirely over the course of this simulation. After the uplift takes place flow accumulation---the amount of upstream water that eventually passes through each cell in the grid-discretized landscape---is calculated.

Next, the landscape is evolved according to the stream power equation:

Stream power equation

where h is the elevation of a cell of the grid-discretized landscape, K is a scalar whose value may be influenced by, for example, lithology (rock type), channel width, and channel hydrology; A is the flow accumulation; and m and n are scaling constants. Put another way, this equation says that erosion is proportional to the type of rock under the stream bed, the amount of water in the stream, and the force with which that water moves (as determined by the local slope of the stream).

The equation is solved using a backward's Euler method. Assuming that the boundary of the landscape is at a fixed elevation provides the boundary condition needed to solve the equation at the edge of the domain. Once the edge is solved, the next set of cells in can be solved, and so on.

1 Answer 1


Probably the best tool for this is Paraview. Once you've downloaded and acquired that you'll want to load your images:

Loading images in Paraview

Use the CSV reader to load them:

Paraview CSV filter

Below, you can find filters faster by pressing Ctrl+Space and starting to type the name of the filter you want.

Next, apply the TableToPoints filter:

Using the TableToPoints filter

Elevation often needs to be exaggerated to get a good looking image, so let's apply a transform filter to do that:

Applying the transform filter

You may also want to clip off elevations which are too high. To do that, use the Clip filter with a "Z Normal" plane and an origin at the height you want to clip from:

Applying the clip filter

To get colours, you'll want to add the Elevation filter so that Paraview interprets the z-axis as elevation (be sure to click the Z axis button):

Using the elevation filter

Using the Delaunay2D filter converts the points into a surface coloured by elevation:

Applying the Delaunay2D filter

The surface is kinda jagged, so let's smooth it (note that I've cranked the "Number of" to 100 here and stepped forward in time to show its effect more clearly):

Using the smoothing filter

After all, this close the Spreadsheet window. Hopefully you're now in a "RenderView". The axes might show that you're not seeing all your data, we need to expand the view so that you do. There's a button labeled "Zoom to Data" that's helpful here:

Zoom to Data in Paraview

After clicking the button, you may need to click the eye next to the "Smooth" filter to re-display your data.

Unchecking the "Axis Grid" and "Orientation Axes Visibility" cleans up the image somewhat while going back to the "Smooth" filter provides an option to disable the colour legend (it's in the bottom right of the following image):

Disable colour legend in Paraview

The "Play" button at the top will show you the animation. To save it, click File->Save Animation and choose appropriate settings. Resizing the window resizes the view in the output animation (whatever the window shows is what the animation gets).

This leaves you with a large number of files of the form, e.g., "landscape.0027.png".

We make them into a video using ffmpeg:

ffmpeg -framerate 6 -pattern_type glob -i 'landscape*.png' -c:v libx264 -pix_fmt yuv420p -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" out.mp4


  • -framerate 6 indicates that 6 images should be shown per second
  • -pattern_type glob -i 'landscape*.png' identifies the files to be used (their lexicographic order should also be their numeric order thanks to the padding zeros)
  • -c:v libx264 -pix_fmt yuv420p some stuff to do with the encoding format
  • -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" shim to deal with images whose dimensions are not evenly divisible by 2
  • out.mp4 output filename

We could also use ImageMagick to create an animated GIF:

convert -delay 10 landscape*png -loop 0 movie.gif


  • -delay 10 Delay 10/100ths of a second between frames
  • landscape*png Images to animate
  • -loop 0 Loop the GIF an infinite number of times
  • movie.gif Output filename

And all of this gives something lovely

Animation of landscape evolution

  • 1
    This is really a phenomenal answer and question. Do you have any background on what this DEM time-series is showing? Erosion of some sort?
    – Aaron
    Commented Jan 24, 2021 at 18:03
  • Thanks, @Aaron: I've edited the question to provide an explanation of the model.
    – Richard
    Commented Jan 24, 2021 at 18:45
  • Wow. Nice answer! Commented Jan 27, 2021 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.