[ltoma@dover:\~] flow test2.asc test2flow.ascThis will read elevation grid test2.asc, compute its flow grid and save it as test2flow.asc; unless a path is specified, files are by default in the current folder.
Below are some suggestions for how to structure the assignment. The idea is to avoid the computation of a direction grid. This is not the most efficient way, and once you understand the problem, you might decide to approach it differently. That's fine!! But I expect everybody will peruse the suggestions below in order to get familiar with the problem.
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.
Your program should have a main() function that takes the name of the elevation grid and the name of the flow grid on the command line, reads the elevation grid in memory, computes a flow grid, and saves the flow grid to a file with the given name. The suggested outline of the main function looks something like this:
int main(char** args, int argc) {
char *elevfname, *flowfname;
//check args and argc to make sure the user entered two files
//names on the command line, and if so, figure out the two 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 flowGrid;
//need to set the rows, columns in this grid and allocate its 2Darray
//compute the flow grid
computeFlow(elevGrid, &flowGrid);
//save the flowGrid to a file
}
One of the formats for grid data is the so called Arcascii format. All datasets provided above are in this format. The first 6 lines in the file represent the header of the terrain; it stores the number of rows in the grid, the number of columns, the geographical coordinates of the lower left corner, the spacing in the grid (assumed to be the same both horizontally and vertically), and a nodata value (this is the value that is used for points where the elevation could not be measured).
ncols 391 nrows 472 xllcorner 271845 yllcorner 3875415 cellsize 30 NODATA_value -9999Following the header there are nrows lines, each line containing ncols values, for a total of nrows * ncols values. These values represent the elevations sampled from the terrain. Thus, a grid terrain is basically a 2D-array of elevation values. You can assume that the elevations are integers (though you can also use floats or doubles if you want). The dimensions of the grid can be found out from the header. The other information stored in the header (xllcorner, yllcorner, cellsize, NODATA_value) you will happily ignore for this lab since we won't care about geographical coordinates---- but you have to read them anyways in order to get to the point in the file where you can start reading the elevations. The NODATA value is the elevation value assigned to points whose elevation is unknown. So a NODATA value in the elevation grid means "unknown".
To store all the information of a grid use a struct:
typedef struct _grid {
int rows, cols; // the size of the grid
....
float** data; //the 2D array of value in the grid
} Grid;
FILE* f;
char s[100];
int nrows;
f=fopen("myfile.asc", "r");
if (f== NULL) {
printf("cannot open file..");
exit(1);
}
fscanf(f, "%s", s);
printf("read %s from file\n", s);
fscanf(f, "%d", &nrows);
printf("read %d from file\n", nrows);
To test that you read the right values, write a printGrid (Grid g) method that prints the grid (header and values); it should print the same information as when you open the grid.asc file. You may want to break this function into two: printHeader and printValues.
In addition to the printGrid method, write a printInfo (Grid g) method that prints the important information about the grid: rows, cols, h_min and h_max.
See some test grids at the top of this page. Some of them (europe.asc) are large, so better not open them in a browser. While testing, you should try getting the small ones to work before you try any of the larger ones.
Knowledge about flow on the terrain turns out to be the fundamental step in many geographical processes. The location of rivers can be measured directly but it is a very expensive process. This is why researchers have looked into ways to model flow based on the elevation of the terrain. Modeling flow opens the avenue not only to knowing the location of rivers, but knowing the river and watershed hierarchy, predicting flooding and so on.
For a grid terrain, the flow direction and flow accumulation are grids of the same size as the elevation grid.
The flow direction at cell (i,j) is assigned such that:
//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(flowGrid, i,j,NODATA);
Here assume you have a set 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.
The first method that you'll write will be:
int flowsInto(Grid grid, int i, int j, int k, int l)which returns true if cell (i,j) flows into cell (k,l), that is, if cell (k, l) is the steepest downslope neighbor of cell (i,j) in grid. The way this method works is the following: it looks at all 8 neighbors of cell (i,j) in grid and computes the lowest neighbor; this is the direction where water would go from cell (i,j). Then it checks whether this neighbor is equal to (k,l), in which case it returns true; otherwise it returns false.
Once you have method flowsInto written and tested, you can start working on computing flow. You will write a method that computes a flow grid based on the (elevation) grid:
void computeFlow(Grid elevgrid, Grid * flowgrid) {
//initialize flow grid
...
//compute flow grid
for (int i=0; i< rows; i++) {
for (int j=0; j < cols; j++) {
//compute flow of point (i,j)
int flow = computeFlow(elevGrid, i, j);
set(flowGrid, i,j, flow);
}
}
}
Naturally, you will write a separate method
computeFlow(Grid elevGrid, int i, int j) that computes the
flow of a point (i,j) in the grid. There are several ways to go about
computing the flow of a point (i,j). Probably the simplest way is to
run a loop through all its neighbors, for every neighbor you will
check whether it flows into the current cell, and if it does, you will
add its value of the flow to the flow of the current cell. Note that
a cell can receive flow from multiple cells. Also note that you will
need to compute the value of the flow of the neighbor, which means
that the problem is probably recursive. When no neighbor flows into
the current cell there are no subsequent recursive calls.
Note that you do not need to specifically compute a flow direction grid. We get around it by computing on the fly flowsInto(...). This may not be the most efficient way, but it saves you from having to worry about another grid.
Two things to think about when writing computeFlow(i,j):
For debugging reasons, you should first get test2.asc to compute flow correctly. Here is how the grid test2.asc looks like:
The grid is: 9 9 9 9 9 8 7 6 7 8 7 6 5 6 7 6 5 4 5 6 5 4 3 4 5Here is what the flow should look like:
The grid is: 1 1 1 1 1 1 2 4 2 1 1 2 9 2 1 1 2 14 2 1 1 3 25 3 1Note 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).
./render2d [your-flow-grid]Note: the grid has to be in arcascii format.