lunes, 5 de noviembre de 2012

Some hydrological analysis

For this entry, we are going to do some hydrological analysis. Starting with a DEM, we are going to extract a channel network, delineate watersheds and calculate some statistics. We will use mostly algorithms from SAGA, but also from other providers, to complement them. And finally we will show how to use one of the powerful tools of SEXTANTE: the iterative execution of algorithms. Using it, we will calculate a histogram to show how elevation is distributed in each subwatershed, and plot it using R.

Let's start.

The first thing is to load this DEM.




And the first module to execute is the Catchment area (Parallel) one (use the search box in the toolbox to find it). You can use anyone of  the others named Catchment area. They have different algorithms underneath, but the results are the same.

Select the DEM in the Elevation field, and leave the default values for the rest of parameters.


This algorithm calculates many layers, but the Catchment Area one is the only one we will be using.


You can get rid of the other ones if you want.

The rendering of the layer is not very informative. To know why, you can have a look at the histogram and you will see that values are not evenly distributed (there are a few cells with very high value, those corresponding to the channel network). Calculating the logarithm of the catchment area value yields a layer that conveys much more information (you can do it using the raster calculator).



The catchment area (also known as flow accumulation), can be used to set a threshold for channel initiation. This can be done using the Channel network algorithm. Here is how you have to set it up.



Use the original catchment area layer, not the logarithm one. That one was just for rendering purposes.

If you increase the Initiation threshold value, you will get a more sparse channel network. If you decrease it, you will get a denser one. With the proposed value, this is what you get.



The image above shows just the resulting vector layer and the DEM, but there should be also a raster one with the same channel network. That raster one will be, in fact, the one we will be using.

Now, we will use the Watersheds basins algorithmtol delineate the subbasins corresponding to that channel network, using as outlet points all the junctions in it. Here is how you have to set the corresponding parameters dialog.



And this is what you will get.


This is a raster result. You can vectorise it using the Vectorising grid classes algorithm.




Now, let's move into the second part of this tutorial. Let's try to create a set of histograms that show how elevation values are distributed within each subbasin. First, let's see how to do it for a single subbasin. The idea is to have a layer that just represents the elevation within that subbasin and then pass that layer to R so we can calculate the histogram and plot it.

First, let's clip the orginal DEM with the polygon representing a subbasin. There is another SAGA algorithm to do so, and it's called Clip Grid with Polygon. Although we are calling an external application, since we are using SEXTANTE, it will be aware of the selected features in the layer we use as input. SEXTANTE will take care of creating an intermediate layer, so we do not have to worry about it and the integration between QGIS and SAGA is more complete. so, if we select a single subbasin polygon and then call the clipping algorithm, we can clip the DEM to the area covered by that polygon.

Select a polygon,




and call the clipping algorithm with the following parameters:


The element selected in the Input field is, or course, the DEM we want to clip.

You will get something like this.



This layer is ready to be used in R, but first we need an R algorithm that can use it. Add a new R script named  Raster Histogram, with the following code:


##R scripts=group
##layer=raster
##showplots
hist(as.matrix(layer))


This R script is, in fact, included as an example script, but you should try to add it yourself and get some practice. If you know R, you can even try to adapt it to plot the accumulated histogram or the values from an hypsometric curve, more commonly used to describe the relief of a watershed.

If you run the script using the clipped DEM as input, you will get a plot of the elevation histogram within the basin used to create that clipped DEM.

We can link this two last algorithm (clipping with SAGA + plotting with R) creating a simple model, as shown below. It takes a polygon layer and a raster layer, and then does the clipping and the plotting.



The intermediate clipped layer should not be a final result, but the output from the R script has to be marked as final, so you should enter a name in the corresponding box. Enter "Plots" as the name of that final output.
Save the model and name it "clipping + plotting". Now if you double click on the new entry corresponding to that model in the toolbox, you will have the following parameters dialog.


If you enter the vector layer as input, it will use the whole layer, and you will get a single histogram. You can select one of its polygons, and you will get the histogram for that basin. however, what we want is to have all the histograms corresponding to all the polygons. Here is where the iterative execution comes in.

Select the watershed vector layer as input, but click on the button on the right-hand side of the row corresponding to that parameter (the one with the green icon). That will cause SEXTANTE to run the algorithm (the model in this case), n times, one for each of the n features in the selected vector layer. That is exactly what we want.

All the resulting plots will be added to the SEXTANTE results window, as usual.


You can also create a model with the algorithms in the first part of the tutorial. That way, only two steps will be needed to get from the DEM to the set of plots: one for creating the layer with watersheds, another one to do the clipping and plotting.



6 comentarios:

  1. Hy!
    Thanks for the wonderful post!
    But, I cannot find the catchment module via search box, (in the sextante toolbox) ... Any recommandations? What am I doing wrong? Any other plugin missing?
    Thanks!

    Kind regards

    ResponderEliminar
    Respuestas
    1. Hy,

      I would like thanks to for the post, but I have the same ploblem with search tool.
      the version of sextane I'm using is 1.0.9, the last version.
      Existing other name for the tool?
      Best regards

      Eliminar
  2. DEM25.zip not found, the link is broken...

    ResponderEliminar
  3. Hi,
    I ´ve been trying to change the colours in my DTM in order to change from the greys to brighter colours as in your DTM.
    Can you explain me how to do it, please?
    Thanks.
    Ana

    ResponderEliminar
  4. Hi, Should we "fill sinks" in the DEM before doing this procedure?

    ResponderEliminar
  5. Thank you so much for this post,

    I hope I can find answer for this problem that I face win I follow the step of your post
    this massage appear:
    "input layers do not have the same grid extent"

    ResponderEliminar