1

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;
}
Improve this question
6
  • 2
    Which MFD? E.g., John Wilson (2002) lists six different ones and then proposes four more! (crwr.utexas.edu/terrainAnalysis/presentations/121802/… , slides 9&10). And what precisely are you asking for: you are asking for "pseudocode" in the title yet in the second line of the question you say you (already) have "pseudocode"! – whuber Mar 6 '11 at 16:17
  • I am looking for a better pseudocode since the one that I have presented does it incorrectly. So there might be someone else who have encountered this problem and have solved it. – teamAngat Mar 7 '11 at 9:27
  • 1
    I don't see any pseudocode in your question. There's a lot of actual code, but it's not terribly helpful (and it's far too long and uncommented to be part of any readable question). One would expect pseudocode to be brief and descriptive. – whuber Mar 7 '11 at 20:00
  • As I said this is already an implementation of the pseudocode that I have read. :) – teamAngat Mar 8 '11 at 17:31
  • Both SAGA GIS and Grass GIS contain implementations of MFD algorithms. You could read their code to see how they implement the algorithm. – johanvdw Aug 30 '11 at 21:15
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.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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