Is there any pseudocode for MFD (multiple flow direction) algorithm available?

Question

Still stuck with the simulation system, is there any good pseudocode for the multiple flow direction (MFD) algorithm? I read a pseudocode example and implemented it, however the simulation was linear: that is the first cell gets the highest amount of rainfall and the last cell does not get any water. I'm resorting to sorting the elevation of the cells from highest to lowest so that the water would spread out properly. I don't want to use r.terraflow or TOPMODEL since this is an independent project and I don't want to use existing simulation systems.

Any ideas as to how to implement MFD or correct my algorithm would be very appreciated.

BTW, null values were set to -500 by the DEM reader. The code is implemented in C# and is not yet optimized, I just want a correct MFD simulation at the moment.

---

static void rainfallAlgo() {
  double total, percentage, transfer;
  double[] tile = new double[10];

  for (int i = 0; i != Convert.ToInt32(totalRow); i++){
    for (int j = 0; j != Convert.ToInt32(totalCol); j++){
      if ((elev[i, j] != -500)&&(elevrain[i,j]>=elev[i,j])){
       Console.WriteLine();
        total = 0;
        for (int k = 1; k <= 9; k++){
          if (k != 5){
            tile[k] = compare(i, j, k);
            total += tile[k];
          }
        }
        for (int k = 1; k <= 9; k++){
          if (total == 0){
            Console.WriteLine("I am ZERO");
            total = 1;
          }
          if (k != 5){
            percentage = tile[k] / total;
            Console.Write(k+ ":");
            Console.WriteLine(percentage);
            transfer = percentage * (elevrain[i,j]-elev[i,j]);
            if (k == 1)
              elevrain[i - 1, j - 1] += transfer;
            else if (k == 2)
              elevrain[i - 1, j] += transfer;
            else if (k == 3)
              elevrain[i - 1, j + 1] += transfer;
            else if (k == 4)
              elevrain[i, j - 1] += transfer;
            else if (k == 6)
              elevrain[i, j + 1] += transfer;
            else if (k == 7)
              elevrain[i + 1, j - 1] += transfer;
            else if (k == 8)
              elevrain[i + 1, j] += transfer;
            else if (k == 9)
              elevrain[i + 1, j + 1] += transfer;
          }
        }
        if(total!=1)
          elevrain[i, j] = elev[i, j];
        Console.WriteLine(elevrain[i,j]);
      }
    }
  }
}


static void addRain(){
  rainfall = new double[Convert.ToInt32(totalRow), Convert.ToInt32(totalCol)];
  elevrain = new double[Convert.ToInt32(totalRow), Convert.ToInt32(totalCol)];

  int i = 0, j = 0;

  while (i != Convert.ToInt32(totalRow)){
    while (j != Convert.ToInt32(totalCol)){
      rainfall[i, j] = 2;
      elevrain[i, j] = elev[i, j] + rainfall[i,j];
      j++;
    }
    i++;
    j = 0;
  }
}

static double compare(int i, int j, int k){
  double total=0;
  if( k==1 ) {
    if ( (elevrain[i, j] > elevrain[i - 1, j - 1] ) && (elev[i-1, j-1]!=-500.0) )
      total = (elevrain[i, j] - elevrain[i - 1, j - 1]) / 1.414213562373095;
    else
      total = 0;
  } else if( k==2 ) {
    if ( (elevrain[i, j] > elevrain[i - 1, j]) && (elev[i - 1, j] != -500.0))
      total = elevrain[i, j] - elevrain[i - 1, j];
    else
      total = 0;
  } else if( k==3 ) {
    if ( (elevrain[i, j] > elevrain[i - 1, j + 1]) && (elev[i - 1, j + 1] != -500.0) )
      total = (elevrain[i, j] - elevrain[i - 1, j + 1]) / 1.414213562373095;
    else
      total = 0;
  } else if( k==4 ) {
    if ((elevrain[i, j] > elevrain[i, j - 1]) && (elev[i, j - 1] != -500.0) )
      total = elevrain[i, j] - elevrain[i, j - 1];
    else
      total = 0;
  } else if( k== 6 ) {
    if ( (elevrain[i, j] > elevrain[i, j + 1]) && (elev[i, j + 1] != -500.0) ) 
      total = elevrain[i, j] - elevrain[i, j + 1];
    else
      total = 0;
  } else if( k== 7 ) {
    if ( (elevrain[i, j] > elevrain[i + 1, j - 1]) && (elev[i + 1, j-1] != -500.0) )
      total = (elevrain[i, j] - elevrain[i + 1, j - 1]) / 1.414213562373095;
    else
      total = 0;
  } else if ( k== 8 ) {
    if ((elevrain[i, j] > elevrain[i + 1, j]) && (elev[i + 1, j] != -500.0) )
      total = elevrain[i, j] - elevrain[i + 1, j];
    else
      total = 0;
  } else if ( k== 9) {
    if ( (elevrain[i, j] > elevrain[i + 1, j + 1]) && (elev[i + 1, j +1] != -500.0) )
      total = (elevrain[i, j] - elevrain[i + 1, j + 1]) / 1.414213562373095;
    else
      total = 0;
  }
  return total;
}

Answer 1

Your question is poorly stated. It as at once general, as it asks if anyone has ideas about how to implement MFDs, and specific, as it seems to ask what is wrong with you code. I'll try to answer this second conceptualization of your question here.

You claim that the code you've posted is not optimized. I have rewritten the code below as I think you should have written it. The changes I have introduced include making good use of loops and offset arrays. These changes are not optimizations, they are helpful abstractions that allow you to work with the code and review it in meaningful ways. Without such abstractions you are not coding, you are producing functional gibberish.

const int di[10]={0,-1,-1,-1, 0,0,0, 1,1,1};
const int dj[10]={0,-1, 0, 1,-1,0,1,-1,0,1};


static void rainfallAlgo() {
  double total, percentage, transfer;
  double[] tile = new double[10];

  //Consider declaring -500 as a const so that it reads intelligbly
  //For instance, you could write "#define NoData -500"
  //Then say, "if elev[i,j]!=NoData"

  for (int i = 0; i != Convert.ToInt32(totalRow); i++){
    for (int j = 0; j != Convert.ToInt32(totalCol); j++){
      if ((elev[i, j] != -500)&&(elevrain[i,j]>=elev[i,j])){
       Console.WriteLine();
        total = 0;
        for (int k = 1; k <= 9; k++){
          if (k != 5){
            tile[k] = compare(i, j, k);
            total += tile[k];
          }
        }
        for (int k = 1; k <= 9; k++){
          if (total == 0){
            Console.WriteLine("I am ZERO");
            total = 1;
          }
          if (k != 5){
            percentage = tile[k] / total;
            Console.Write(k+ ":");
            Console.WriteLine(percentage);
            transfer = percentage * (elevrain[i,j]-elev[i,j]);
            elevrain[i+di[k],j+dj[k]]+=transfer;
          }
        }
        if(total!=1)
          elevrain[i, j] = elev[i, j];
        Console.WriteLine(elevrain[i,j]);
      }
    }
  }
}

//This function isn't called from anywhere; I don't know what it does
static void addRain(){
  rainfall = new double[Convert.ToInt32(totalRow), Convert.ToInt32(totalCol)];
  elevrain = new double[Convert.ToInt32(totalRow), Convert.ToInt32(totalCol)];

  for(int i=0;i != Convert.ToInt32(totalRow);++i)
  for(int j=0;j != Convert.ToInt32(totalCol);++j){
    rainfall[i, j] = 2;
    elevrain[i, j] = elev[i, j] + rainfall[i,j];
}

static double compare(int i, int j, int k){
  const int dr[9]={0,1.414,0,1.414,0,0,1.414,0,1.414,0,1.414,0};
  double total=0;
  if ( (elevrain[i, j] > elevrain[ i+di[k], j+dj[k] ] ) && (elev[ i+di[k], j+dj[k] ]!=-500.0) )
    total = (elevrain[i, j] - elevrain[ i+di[k], j+dj[k] ]) / dr[k];
  return total;
}

With the introduction of the arrays di and dj it becomes more reasonable to analyze your code.

I am left wondering why you have chosen to label your central cell as 5 as this makes your loops much more painful. It is better to use 0 as the central label so that you can loop sequentially over all of the neighbouring cells.

I have also used an array dr to store the offset distances between the centers of the central cell and its neighbours.

I think you will find that the judicious use of arrays such as those I've introduced here will lead to better code and that better code will lead you more quickly to the solutions you seek.

You can also wrap your grids into classes which provide meaningful accessor names to help you follow what's going on. You can take a look at some of my code here and here to get ideas.

Happy coding.