6 minute readRainwater harvesting (RWH): potential sites

In the first part of this RWH series I explained what kind of techniques there were for rainwater harvesting. In this second part I will show you a basic method using RSAGA to get a first idea of suitable locations. This is based solely on the digital elevation model, or DEM, which is a digital representation of a terrain’s surface. There are only a few steps to this method, of which the fist is data collection.

Finding open sourced data

I am a huge fan of open source data, as it gives research and projects more possibilities (because it is free). A downside can be the quality of programmes or data, but usually there is a big community to fix bugs, and commercial enterprises are opening their sources as well (which guarantees some quality at least).

Digital elevation model (DEM)

Finding data can be difficult, luckily this article that can help. The article shares five global DEM sources. Although the data is free, the resolution may not be suitable in all cases. One of their suggestions is from the United States Geological Survey (USGS). There I downloaded a DEM for an area in the south of Chile (tile name: ASTGTM2_S41W073), with a resolution of 30 meters (so each pixel is 30×30 meters).

Here you see where the tile is situated in Chile (near Valdivia/Osorno)

Example digital elevation model (DEM) in Chile
Example digital elevation model (DEM) in Chile

The DEM downloaded from USGS looks like this (white areas have higher elevation than black):

Close-up of the digital elevaiton model (DEM)
Close-up of the digital elevaiton model (DEM)

And here you see the coloured version. There is a mix of mountains, lakes, rivers and lowlands.

Coloured version of the study area
Coloured version of the study area



All the steps in making the eventual map of all potential locations for RWH are done in R, with some packages in addition. These are mentioned first, after which we will discuss the different processes.


We will use the following packages:

  1. RSAGA is a package providing access to geoprocessing and terrain analysis functions of SAGA-GIS (System for Automated Geospatial Analysis; Conrad et al., 20151), and R functions for spatial data manipulation and applying functions to stacks of grids;
  2. Rgdal is used to import the DEM file; and
  3. The raster package is used for reading, writing, manipulating, analysing and modelling of gridded spatial data.

Filling the sinks in the DEM

The DEM data may include some ‘sinks’, or depressions in the terrain that could influence the hydrology. Water may not continue to flow downstream once it flows into a sink. So, this is the first step: using rsaga.fill.sinks. Fill in the in.dem, out.dem and the method (= “wang.liu.2006“).

Using RSAGA means files are no longer .tif files but .sdat. To still be able to read them use plot(raster(readGDAL(“filledDEM.sdat”))). This is the result:

This is the filled DEM
This is the filled DEM

With the raster package it is possible to subtract one DEM from another, of which the result is:

The green spot is a volcano which has been filled up. There seem to be quite some sinks in the original DEM. The lakes have also been filled up.

Contributing area

The next step is to find out where pixels drain their water into, and we need to know where the pictures received the water from (contributing area). Using rsaga.topdown.processing will calculate this from the highest to the lowest cells.

The next two processes are calculating the total and specific catchment areas, which is in effect the contributing area. For each cell is calculated how many cells drain into it.

(The images are not that exciting so I did not post them, but don’t forget to check out what they look like for yourselves!)

Topographical wetness index (TWI)

Now that we know where water will flow, it is nice to know where saturation may occur. This is where the topographical wetness index (TWI) comes in, which is another way of representing the accumulation of water. Higher values indicate drainage depressions while lower values represent crests and ridges. This is what the map looks like:

Topographical wetness index (TWI) of the study area
Topographical wetness index (TWI) of the study area

With green areas being wetter and white being dryer. You can see that in the mountains there are a few spots that are wetter, the lakes are visible and the flatter area on the left has many depressions.

If you want to put rainwater harvesting measures on all green spots you’d need a lot of time, so let’s decide which areas are more important.


First, lets see the histogram of this map to see what TWI values are dominant:

histogram of TWI

It appears there is one value that is extremely high (6.5). This is probably the lakes. We don’t want to implement anything in the lake, so we will reclassify this value to lake. The values below 0 are unsuitable, leaving the values above the lake value for suitable. We can split these up into main (from 99th quartile up) and secondary priority (from lake value to 99th quartile):

RWH locations based on TWI (from DEM)

Although this is a crude method, it does give some basic insights into where you can start. Field validation of these locations will show if they are really suitable locations, or not. Other methods also take other variables into account such as land use, or distance from the road, which will further determine the suitable locations. The number of red pixels in the above image will then decrease of course, making it easier to do field validation.

The resolution of the data is 30 meters, meaning that one pixel has quite some area. A next step could be converting this raster dataset to polygon and then selecting the biggest polygons. These will be the sites with the biggest impact for example. But this all depends on what you as a user want of course.

Rounding up

This post shows a quick way of determining some potential rainwater harvesting locations, based solely on the elevation. There are other methods, which we will maybe explore in the future. The criteria for the different techniques needs to be taken into account then (maximum slope for example).

For these maps I used Rstudio, but there are other programmes available, such as QGIS and ArcGIS (paid service), which produce similar maps with the same method.

Conrad O, Bechtel B, Bock M, et al. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. G. 2015;8(7):1991-2007. doi:10.5194/gmd-8-1991-2015

Add a Comment

Your email address will not be published. Required fields are marked *