Assignment 4: FD and FA

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 like this:

[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.

Test grids

Test grids are available here.

Another 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 no data voids through interpolation techniques. The result is seamless, complete coverage of elevation for the globe.

Another site is USGS geodata. This provides SRTM 30m data. It comes in tiles. From this site you can select "Digital elevation" and you'll get GTOPO30 data, which contains Global 30 Arc-Second Elevation Data Set global DEM with a horizontal grid spacing of 30 arc seconds (approximately 1 kilometer). These files might come in a different format, so if you want to use one of them, I can convert them for you. The dataset europe.asc is from this site, and contains 1km grid data for Europe.

General outline

The suggested outline of the main function looks something like this:

nt 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;
     //need to set the rows, columns in this grid and allocate its 2Darray 


     //compute the FD grid 
     computeFD(elevGrid, &FDGrid);

      //compute the FA grid 
     computeFA(elevGrid, &FAgrid);

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

}

For a grid terrain, the flow direction and flow accumulation are grids of the same size as the elevation grid.

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

The flow direction at cell (i,j) is assigned such that:

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 there are ties, they are broken arbitrarily.

To compute FA, I suggest that you 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(flowGrid, 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.

Then you would follow with something like this:

//assume FD rgrid has already been computed, and is given as argument
//assume FAgrid is allocated and its 2Darray is allocated
void computeFlow(Grid  elevgrid, Grid FDgrid, Grid * FAgrid) {


   //initialize  FA grid 	
   //set everything to 1
   ....

		
    //compute FA 
    for (int i=0; i< rows; i++) {
	for (int j=0; j < cols; j++) { 

           //compute FA of point (i,j)
	   int flow = computeFlow(elevGrid, FDgrid, i, j);

           set(FAGrid, i,j, flow);
	}
    }
}

Note that I suggest separating the code into a function computeFlow(elevGrid, FDgrid, i, j) that computes FA for a specific point (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). Feel free to move away from this outline, but you will need to argue that your end code is more versatile/readable.

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

Infinite recursion: if there is a bug in your flow direction logic that allows cycles, you'll get 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.

For debugging reasons, you should first get test2.asc to compute flow correctly. 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).

Testing your output

To test your output, you can try vizualizing your FD and FA grids using your 2d renderer.

Submitting your work

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

Enjoy!