[GIS] Best flow routing algorithm for determining location where spill enters stream

algorithmroutingtaudem

What would be the best flow routing algorithm for determining the quickest (and by this I suppose I mean shortest) path that a spill would take when traveling to a stream?

My best educated guess would be that the D8 would provide the quickest path for a spill traveling from land to the stream. It minimizes dispersion (which would most likely slow the spill down) and maximizes velocity (ignoring surface friction due to different landcover types – if this can be ignored) by choosing the steepest path.

The end application would be in determining the overland flow path from a spill site to a stream/river reach that provides the earliest possible time for entry of the spill into the stream/river. There are different assumptions or starting datasets that could change the answer. For instance, a 30m digital elevation model (DEM) cold prove useful in a forested or agricultural area, but a 1m or sub-1m DEM might be necessary for sorting out the true flow path (and therefore the true shortest path) in a highly urbanized area. In either of these cases, however, once one has a DEM that he/she feels represents the terrain well enough for an analysis, which algorithm should they choose?

This query is borne out of trying to avoid calculating the distance from the spill to the stream/river reach using simple straight line geometry in a GIS tool, such as the proximity toolset in ArcMap.

Does anyone have any experience with this?

EDIT: Based on the responses so far, it seems that using a routing algorithm that allows for dispersion, and then calculating the shortest path out of the possible paths predicted by the algorithm, might be the solution I want. In the answer suggesting the use of TauDEM, the journal article that is linked to contains this image:

enter image description here

(image from Tesfa et al., 2011 – http://www.neng.usu.edu/cee/faculty/dtarb/Tesfa_EMS2011.pdf)

Examining the predicted flow path of the D-Inf algorithm, we see that several different paths are possible. Using such an algorithm, I suppose on could calculate the shortest path, the longest path, and a range of other statistics on the predicted pathways. If similar or near identical ground cover is found across the different pathways, then some combination of steepness and path length may dictate the quickest path? Which might be different from the shortest path? I'm starting to think I have to review my question and/or thought process – it may be necessary to incorporate equations/models governing the speed of overland flow along side the algorithms that examine possible flow pathways. I need to keep reading and ingesting the article, but I think that the answer lies within the ideas presented in the article – so for now I will mark the suggesting answer as correct.

Best Answer

Some good software for this purpose is TauDEM (Terrain Analysis Using Digital Elevation Models).

Main page, with info, publications, and downloads for command tools and optional ArcGIS Toolbox interface
Source code on GitHub

The software uses both D8 and so-called D-infinity (D∞) flow models, which is best illustrated on Fig. 1 of Tesfa et al. (2011). You can test which flow path method is quicker or more accurate for your area. Furthermore, both serial and parallel processing with any number of CPUs is supported. Section 2.2 of Tesfa et al. (2011) linked previously covers much of your questions of flow routing distances down to a stream (end point cells).

Related Solutions

[GIS] ny pseudocode for MFD (multiple flow direction) algorithm available

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.

[GIS] Should I implement TSP or Dijkstra

If I understand you want to implement a nearest-neighbor algorithm in java(?) Nearest-neighbor is a well known not very efficient tsp (or vrp) solving approach... And you have to (is it really a constraint?) use OSRM to calculate the distance between all your points. OSRM is a c++ program which uses boost and is osm compliant. Maybe you should consider using something else, like pgrouting for example, using the same (openstreetmap) data and contains various shortest-path options; it would save you the java native interface implementation (and be more easily cross-platform able).

Dijkstra creates a tree starting from a given point. This tree can be used to calculate the distances between your starting point and all the other vertices you want (if a route exists in your graph) by adding the weigth of each edge. It seems a fair way to generate a distance matrix.

I think you should have a look at some VRP algorithms... Good luck.

Best Answer

Related Solutions

[GIS] ny pseudocode for MFD (multiple flow direction) algorithm available

[GIS] Should I implement TSP or Dijkstra

Related Question