Terrain elevation models are used to model and understand a wealth of natural processes in geosciences and beyond. Like in the previous labs, in this lab you will work with grid terrains. You will extend your GridGIS by adding functionality to compute the flow direction grid, and to compute and classify the flat areas of the terrain.
In this lab you'll build on the previous methods and write code to compute the flow direction of each point in a grid. Naturally, this defines a grid called the FD grid. The FD at a point is defined as the direction of its steepest downslope neighbor. The neighbors are coded as follows:
6 7 8 5 x 1 4 3 2
That is, if the FD of the cell marked with an x is towards its East neighbor, its flow direction should be assigned to 1; if it's toward its SE neighbor, its FD is 2; and so on.
When all neighbors of p have elevation greater or equal to the elevation of p, the FD of p cannot be set. This includes two cases: when p is part of a pit, and when p is part of a plateau area.
To summarize, the FD of a point p in the grid is defined as follows:
Your first task is to write a method computeFlowDirection() which computes and returns a FD grid.
Grid computeFlowDirection()As a helper to this function, you will write a function that computes the FD for a specified position (i,j), and you will call it in a loop for all i and j.
For example, for test2.asc:
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 5 The FD grid is: 2 2 3 4 4 2 2 3 4 4 2 2 3 4 4 ...
Once you compute the FD grid, try to visualize it using your method that renders elevation grids.
Note that the 0 values in the FD grid represent the flat areas in the terrain. Remember the definition of blobs from class --- in this case we are interested in blobs of 0's, not of 1's. The first step is to identify the flat areas, and classify them.
We distinguish two types of flat areas: plateaus and sinks.
|
|
|
Your second task is to classify each flat area as being either a plateau or sink, and re-label each point with FD=0 accordingly.
final static int SINK = -1; final static int PLATEAU = -2; /* this function re-labels all the points with FD=0 to FD=SINK or FD=PLATEAU, depending on whether the respective blob (flat area) is a sink or a plateau */ void labelFlat(Grid fd)You'll break this into a couple of smaller functions. The first one is findBlobType.
/*
When this function is called, the point (i,j) must be inside the grid
and it must be part of a blob, i.e. its FD=0. This function finds the
type of flat area and returns its type (SINK or PLATEAU). As a side-effect, part of blob may be marked.
*/
int findBlobType(Grid fd, Grid elevation, i,j) {
if (i,j) is outside the grid or if it's FD is not 0 ==> error
//the blob containing (i,j) is a sink if ALL the neighbors of the
//blob which are not part of the blob (i.e. their FD is a value 1 through 8) have
//elevations higher than the blob. Otherwise it's a plateau.
mark (i,j) to avoid infinite recursion
for each neighbor (i+k, j+l) of (i,j)
if the neighbor is outside the grid, skip it
if the neighbor has FD = NODATA, skip it
if the neighbor has a valid FD (a value 1 through 8)
then we found a neighbor of this blob; check its elevation,
and if elevation(i+k, j+l) is not higher than elevation(i,j),
then return PLATEAU
if the neighhbor has FD=0 and is marked, skip it
if the neighbot has FD=0 and is unmarked
call findBlobTYpe recursively on it
if this returns PLATEAU, return PLATEAU immediately
//if we are here, all neighbors are higher or none of the recursive
//calls returned a plateau. so return SINK
}
The second one is labelBlob.
/*
The blob containing (i,j) is partially marked. This function labels all the
points in the blob at position (i,j) with the type given as argument.
*/
void labelBlob(Grid fd, int i, int j, int type) {
//recursive formulation
basecases: if (i,j) is ouside grid or fd (i,j) is not 0 or marked, return
if fd(i,j) is 0 or marked, set it as type and recurse on
all its neighbors
}
So far, you have findBlobType and labelBlob which
must be called on a point (i,j). Put everything together in a function
that re-labels the points with FD=0 in the FD grid, as follows:
void labelFlat(Grid fd) {
for i = 0; i< nrows; i++)
for (j=0; j< ncols; i++)
if FD(i,j)=0 then findBlobType(fd, i,j) labelBlob(fd, i, j) with the resulting type. As a result, the entire blob containing (i,j) will be labeled.
else: this point is either not a flat area, or its been
already labeled, so continue
}
Example
The grid is: ... The FD grid: ... The re-labeled FD grid is: ...