Assignment: Computing FD and FA on a grid terrain

Develop C/C++ code to compute the FD and FA grids of an elevation grid. The names of the grid should be specified by the user on the commans line, something along these lines:

[ltoma@dover:\~] flow test2.asc test2FD.asc test2FA.asc

This will read elevation grid test2.asc, compute its FD and FA grids and save them as test2FD.asc and test2FA.asc. Unless a path is specified, files are by default in the current folder.

Flow direction and flow accumulation

Flow is one of the fundamental processes modeled on digital terrain models.The goal of modeling flow is to be able to answer questions such as: where does the water go (when it rains)? What are the areas susceptible to flooding? What is the part of that terrain that flows through a given point? What is the impact area when poluting a part of the terrain? What is the impact of storms, waves, and sea level rise?

Flow on grids is modeled using two basic parameters: flow direction and flow accumulation.

Flow direction: The flow direction at a grid cell (i,j), FD(i,j), represents the direction water flows at that cell. That is, if we drop a bucket of water on cell (i,j), it will all flow in the direction FD(i,j). The simplest model for FD is to assume that:

water always flows downhill
water flows towards the steepest downhill neighbor

Granted in reality flow direction may be more complicated (like water flowing in multiple direction, or on flat areas), but all in all this is a good model that gives good results, is easy to compute, and has nice properties.

Flow accumulation: Flow accumulation represents the accumulation of water, assuming that it flows using the computed FD. More precisely, assume that every cell in the grid is initialized with 1 unit of water, and that every cell in the grid sends the water it receives to the neighbor pointed to by its FD. That is, when a cell (i,j) receives 1 unit of water from a neighbor, it will send it along to the neighbor pointed to by FD(i,j). FA(i,j) represents how much total water flows through cell (i,j).

Given a terrain, the goal is to compute FD and FA for the entire terrain. For a grid terrain, flow direction and flow accumulation are grids of the same size as the elevation grid defined in the natural way:

Flow direction grid: For all (i,j), the value at (i,j) represents the FD of cell (i,j).
Flow accumulation grid: For all (i,j), the value at (i,j) represents the FA of cell (i,j).

Dealing with real data

Real grids will have NODATA values, flat areas, and cells of equal elevation. A NODATA value means that the value of the elevation at the point is unknown. On flat areas there are no downhill neighbors, so flow is technically undefined using the definition above. And, a point may have several steepest neighbors, not just one. We'll handle these cases with the following assumptions:

If a point is NODATA (undefined), its flow direction and flow accumulation should be NODATA (undefined); A NODATA point does not send its flow to any neighbor, so all NODATA points should end up with flow = NODATA.
If a point has no neighbors that are lower than it, then it does not send its flow anywhere. This point is (part of) a pit.
If when determining the steepest downslope neighbors there are ties, they are broken arbitrarily.

General outline

The goal is to design (structure) the code practicing as much as possible data abstraction and encapsulation. There should be a grid module that implements all functions on grids, such as reading grid from files, writing grid to files, initializing, computing flow and accumulation. Access the grid information (things like rows, columns, element[i][j]) via getters and setters. This way, if at any point in time you decide to go from a 2D array of values to a 1D array of values, the only part of the code that will change is the getter and setter that get/set the data value at (i,j).

Below is a possible outline of the main function, just to give you an idea. Feel free to change to match your preference and style (for e.g. you may use classes instead of plain C structures).

 int main(char** args, int argc) {


      char *elevname, *fdname, *faname; 
     //check args and argc to make sure the user entered  files
     //names on the command line, and if so, figure out the  file
     //names


     Grid elevgrid; 
     / * read the elevation grid from file into this
     structure. Note that first you have to read the number of rows
     and columns, then you have to allocate the 2D-array in the grid,
     then you have to fill it with values from the file
     */


     Grid FDgrid, FAgrid;
   

     //compute the FD grid; note that this needs to initialize the grid,
     //and allocate the data
     computeFD(elevgrid, &FDGrid);

      //compute the FA grid; note that this needs to initialize the grid,
      //and allocate the data
     computeFA(elevgrid, &FAgrid);

      //save grids to file
      gridtoFile(FDgrid, fdname);
      gridtoFile(FAgrid, faname);


     //cleanup memory
     freeGridData(elevgrid);
     freeGridData(FDgrid);
     freeGridData(FAgrid);
}

To compute the FA grid, a good idea is to separate the code into a function computeFlow(elevgrid, FDgrid, i, j) that computes FA for a specific cell (i,j) in the grid; and a function that computes FA for the whole grid by calling the previous function in a loop for all points. To compute FA, you will first initialize the FA to 1.

//initialize flow grid 
for (i=0; i < rows; i++) 
   for (j=0; j < cols; j++) 
     if point is not NODATA then set(FAgrid, i,j, 1);
     else set(FAgrid, i,j,NODATA);

Here I assumed a function on a grid set(Grid grid, int i,int j, float val) that sets the value of an element (i,j) in grid to the desired value.

Two things to think about when writing computeFlow(elevgrid, FDgrid, i,j):

Infinite recursion: if there is a bug in your flow direction logic that allows cycles, you'll get an infinite loop or infinite recursion. Basically, make sure that flowsInto() cannot create cycle. Special attention to NODATA and boundary cases (on the edge).
Efficiency: Try running your algorithm on the large grids. If you end up implementing a quadratic algorithm, you'll feel it. A good algorithm for computing the FA grid will be sub-quadratic.

Testing your output

A couple of small grids are provided to help with the debugging.

Here is how the grid test2.asc looks like:

Here is what the flow should look like:

Note that point (4, 2) has height 3, and receives flow from all its neighbors, as it is the lowest neighbor for each of them. So 2 + 3 + 14 + 2 + 3 = 24. Add to this the initial 1 unit of flow at point (4, 2) and you got 25.

If you got test2.asc to run correctly, then your program is almost correct (except maybe nodata values, and larger grids).

Test grids

A variety of real life grids are available here.

You are encouraged to search and download your own grids. Grid elevation data is available for the entire world. It would be nice if you picked an area that you are interested in and you want to learn more about.

Here is a good description of elevation products from USGS.
A site where you can download data is CGIAR-CSI site. The CGIAR-CSI GeoPortal provides SRTM 90m Digital Elevation Data for the entire world. The CGIAR-CSI SRTM dataset has undergone post-processing of the NASA data to fill in the nodata voids through interpolation techniques. The result is seamless, complete coverage of elevation for the globe.
In 2010, SRTM (Shuttle Radar Topography Mission) flew and collected data for the whoel world. This is referred to as the SRTM dataset. Data comes in tiles. The resolution of the data may vary from 30-arc seconds (1 kilometer horizontal spacing) to 1 arc-second (30m). You can download SRTM data from here. Check out an image sample. In 2014, the White House announced that the highest-resolution topographic data generated from NASA's Shuttle Radar Topography Mission (SRTM) in 2000 was to be released globally in the coming years.
Here is the EarthExplorer from USGS: Earth Explorer.

Submitting your work

Make a folder called flow in your svn folder on microwave, and update it with your work.

Enjoy!