Raingauge Network Optimization and Gis: A Case Study of the Mananga Basin

(1)

Report no. 58

October 1993

Raingauge Network

Optimization and Gis

A Case Study of the Mananga Basin

Communications of the Sanitary Engineering &Water Management Division

M.J. vanDijk/T.H.M.Rientjes

.~

..

Mananga basin Rapp

T

U

Delft

CT

Faculty of Civil Engineering

Sanitary Engineering&Water Management Division Hydrology Section

(2)

RAINGAUGE NETWORK OPTIMIZATION AND GIS

A CASE STUDY OF THE MANANGA BASIN

CT Mananga basin

.

October 1993

M.l. van Dijk & T.H.M. Rientjes

., . ecornsche UniversiteitDelft • bhotheek Faculteit der Civiele Tsch .

(

I nl9.

kaores Stevinweg 1)

P tbus 5048 .

":!

eoo

GA L· . /

Delft University of Technology Faculty of Civil Engineering

(3)

1 Introduetion 1

1.1 The Mananga basin . . . .. 1

1.2 Rainfall analyses. . . .. 1

1.3 The Mananga river 1 1.4 Cases 1 1.5 Ilwis 2 1.5.1 Input module . . . .. 2 1.5.2 Vector module . . . .. 2 1.5.3 Raster module . . . .. 2 1.5.4 Tables module . . . .. 3 1.5.5 Points module 3 1.5.6 Output module . . . .. 3 1.5.7 Command option 3 1.5.8 Dos option 3 2 Creating the digital elevation model (DEM) . . . .. 3

2.1 Rasterize 4 2.2 Interpolation 5 2.3 Display&Store .. . . .. 5 2.4 Pixel information 6 2.5 Classification . . . .. 7 2.6 Colors . . . .. 8

3 Creating the catchment map . . . .. 9

4 Creating the river map . . . .. 10

5 Digital elevation modelwith river at catchment level 11 6 Yearly rainfall of the Mananga basin . . . .. 15

6.1 The semivariogram 17

6.2 Kriging 18

6.3 Detrending Kriging 19

7 Network optimization 22

(4)

A casestudy of the Mananga Basin

1 Introduction

1.1 The Mananga basin

The study area of this case study is the river basin of the Mananga river, one of the main rivers on Cebu-island, the Philippines. Today the Mananga is used as discharge river, in the near future this river may be used for drinking water purposes of the region around Cebu-city, some 25 kilometres east of the Mananga catchment. This catchment area comprises about 102 square kilometres, and is surrounded by a chain of sharp-edged mountains which reach to a maximum height of 700 meter. The total length of the Mananga-river is 24.5 kilometres, the drainage area upstream Camp-4 is about 64km2 (see the front page figure). According to the Köppen classification the elimate of Cebu is classified as TC; tropical characterized by high temperatures and humidity, and heavy rainfall. There are, however, pronounced regional and, in most parts, seasonal variations. Average annual rainfall ranges from 1300 to 1700 mm; average annual temperatures ranges from 26°C to 29°C.

1.2 Rainfall analyses

The characteristics of the rainfall-processes of the basin can be analyzed with the raingauges in the basin. However, in the Mananga basin there are only four operational raingauge stations, a too small number to describe the characteristics from. Therefore, eight other stations lying in the surrounding mountain ranges are included in the analyses of the rainfall process.

The yearly rainfall distribution over the basin exists through a variety of weather patterns inducted by different oceanic and continental air streams. Dominating weather patterns are the a) Inter Tropical Convergence Zone (!TeZ), b) easterly waves, c) tropical cyclones and, d) convective eells.

1.3 The Mananga river

One can classify the discharge pattern of the Mananga river as a braided stream (anastomosed river) . This means that the stream is divided into an interlacing system of channels, which sediment contain a heavy load due to accelerated erosion. The rock situation in the bed itself, prevents the river from incision and promotes lateral current. Inspite of the relatively young geological age of the Philippines, the width of the Mananga riverbed is considerable (more than 100 m in some parts). In those parts, the river can develop and really does develop the channel system. During high discharges, all these channels will be filled with water and form one complete water body.

1.4 Cases

The Mananga basin will be used to demonstrate the construction of a Digital Elevation Model, some applications of geostatistics, and raingauge network optimization in combination with GIS. Use is made of the GIS software package ILWIS.

(5)

- an overall rainfall map (annual values) will be made with the use of detrending Kriging, - an optimization of the raingauge network will be performed.

The raster maps of the different phenomena will be made on a grid with a grid spacing of 100*100 meter, resulting in a raster of 106 columns and 136 rows (pixel size = 10000 m').

1.5 Dwis

llwis (the Integrated Land and Water Information System) is a GIS that integrates image processing capabilities, tabular databases and conventional GIS characteristics. The program was developed by the ITC (International Institute for Aerospace Survey andEarth Sciences) at Enschede in 1985. The conceptualization of the system takes into account that not all GIS users have a thorough knowledge of computers. Therefore all operations are performed through a user friendly menu. Experienced users can, on the other hand, perform operations directly through commands and/or command files.

Figure 1 presents an overview of the llwis main menu, this menu is subdivided in several modules and submodules. A short overview of the menus will be presented.

ILWIS GEOGRAPHIC INFORMATION SYSTEM INPUT POINTS Dos VECTOR OUTPUT Return RASTER User Change TABLES Command

Figure 1: Ilwis main menu 1.5.1 Input module

The input module enables the user to gather geographical information (in vector or raster form) and attribute data. This module also enables the user to convert files with other formats to files in llwis format.

1.5.2 Vectór module

The vector module enables the user to digitize, display, update and rasterize vector data. A raster to vector conversion option and a stream ordering analysis facility are provided.

1.5.3 Raster module

The raster module enables the user to process, analyze and visualize the geographical information stored in raster format. The raster module also enables the user to georeference raster maps and link raster maps with the attribute data stored in the internal llwis database. The raster module consists of the sub-modules Visualization, Spatial Modelling and Image Processing.

(6)

visualizaiion

This module enables the visualization and storage of raster maps. Spatial modelling

This module enables interactive modelling and analysis using one or more raster maps and/or tables.

Image Processing

This module contains all the standard digital image processing facilities and some extra statistical features.Itcao be used to process remotely sensed data, but also data derived from other sourees cao be processed digitally

1.5.4 Tables module

The tabIe manipulation module cao be used to manipulate attribute data, i.e. the non-spatial information used in a geographic information system. The module allows the user to create database queries and perform standard database operations.

1.5.5 Points module

The points module provides facilities to convert point data to vector or raster maps. 1.5.6 Output module

The output module provides possibilities to plot or print maps, graphs and tables. Itis also possible to create, display or update annotation such as legends, symbols and text.

1.5.7 Command option

This option proceeds to the command line. Inthe command line mode llwis cao be executed without intervenienee of menus.

1.5.8 Dos option

This option ca1ls Dos and allows the user to execute Dos commands. Return to llwis by typing Exit.'

2 Creating tbe digital elevation model (DEM)

With a digitizer the hundred meters isolines of an analogue topographic map (1:50.000) are digitized and stored in three vector files:

altitude. crd(coordinates file)

altitude.seg (file with line segments)

(7)

same altitudes are connected and are called segments. These segments have to be converted into a raster model. A few procedures must be performed before a nice DEM is constructed; a sequence of the procedures that will be discussed are;

- Rasterize, - Interpolation, - Display&Store, - PixelInfo, - Classification, - Colorlut.

The first procedure is to rasterize the points of the different segments. After this estimations of the altitudes have to be made for pixels where no value is present, this will be done by interpolating these values out of the known altitude values.

2.1 Rasterize

When converting a segment map into a raster map the pixel size and the coordinate system are user defined (in practice the coordinate system of the topographic map is adopted). Also the input- and output-filenames are asked. The output map, a raster map, contains two files; a MaP Data file with extension .MPD, and a MaP Information file with extension .MP!. Go to theVector module and make a raster map of the segments with the following scheme:

VECTOR - RASTERIZE - SegmentToRaster .. Segment map: altitude

• create an attribute map: y

.. enter values manual: y

.. enter the missing values, 100 for _100, etc. • Raster map: altitude

• coordinates oom existing map: n

• meters per pixel: 100

• Enter minimum and maximum X and Y: LowerLeft Corner

X: 584500 Y: 1139000

Upper Right Corner

X: 595000

Y: 1152500

name of segment file to he converted accept altitude as attribute values

type an attribute value for each code name of output file

create raster coordinate system pixel size defined.

Grid minand max.

A small map with the isolines is displayed on the high resolution screen. The map

altitude.mpd is a map with altitude values on the pixels nearest to the isolines (figure 2).

If you have a sub-menu currently active you can exit it either by moving to the 'return' option in the menu and pressing return, or by pressing the 'escape' [ESC] key.

(8)

Segments user defined grid

400

···,

···r····

···

r

· ····

··

r ·

····

··i···

····r··

·

· ··

·r

··

·

· ·1

·

· ··

·

· ~

soo

j

l

l[

1-1l1

·

· ·j

···

···:·

····

··

··1

··

· ····

··

· :

···

··

··t··

· ···

·i

···

···:···

···

·t

·

· ···

1 .

:

:::

r

:

~:

:I::

~

::

:

r-

· ~

· ···

r

·

· ;·

· ··

t

···?

·

· ::.

::'

1 :

:.

::1

:

J

:

.~

~

~"

.

•••••

• .

...

.

.l

ll

r···~

·

· T·

..:.

···

!

,

:

I1

:

,

'

::=

:

..

.

:.

.

t-

...

~

.

..L.

.

Figure

2:

Rasterize procedure

2.2 Interpolation

The second step is to interpolate the isolines(IOQ meters altitude lines) to obtain an altitude va1ue at every pixel. The interpolation will be done with the procedure FromIsolines. This procedure performs a linear interpolation; for each pixel the distances to two nearest isolines are calculated. Then a linear interpolation is performed for the unknown pixel. Because we have created a raster map we willieave the Vector module and go to the Raster module.

RASTER - SPATIAL MODELLING - INTERPOLATION - FromIsolines - input map: altitude

- four times enter - output map: altlOO

name of input raster file

accept number of rows and columns output raster file

The raster map altlOO can be displayed on the high resolution screen, this is performed in the Visualization module.

2.3 Display&Store

RASTER

=-

VISUALIZATION

=-

Display&Store

- input map name: altlOO - 4times enter

- output map name:$

- stretching: n ~enlarge factor: 4

high resolution screen

useDoemalattribute values inflate map four times

(9)

What we see is only a white mass with some colors at the bottom! This has to do with the way in which the high resolution screen is controlled. The high resolution screen can only display pixels with values within the I byte range (integer values 0-255). Every screen pixel will obtain a byte value. Our altitude map has attribute values (altitudes) which range between 95 and 805 meters, so the high resolution screen only recognises pixels with altitude values with a maximum of 255 meters. Pixels with values higher than 255 willautomatically obtain the largest possible screen value, namely 255. The white color of most of the pixels depends on the color given to every possible byte value. There are 256 byte values, so also 256 possible colorluts. The default color for the byte value 255 is white, this causes an almost complete white screen.

When for stretching the answer would be y, stretching of the altitude values will be performed. A histogram of the altitude values will be calculated and stored on disk. The raster map will now be stretched according to the 1 and 99% bounds of the histogram. Values within the 1 byte range will be assigned to the altitude values (see example below). The values of file altlOO range between 95,0 and 805,0. This means that pixels with values 95,0 obtain values 0, and pixels with values 805,0 become 255. All other pixels

will obtain new recalculated values ranging between 0 and 255. Altitudes 95,0 805,0

-New Values

o

255

For example a pixel with an altitude of 300,0 meters will get a value of 65, 400 becomes 103, 787 becomes 252, etc.

For the enlargement factor we chose the value 4, this means that the image will be displayed 4 times larger than the original size (106*136). The high resolution screen has aresolution of 800*600 screen pixels (the resolution depends on thetype of screen; VGA, SVGA, etc.). So the largest possible displayed map cao have 800 columns and 600 rows. Ifthe columns or rows of a raster map are larger, the enlargement factor should be chosen smaller than 1.

In our case we have a raster map with 106 columns and 136 rows, this makesitpossible to enlarge the "image to obtain a better view. Enlargement factor 4 means that for every pixel of the raster map 16 screen pixels will be taken (4*4).

2.4 Pixel information

There is an option to retrieve information of the pixels. Choose the procedure PixelInfo in the Visualization (Raster) menu. We see in the top left corner three columns, the fust two show the row and column number of the pixel which is depicted on the screen, the third column shows the byte value of the pixel. Below these values the X-coordinate and Y-coordinate of the pixel are displayed. When we want to see more information of the pixels, choose R (raster map). Type altlOO when the raster map is asked. Choose N when a table is asked (we don't have a table with values). Press enter if a second raster map is asked, if information of another (second or third) raster map is wanted the name of this map should

(10)

-Acasestudyof the Mananga Basin

be entered. When moving to a pixel, the available information of this pixel is displayed. Note the difference between the pixel values from the top left corner and the raster map values column. The primer values are the byte values (range from 0 to 255).

It is also possible to view values in raster maps on the monochrome screen:

RASTER - VISUALIZATION - ViewValues - map name: alllOO(.mpd)

- select window: 4times enter

- OutputDevice: screen

2.5 Classification

name of input raster file

accept number of rows and columns monochrome screen

The problem that raster values greater than 255 can not be shown on the screen (except with stretching) can be solved by classification of the attribute values. Therefore the altitude values have to be classified (indexed) and provided with a color. ILWIS has a 'color set' with color index values ranging 0 to 255. Every index value can be provided with a color, build by the user in the 'color set' option. Ifa raster map pixel has a value of 20, ILWIS will

look in the active color palette for the color assigned at the index value 20. If the value of a pixel in the raster map goes beyond the 1 byte range (integer value between 0-255), ILWIS will assign a color to this pixel with index number 255 (default value for 255 is white). Since most of the altitude values are higher than 255, the resulting raster map is mainly white. To systematize this we will introduce altitude classes, in such a way that the raster map has values withinthe 1 byte range.

Before a map can be classified we must make a classification table, thiswill be done in the Table module.

TABLES - SpecialTables - ClassifyTable - input classify table name: altclass

- create this tabie: y

The proposed classification table is:

Bound Class 50 0 100 1 150 2 200 3 250 4 300 5 350 6 400 7 450 8

table doesn't exist yet

Bound Class 500 9 550 10 600 11 650 12 700 13 750 14 800 15 850 16 2000 17

(11)

The first column of the table represent the bounds of the raster values. The second column represent the new raster values. This means that altitudes with a value between 50 and 100 belong to class 1, values between 100 and 150 to class 2, etc.. An example of the classify procedure is given in figure 3.

·

· ··

r~~;

··

r~;;

· ·T··~~

·~

··

· T·

·;~

·~

·

· r·

·~

·~~·

·

· r···

.

---r

~t~t~t~I~~

t

-I

510

I

527

i

563

i

593

I

603

I

r

~~~r~

r

;~~] ~

r

~~r-···· ··· -r- ··· ··· ··· ·· t··· · ···· ··· ··· ·t···t· ···_· ·_·· · _~··_···_···· ···· ~··· · · ·

.

_

~

_

~

L. ~

t

~~ ~

.~

.L .

.

..

.

..

.

...

..

.

..

.

~ ~

..l.

~.~

.L

~

.~

_

~.~

..

.

!

10

i

10

I

11

I

11

i

12 1

J

l~

r

~~I ;;I ;~

r

;;r

Figure 3: Classification

Now we can classify the map altloowith this classification tabie.

RASTER - SPATlAL MODEWNG - CALCULATION

- a1l101:=clfy(altlOO,altclass); create classified raster map

clfy is a Ilwis function which links a raster map to a classification tabie. The new map altl01 wi1l be calculated by classifying mapaltl00 with the classification tablealtclass. Notice that fust the name of the raster map should be given and second the name of the classification tabie.

Dur new raster map altl01 is now ready to display on the high resolution screen. 2.6 Colors

To assign colors to the different classes (index values) the Colorlut module will be used. In

this module one can choose between a standard palette, a predefined palette or a user defined palette. The·changingof colors is an interactive operation you can perform while looking at an image at the high resolution screen. For the altitudes a color palette is already defined.

RASTER - VISUALIZATION - Colorlut - Laad

- color file: altcolor (,col) file with color palette, 18 colors defined

The colors of the raster altitude map will now change from green to red in ascending altitude. With the classification and color procedures we achieved the following:

(12)

A case study of the Mananga Basin

Before classification After classification (indexed) value index value index

0-50 0-50 0

o

{black} 50-100 50-100 1 1 {green} 100-150 100-150 2 2 { 4 } 150-200 150-200 3 3 { 4 } 200-250 200-250 4 4 {

a }

250-300 250-255 5 5 { • } 300-350 255 {white} 6 6 {

a }

350-400 255 {white} 7 7 {

a }

400-450 255 {white} 8 8 {yeUow} 450-500 255 {white} 9 9 { 4 } 500-550 255 {white} 10 10 {

a }

550-600 255 {white} 11 11 { 4 } 600-650 255 {white} 12 12 { 4 } 650-700 255 {white} 13 13 { U } 700-750 255 {white} 14 14 {

a }

750-800 255 {white} 15 15 {

a }

800-850 255 {white} 16 16 {red } 850-2000 255 {white} 17 17 {blue}

3 Creating the catchment map

The next map we will create is a raster map of the catchment. With a digitizer the coordinates of the boundary (watershed) of the Mananga basin are digitized and stored in three vector files:

boundary.crd boundary.seg

boundary.slg

With the watershed segments in these files we will make a catchment raster map in order to reduce the full sized raster maps to only catchment area raster maps. For this purpose we will make a map with values 1 for points inside the basin and values 0 for points outside the basin. To make such a map we fust have to construct a polygon from the segments.

VECTOR - DISPLAY &CHANGE - Segments

- Segment file name: boundary

1----...;;..---POLYGONIZE

.. closed polygon:y

- polygon name: boundary

.. mask: >Ic

only one polygon of the total catchment name output file

take all segments

We will display the polygon to show the shape of the Mananga basin. Display fust the altitude map of the Mananga basin altlOl (this will be the coordinate reference basis of the screen), next the polygon will be displayed on top of this altitude raster map.

(13)

VECTOR""" DISPLAY&CHANGE - - Polygons

- - - -..

-Po-I-yg-o-n-ma-p-n-a-m-e-:-bO-u-nd-a-ry----I---,;...;;..---Display Polygons .. Window name: $

... use existing coordinate system: y ... mapping unit name: *

- clear window:y

high resolution screen only if a raster map is visible all mapping units

clear high resolution screen

To combine the DEM with the catchment, the catchment polygon has to be rasterized. The pixels covering the polygon on the screen will be given the values 1,the other pixels will get values O. For this we stay in the vector module and make a raster file of the polygon. This is done the same way as with segments, the difference is that we use thePolygonToRaster

module.

VECTOR""" RASTERIZE

=-

PolygonToRaster .. polygon file :boundary

.. use polygons :n

.. create an attribute map:y .. enter values manual: y .. enter value: 1

.. raster file: boundary

.. coordinates oom existing map: y - map name: altitude

Use an existing reference system

reference system of the altitude raster map

The created raster map has to be reformatted to maintain a map with the values 0 and 1. The just created map only has values for points within the catchment, for points outside the catchment no value is defined. In theRastermodule one can make calculations with different maps, also conditional statements can be made. We want pixels only to have the values 0 or 1, so we make the conditional statement that if values inside the catchment map are higher than zero (they have a value), they must take the value 1, else they have the value O.

RASTER - SPATIAL MODELLING - Calculation

.. bound1OO: =ij(boundary> 0,1,0),' ,

-Finally we have created the file we wanted, boundlOO (.mpd, .mpi). All the other files we created to realize this map can be deleted from the directory, we will not use them anymore. Display the catchment raster map on the screen and look at some of the values.

4 Creating the river map

As alike as making the catchment map, the coordinates of the river are digitized and stored in vector format. With the segments of the river we will make a raster map of the river. Display the DEM and display the river segments on top of this DEM (use thevector module with segments). The name of the river segments file is River.

(14)

A casestudy oftheMananga Basin

VECTOR

=-

Display&Change=coSegments

=-

Display Segments

The screen shows a raster map of the altitudes overlaid with a vector map of the river. The segments of the river on the screen will he rasterized in order to merge this river map with other maps.

VECTOR - RASTERIZE - SegmentsToRaster - segment file: river

- create an attribute map: y

.. enter values manual: y

.. enter value: 1()()()

.. raster file: river

.. coordinates from existing map:y

.. map name: altitude

all river pixels have value 1000

use anexisting reference system

reference system of the altitude raster map

To demonstrate this raster map of the river, display it on the screen in the raster module.

Figure 4: Raster map of the river, with the watershed displayed in vector mode

5 Digital elevation model with river at catchment level

With our created raster maps of the altitudes, catchment and river a nice representation of these topographical features of the Mananga basin will be made. Our basic altitude map

altlOOwith the altitude at all pixels will be used to start with.

First we will see what the altitudes are in the catchment; we are not interested anymore in

values outside the Mananga basin. The only thing we have to do is multiply the altitude raster map (altlOO) with the boundary raster map (boundlOO). Because the values of the boundary map are 1 inside the catchment and 0 outside the catchment,the resulting map will

(15)

be the altitude values inside the catchment. After multiplication this raster map has to be classified.

RASTER - - SPATIAL MODELLING - Calculation

... altI02:=altlOO*boundl00;

... alt103: =clfy(all102,allclass);

altitude at catchment level classified map at catchment level

The next map we will make is a combined raster map of the river and the altitudes, these two maps will be merged. The second map is glued to the fust one. Pixel values of the altitude map (fust map) will be assigned to all those pixels of the river map (second map) where no value is defined (for an integer map), or where a value zero exists (for a byte map).

RASTER - SPATIAL MODELLING - GEOREFERENCE - - Merge .. first map : altIOO

.. 4 times enter to accept defaults ... add borders: n

... second map: river

... 4 times enter to accept defaults ... output map :altI04

... enlarge factor: 1

... ignorezero values:y

altitude raster map accept columns and rows

river raster map

accept columns and rows merged raster map

This map will now be classified using the classify table altclass, and reduced to catchment sca1e using the file boundlïû.

RASTER - SPATIAL MODELLING - Calculation

- altI05:=clfy(allI04*boundl00,allclass);

I

Display this map (altlO5) on the high resolution screen.

Our final step is to add a legend to this map with altitudes. Therefore we go to the output module and make a legend.

Legend table

Text Fill Text Fill

so-iooe

1 SOO-SSOm 10 1OO-lSOm 2 SSO-6oom 11 lSO-ZOOm 3 600-6S0m 1Z Z()()..25Om 4 6S0-700m 13 25o-3OOm S 700-7S0m 14 3OO-3S0m 6 7S0-Soom IS 3S0-4oom 7 SOO-SSOm 16 4oo-4S0m 8 river 17 450-S00m 9

(16)

A casestudy of the Mananga Basin - enter window name: $

OUTPUT

=-

Annotation

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _I-hi-

·-g-h-reso--Iu-ti-o-n-sc-ree-n---Legend

.. enter input legend filename: (enter) no legend exists yet Edit - Title - title: altitude - X: 598000 - Y: 115{)()()() .. size: 5 zojust:M.C. • outline: no - (escape) text of legend justification

go10the legendmainmenu Boxes - X: 598000 - Y: 1145000 - size: 4,4 -just: M.C. - nrcol: 1 - outline:yes .. edit boxes size of boxes justification number of columns enter legend table

Enter the legend table in the edit mode. Move to the column named

Text,

press

ins

to insert

-a row -and type the fust row of the legend t-able (50-loom). Move to Fill -and enter v-alue 1 for the class assigned 10 this text line. Continue this procedure forallclasses. After this press escape, save the legend as altleg, and display the legend on the screen (figure 5).

Altitude

o

50-1_

o

10&-1SOoo

o

15e-Z_

o

ZG&-Z5e. 0%50-3_

o

38&-35e. El35&-4_ Q4e&-45e. 450-500. se&-55e. . ~ .688-65e. • &50-7_ .708-75e.

.750-8_

.8G&-8SOoo . . .1_..

(17)

What we have accomplished so far is the generation of a few maps, altitude ,river and catchment. A resume of the maps is given below:

Altitude raster maps

.. altitude : isolines map of the altitudes

- altlOO : interpolated map of altitude isolines

- alslûl :classified map ofaltlOO

.. alt102 : altitudes in the mananga basin

.. altlD3 : classified map ofaltlD2

.. altl04 : altitudes and river

.. altlDS : classified map ofaltl04in the mananga basin Catchment raster maps

.. boundlOO : catchment map with values 0 and 1 River raster maps

- river : river map with values of 1000

After these fust exercises it is assumed one can rasterize vector maps, merge two maps, classify these maps, use colors and display the maps on the screen.

(18)

6 Yearly rainfaU of the Mananga basin

With a Geographical Information System we can undertake some statistical operations on data sets. There is also some mathematical modelling capacity, among which regression modeis. This tooI will be used to estimate the correlation between altitude and precipitation in the Mananga basin.

For this exercise two data sets are available:

- the mean annual precipitation values at 12 raingauge stations in and around the Mananga basin,

- the altitudes of those 12 stations.

In the previous exercises we constructed already a raster map containing the altitude values in the Mananga basin, named altl02.

The goal of this exercise is to determine the relation between altitude and precipitation, and to construct, with this information, a mean annual rainfall map of the Mananga basin. Data from 12 raingauges are present to establish this relation.

Station Altitude Rainfall Station Altitude Rainfall

(m) (mm) (m) (mm) 1 395 1666 7 250 1345 2 225 1514 8 154 1146 3 220 1695 9 730 1589 4 190 1412 10 460 1636 5 420 1534 11 540 1478 6 370 1399 12 50 1371

From these data we can in genera! say, "the higher the raingauge, the higher the rainfall value". However, we need a more precise technique for specifying this relation. In statistics a procedureca1ledregression analysiswiU do this, for which llwis provides a module named Statistics. This procedure analyses the relation between two raster maps or two data sets. The data set with the altitudes and precipitation data is stored in a (tabie) file

stations.

dat. The rainfall valuesZ(x) can now be thought of as the sum of two components, a deterministic component '(trend, a_o

+

aI*A(x) calculated with linear regression, and a variability

component (a(x):

Z(x) - Qo+al*A(x) + cx(x)

where _aoand al are the parameters of the altitude effect, A(x) are the altitudes, and a(x) are the residuals.

(19)

TABLES - Table Calculation - Enter Table name: stations.dat

I

Special Features - Statistics- Least Squares - Select X-axis: altit&

- Select Y-axis:prec&

- order: 1

- Display:y

.. Window name: $

- Clear window: Y

.. 7 times enter

altitudevalues are the independent part precipitation values are thedependent part

Polynomial

linear regression

display on high resolution screen

accept defaults

After running this regression procedure we see that there is a distinct relation between the altitude and the rainfall values, we have a a_oof 1346.77 mm, and a al of 0.455 mm/m.

2000 1800 1600

•

• _•

1400

---

_•

1200

_•

1000 SOO 600 400 200 0 0 200 400 600 800

Figure 6: Linear regression

With the parameters _aoand ah and the altitude valuesA(X) the deterministic component of the rainfall values are calculated. The differences between the deterministic rainfall values and the meàsured values, the residuals atx) can be used to estimate the variability of the rainfall process in space (variability around the deterministic trend). This variability is caused by scale effects, it can be decomposed in a large-scale, a micro-scale, and measurement errors. The assumption is made that the closer the stations, the more rainfall values become alike; some kind of correlation is expected between the rainfall values depending only on the distance, the lag vector h, between two measurement stations. The basic theory by which such relations are explained comes from Geostatistics. For this exerciseitmust be understood that the rainfall values are not purely random; rainfall values consists of a certain spatial relation.

(20)

c

6.1 The semivariogram

When a graph is made of the distance between two stations on the X-axis and the rainfall varianee differences between those stations on the Y-axis, the relation becomes clear. Such a graph is called a semivariogram. Different mathematical models are made to fit the experimental semivariogram (i.e. the calculated semivariances from the measuring points).

In figure 7 a spherical model is drawn, with three parameters determining the shape of it. Co is the nugget effect, representing the random varianee, a varianee caused by measuring errors and micro-scale variations (variations over a distanee smaller than the average distanee between measuring stations). a Is the range, representing the maximum distance when two stations still have some correlation. The stations with a distance larger than the range are assumed to be independent. The maximum semivarianee of the phenomenon, reached at the range distanee, is called the sill Co

+

C.

y(h) c:.+c ""..~. -,

I

i i

i

Co...•...

!

~

!

:

a h

Figure

7:

Spherical model

In the following table the measured rainfall values, the deterministic (altitude depending) components, and the variability (residuals) are presented.

Station number Measured value Deterministic comp. Variability comp.

(ao+a1*A(x) (a(x) 1 1666 1527 139 2. 1514 1499 15 3 1695 1445 250 4 1412 1433 -21 5 1534 1538 -4 6 1399 1515 -116 7 1345 1461 -116 8 1146 1417 -271 9 1589 1679 -90 10 1636 1556 80 11 1478 1593 -115 12 1371 1370 1

A semivariogram of the residuals is shown in figure 8. The model of the semivariogram, fitted through the calculated semivariances, is used in geostatistics to show the spatial

(21)

...

-

..

-A casestudy of the Mananga Basin

25000 20000 15000 ••• • •••• • ••••••• •••••••• • • • ••••••;;.0•••..-...-...----1

.

... 10000 5000 o 0 12000 24000 36000

Figure 8: Semivariogram, Co=3000,

C=17000, a=2()()()()

dependency of a phenomenon. This spatial dependency on its turn is used for interpolation with the Kriging method.

6.2 Kriging

The Kriging method gives estimates of a phenomenon at non-measured 1ocations with a prediction of the reliability of the estimation. In our case we can calculate the expected rainfal1 values on a grid of 100*100 meters, with for each estimate a standard deviation representing the accuracy of the estimate.

Kriging is a 'linear' type of estimator, that is , it gives a weighted average from the measuring points. This estimator is denoted byZ and is equal to:

"

Z· - À_{1* Z1} +~*~ + .. +À"*Z" -

L

Àj*Z,

i-i

where the À1,ÀZ' ••'À"are the weights assigned to each measuring point.

The weights are calcu1ated using the semivariogram, this makes that the assigned weights depend on the distance between an estimated point and a measuring point, and the distance between the measuring points themse1ves, this to prevent the overweighing of c1ustered data. An example of the distribution of the weights is given in figure 9.

0.096 .0 0.307

~

\..

..

·0.001

.

...

..

.

..

.

...

....

...

.

e

. . . ...

:~~~~~~~~~~~~::::::::

::

(8)0.159 ,1.Cl'.I -0.035 point -(8)

e

0.474

Figure 9: weights with ordinary Kriging

(22)

6.3 Detrending Kriging

In the original Kriging equations, there is a restrietion for the phenomenon under study; no trend is allowed. In our case we assume the existence of a trend of the yearly precipitation depending on the altitude values. We have to subtract the trend from the original values and interpolate the residuals over the whole domain. In this case we speak about detrending (or residual) Kriging.

With the spherical model of the semivariogram, Kriging estimates and their standard deviation are calculated and stored in a table file:detr.pnt. This file containing four columns (X-coordinates, Y-eoordinates, estimates of the rainfall variabilities, and standard deviations of the variability estimates) cao be read by ILWIS in the Point module.

POINT - PointToRaster - Table name: detr.pnt

- Use column of table as attr. value?:y

- Select attribute column:prec&

- Raster File: detr

- Copy transf. from existing map?:y - Enter map name: altitude

variability estimates output raster map name

use reference system of this file

This raster file must be classified with the underlying table (detrprec.clt), fitted to the watershed and stored in the file detrprec.

Bound Class Bound Class

-50.0000 1 75.0000 6 -25.0000 2 100.0000 7 -0.0001 3 125.0000 8 0.0001 0 150.0000 9 25.000 4 175.0000 10 50.000 5

Show the map of the interpolated residuals(detrprec.mpd) with the color filedetrprec. coland legend detrprec.leg on the high resolution screen (figure 10).

With these residuals the rainfall values at every pixel cao be calculated using the altitude map and the parameters _aoand al.

Z(x) - 1347+0.455

*

A(x) +a(x)

In the calculation function of the spatial modelling module, the new rainfall map, calculated from the altitudes and the residuals can be made. The following function can be applied:

RASTER - SPATIAL MODELLING =-Calculation precyear: =if(altl02>0, 1347 +0.455*altl02 +detr, 0)

(23)

yeal'ly I'alnfall De~rendlng Krlglng residuals

D

-?S - so -Q -sa - -25 _

D

-25 - e_

a

9 - 2 5 _

a

25- 5 9 _

a

59- ? S -m!! ?S - 199 _

lil

199 - 125

-11125 -

159 -.159 - 175 _

Figure 10: Residuals interpolated with detrending Kriging

performed, otherwise the rainfall values should be set to zero. The rainfall raster map can be c1assified with table precyear. clt and displayed with legend precyear.leg, The c1assification table

precyear.

cltlooks like this:

Bound Class Bound Class Bound Class

1450 0 1575 5 1675 9

1475 1 1600 6 1700 10

1500 2 1625 7 1725 11

1525 3 1650 8 1750 12

1550 4

We now have a map incorporating the mean annual rainfall values with altitude effect over the Mananga basin (figure 11).

Together with the interpolated rainfall values a plot of the variances of the rainfall estimates can be made; this consists of the Kriging standard deviations, a measure of accuracy of the estimates. A pixel with a high value for the standard deviation represent a pixel where the estimate must be interpreted with caution. These standard deviations are stored in the fourth column of the file detr.pnt. From this column a raster map must be made using the same

procedures as with the variability estimates, store this map in a file (detrvar).

This raster file must be c1assified with table detrvar.clt, fitted to the watershed and stored in the file varrain.

(24)

A case study of the Mananga Basin yeal'ly I'alnfall Detl'endlng Krlglng

o

145&-1475 _

o

1475-1500 _ EI1500-1525 _ EI1525-1550 _

o

1SS&-1575 _

m

1575-1680 _ 8 1fo09-1625 _ 8 162S-165O _ lil165&-1675 _ 1116?5-1700 _ • 1700-1725 _ • 1725-1750 _

Figure 11: Rainfall values calculated with deirending Kriging

The classification table detrvar. cltlooks like this:

Bound Class Bound Class Bound Class

70.0 0 90.0 4 110.0 8

75.0 1 95.0 5 115.0 9

80.0 2 100.0 6 120.0 10

85.0 3 105.0 7 125.0 11

A legend is made, named varrain.leg. Show the map of the Kriging standard deviations

(varrain.mpd) with the colorfile varrain.col and legend on the high resolution screen.

~

;:-.'t~'_

Figure 12: Kriging estimation varianees

Estl.atlon UaI'lances (standal'd devlation) D7075 -IE7598 -1illI98-95 _ Irn8590

-tm

9095 -II1II95-188 •• • 198-185 _

(25)

From this image (figure 12) one can clearly see the locations where the raingauges are stationed, this by the much lower standard deviations. The Northern part of the area, without raingauges, is less trustworthy than the middle part, resulting in an error standard deviation of above 100 mmo

7 Network optimization

Itis clear that the estirriated varianee is an important factor for determining the location of raingauges, the varianee depends only on the geometrical location of the raingauges. Obviously, the choice of the variogram model and its parameters is conditioned by the particular set of available data. But once the variogram model is chosen, the varianee can be viewed as depending exclusively on the location of the raingauges. Hence it becomes possible to compute the error varianee associated with any set of hypothetica1 data points without getting actual data at these points.

The estimation varianee is thus a very suitable tool for network optimization. The existing network of raingauges can be replaced or supplemented by the best representative new set of raingauges at specified locations.

In our case the standard deviations are highest in the northern part of the area. The best places to add a new raingauge would be on that particular place where the estimation varianee is highest. If two stations are added, the estimation varianee in the northern part is reduced considerable (figure 13).

Estl.atlon Uarlances Cstandard deviatlDn) [ ] 78-75 _ 1illI75-Be _ lEIBe-lIS _ fj]lIS-ge _ im98-95 _ 1195-188 •• • 18e-185_

Figure

13:

Kriging estimation variances with two additional stations

The two added stations are only known by coordinates from the rainfall map. A topographic map should be analyzed to check if it is possible to place the new stations at or close to the chosen locations in order to satisfy other criteria like reachability. For this exercise we have to our disposal a road map of the Mananga basin from which the distances to the new locations can be calculated. The road map is stored in a file

road.mpd.

From this me a distance map will be calculated and classified:

(26)

Acasestudy ofthe Mananga Basin

RASTER=-SPATlAL MODELLING=- Distanees .. Compute also Thiessen map:n

.. Input map name: road

~ Output distanee map name: Dist

.- scale [0]: enter only integer values

Figure 14 shows a distance map with a vector model of the watershed and a legend. Classes of distances are: 0 =>100, 100=> 250, 250=> 500, 500 =>1000, and 1000=>00 meters.

If the distanee map and the varianee map (without extra stations) are combined, the best locations for additional stations become visible. The problem is now what weight to assign to the distanee and what weight to assign to the estimation varianee for optimization. This is very subjective and will be different for every hydrologist.

dlctance classes 0&-188 • Elt...-Z58 • [Jzso..;w • mlSGO-t_ • • >

1_.

Figure

14:

Distances

from

the

roods in classes

Make a map which calculates the best place for additional raingauges. An example is given in figure 15, which has been made by the calculation map:

=

(maxdist-distt/Iû

+

detrvar5,*

and classified. In the distanee map the maximum distanee is 5885 meters.

(27)

Figure 15: Optimallocations for additional raingauges

Possible loca~ions

fo~ new ~aingauges

o~ 0 .... [] noNlu. minoNlu.-ce.cN

.c:aeN

• ~ c:aeN

(28)

8 References

Chua, S.H., Bras, R.L.

Optima! estimators of mean area! precipitation in regions of orographic influence Journal of Hydrology, vol 57 no 1/2, May 1982

Delhomme, J.P.

Kriging in the hydrosciences.

Advances in Water Resources, vol 1 no 5, 1978.

Dijk, M.J. van, Kappel, R.R. van

Optimization of the Cebu rainfall network.

M.sc Thesis, Hydrology section, Technical University of Delft, 1992. Dijk, M.J. van, Rientjes, T.H.M.

Geostatistics and Hydrology, student syllabus, 1993,in publication. International Institute for Aerospace Survey andEarth Sciences

ILWIS 1.3, User's manual, volumes 1+2, may 1992.

International Institute for Aerospace Survey andEarth Sciences ILWIS 1.3, computer program, may 1992.

(29)

In the series "Communication of tbe Department of Sanitary Engineering and Water Management"are edited:

1. Sieben, H.H.: Patterns and variability of phosphate and heavy metals in sediments of two shallow lakes.

(November 1985)

2. Flipse, M.J. en Heide, J. van der:Ontwikkelingen met betrekking tot vaste afvalstoffen ex art.4, 17, 25, 26

van de Afvalstoffenwetin periode van ca.1980 tot 1985.

3. Kop, J.H.: Planvorming voor de drinkwatervoorziening. (februari 1986)

4. Blanken, J.G. den en Hoogb, M.P.A.J. de: Modellen voor desinfectie van gezuiverd afvalwater met chloor en

ozon.

5. Kop, J.H.: Het probleem van de wederzijdse afstemming van de belangen van drinkwatervoorziening en

milieubescherming bij de planning voor de winning van zoet grondwater. (augustus 1986)

6. Boekelman, R.H. enNiet,H. de: Het berekenen van modelkrommen voor Geo-elektrlsche metingen.

7. Vos, W.L., Donze, M. and Buiteveld, H.: On the reflectance spectrum of algae in water: the nature of the

peakat 700 nm and its shift with varying algal concentration.(October 1986)

8. Smit, D., Mameren, H.J. van en Veldkamp, R.G.: De zuurstofhuishouding van de Utrechtse Vecht.

(november 1986)

9. Heide, J. van der: Kinetische modellen voor ontwerp en beheer van actief-slib-installaties deel 1 en 2.

(februari 1987)

10. 8oulan, R.P., Donze, M. en Klapwijk Sj.P.: Fosfaatbalans van de polder Reeuwijk en een aantal

deelgebieden.(maart 1987)

11. Groot, C.P.M. de en Breemen, A.N. van: Ontspanningsflotatie en de bereiding van drinkwater. (juli 1987)

12. Blanken, J.G. den en Hoogb, M.P.A.J. de: Modelvorming voor verwijdering van indicatororganismen in het

actief-slibproces.

13. Misbra, K.K. and Breemen, A.N. van:Gravel-bed flocculation.

14. Vlis, E. van der: De filtratietheorie.(maart 1988)

IS. Koreman, E.A. en Breemen, A.N. van: Toepassing van het vriesdooiproces bij de ontwatering van

coagulatieslib. (mei 1988)

16. Ganzevles, P.P.G., Kop, J.H. en Ywema, R.:Materiaalkeuzeafvalwaterleidingen. (juni 1988)

17. Nieuwenbuyze, R.F. van, Stokman, G.N.M., Kuijper, R., Gerritsen, J.J. en Donze, M.: Detectie van

proceswater met behulp van thermischeremote-sensing, (juni 1988)

18. Blanken, J.G. den en Hoogb, M.P.A,J. de: Modelvorming voor een goede procesregeling van de desinfectie

met chloor c.q. ozon aan de hand van instelbareenlofdirect meetbare variabelen.(augustus 1988)

19. Noppeney, R.M.:De invloed van stagnante zones op dispersie. (november 1988)

20. Noppeney, R.M.: Gevoeligheidsonderzoek Alarmmodel Rijn; De invloedslengte van samenvloeiingen bij

dispersie.(november 1988)

21. Noppeney, R.M.: De verspreiding van olie op rivieren benaderd met het Taylor-rnodel.(november 1988)

(30)

24. Blanken, J.G. den:Afscheidssymposium prof.ir, A.C.J. Koot. (januari 1989)

25. Hoojkaas, LJ., Donze, M. en Klapwijk, Sj.P. : Fosfaatbalans van de polder Reeuwijk en de Reeuwijkse plassen.(januari 1989)

26. Verwoerdt, P. en Mazijk, A. van: De één-dimensionaledispersievergelijking van Taylor bij een opdeling van de rivier in vakken.(maart 1989)

27. Mazijk, A. van:GevoeligheidsonderzoekAlarmmodel Rijn;eindrapportage. (mei 1989)

28. Blanken, J.G. den en Hoogb, M.P.A.J. de:Desinfectie van behandeld afvalwater met chloor: vergelijking van eenpunts-en tweepuntsdosering;

deel 1:Tekst,bijlage A, B en C.

deel 2:Bijlage D,E, Fen G. (mei 1989)

29A. Verstappen, G.G.C.: Gedrag van organische micro-verontreinigingen in rivieren. (juli 1989)

298. Mooren, J.J.M. en Heide, J. van der: Leaching of heavy metals from thermally decontaminated soils.

(maart 1989)

30. Nieuwstad, ThJ., Wortel, N.C., Bout, F.N. van den en Alting, B.J.: Een vergelijking tussen ladingsgewijze en continue zuivering van afvalwater.(juni 1989)

31. Kramer, J.P., Wouters, J.W. en Kop, J.H.: Dynasand Filtratie.(July 1989)

32. Nieuwstad, ThJ.: Treatment of municipal wastewater in a pilot-scale airlift-loop reactor.(December 1989) 33. Ankum, P.: Polders; achtergronden,ontwerp en toekomstige ontwikkelingen.(juni 1990)

34. Brandsma, T.: Evaporation in Hydrology and Meteorology. (July 1990)

35. Mooren, J.J.M.: Het uitlooggedrag van kunstmatig samengesteldeen verontreinigde grond.(2 delen) (augustus 1990)

36. Singb, S.N., Boekelman, R.H., Rientjes, T.H.M. en Dam, J.C. van: Behaviour of groundwater of the polder Groot-Mijdrecht.

37. Boekelman, R.H. en Rientjes, T.H.M.: Workshop hydrological models.(September 1990)

38. Stavrides, N., Rientjes, T.H.M. en Dam, J.C. van: Network optimization, a simple approach applying GIS and MLR.

39. Duindam, P., Morales, C. and Heide, J. van der: Investigacion sobre los desechos solidos de laciudad de Masaya, Nicaragua.(enero 1991)

40. Heide, J. van der: Evaluacion hidraulica de plantas potabilizadoras de filtracion rapida en Nicaragua. (junio 1991)

41. Heide, J. van der: Metodologia de potabilizacionde agua superficial en Nicaragua.(abril 1992) 42. Veldkamp, R.G.: Randvoorzieningen van rioolstelsels kritisch beschouwd. (oktober 1991)

43. Kruithor, J.C., Schippers, J.C. aod Dijk, J.C. vao: Abastecimiento de Agua Potable de Agua Superficial en los aiios Noventa. (febrero 1992)

44. Rietveld, L.C. aod Matsinbe, N.P.: Pilot Plant Studies on Slow-Sand-Filtration and Up-Flow-Roughing Filtration in Mozambique,(January 1993)

(31)

airlift-loop reactor. (June 1993)

48. Branclsma, T.: Sewer Systems and Climate Change.(September 1993)

49. Ankum, P. and Brouwer, R.: Nevengeulen en Sedimentverdeling.(November 1993)

50. Brouwer, R.: Review of the water management systems in the Gujarat Medium Irrigation 11 Project. (August 1993)

51. Hoornstra, J.S. and Jong, J. de: Water quality management in the Netherlands. (October 1993) 52. Jong, J. de: Thirdpolicy document on water management (June 1993)

53. Jong, J. de and Hoornstra, J.S.: Analyse der anforderungen im EG-bereich (October 1993)

54. Jong, J. de, Baarsma, J.P., Frintrop, P.C.M. and Mulder, W.H.: Screening strategies for organic micropollutants in context with water pollution control. (October 1993)

55. Jong, J. de: The method and mechanisms of establishing consensus on water management policy in the Netherlands. (October 1993)

56. Dijk, M..J. van and Rienijes, T.H.M.: Geostatistics and Hydrology; part 1:Spatial and Temporal Variability. (September 1993)

57. Dijk, M..J. van and Rienijes, T.H.M.: Geostatistics and Hydrology; part 2: Estimation Techniques.(November 1993)

(32)

· Bibliotheek TU Delft

Fac. CiTG subfac. Civiele Techniek