Chapter 1

Case study GIS modelling

This case study deals with methods of combining multiple maps for empirical modelling in a GIS. Several approaches to analyze multiple maps are introduced by means of a number of basic exercises that follow closely the chapter "Tools for map analysis: multiple maps" from the book Geographic Information Systems for Geoscientists, by Graeme Bonham-Carter (Pergamon,1994, 398 p.). The various methods, in the book illustrated with a hypothetical landfill site selection problem and mineral potential mapping, will be treated in our case study with a real world data set from the Chinchina Area , Central Colombia, South America and comprise the following techniques:

Boolean logic models,

Binary evidence maps,

Index overlay with multi-class maps,

Fuzzy logic method,

Bayesian methods; weights of evidence modelling

.

For the first four of the above techniques, A landfill site selection problem is used for the first four of the above mentioned techniques. A landslide hazard zonation problem is used to explain the Bayesian methods.

1.1	Problem definition: landslide hazard zonation using weights of evidence modelling

In the methods we have used so far, we were interested in finding the optimal location of a landfill site. We used subjective weighing, based on expert opinion, so called "knowledge-driven" modelling.

In this exercise we will assign weights by statistical procedures. The weights of different classes are determined by calculating densities. The method is therefore "data-driven".

The landfill study case is not useful to demonstrate these "data-driven" techniques, since we cannot evaluate the best site for a landfill by evaluating the conditions at existing landfills (there is only one in the area).

Therefore we will use another problem for the statistical analysis: the prediction of landslides. In the statistical modelling of landslides we use the rule "the present and past are the key to the future". This means we can look at the conditions under which landslides have occurred in the past, and use the critical combinations for predicting the possible occurrence of landslides where the same conditions prevail, but which are still landslide free.

We will use the same data set as for the landfill case study, with the exception that we will only work with the maps slope, geol, landuse, and slide. Apart from these we will use the map geom, which contains a very coarse subdivision of the study area in three geomorphological units:

-	The western hills; a terrain unit with rounded hills, mostly in schists, with relatively gentle slopes,

-	The terraces; a terrain unit consisting of the terraces of the Rio Chinchina, and,

-	The Romeral fault zone; the terrain unit in the eastern part of the study area, with steep slopes in a fault zone.

The map slide contains the landslides that have occurred in the area, it is an ID map, in which each landslide has a different code. The explanation of the landslides can be found in the table slide. Information is available on the following items:

The landslide code is composed of two parts: first the landslide ID number is given, followed by a - and after that a six-digit code for [TYPE], [SUBTYPE], [ACTIVITY], [DEPTH], [VEGETATION] and [SCARP]. The explanations for these are given below:

TYPE	SUBTYPE	ACTIVIT
Y	DEPTH	VEGETATION	SCARP
1	Slide	Rotational	Stable	Shallow	Bare	Scarp
2	Flowslide	Translational	Dormant	Deep	Low vegetation	Body
3	Flow	Complex	Active	High vegetation
4	Derrumb e
5	Creep

Table: Example checklist used in photo-interpretation of mass movements.

In this exercise we will use the weights of evidence method, as explained in the book from Bonham-Carter.

1.1.1	Prior probability

Consider that you would only have the landslide map of the area, without any additional data. How to answer the question: what is the chance of having a landslide in the area?

If there is no additional information, i.e. if you do not know whether you are in a geological unit with many landslide, or on a steep slope, we speak of the prior probability.

The prior probability that a certain pixel in the map will have a landslide is simply the density of landslides in the entire area. In other words, the number of pixels with landslides, divided by the total number of pixels in the map.

In a formula that is:

where:

Pprior = P{S}	Conditional probability of having a landslide S

Npix(Slide)	Number of pixels with landslides in the map

Npix(Total)	Total number of pixels in the map

1.1.2	Conditional probability

If we have more information available than just the landslide map, we can see what will be the probability that we may have a landslide, given that we are in a certain geomorphological unit, for example. If we are in the geomorphological unit "Romeral fault zone" the probability that there may be a landslide is obviously much larger, then when we are on flat river terraces. This is what we call conditional probability. Considering the relationship between a binary variable map (B) and a landslide map (S) The probability given a certain condition. It can be expressed as:

P{S B}	the conditional probability of having a landslide while you are in unit B.

In other words, the conditional probability of having landslides, given that you are in a certain unit, is the density of landslides within that unit, calculated as the number of pixels with landslides in the unit, divided by the total number of pixels in the unit.

1.1.3	Calculation of weights

As can be seen in the text of Bonham-Carter the conditional probability of having a landslides given that we are in a certain terrain unit can also be written as:

This states that the conditional (posterior) probability of a landslide, given the presence of the terrain unit equals the prior probability of the deposit multiplied by a certain factor. This factor is treated as a weight and express to what extent the probability of landslide occurrence increases or decreases with respect to the prior probability for the slope class under consideration.

In a similar way you can derive the posterior probability of a landslide occurring given the absence of the terrain unit:

P{B}	= Npix(B)/Npix(total)

P{}	= Npix(total)-Npix(B)/Npix(total)

In this equation the multiplication factor expresses the amount to which the posterior probability increases or decreases with respect to the prior probability given the absence of that particular terrain unit.

The same model can be expressed in the so called "odds formulation" in which the prior and posterior probabilities are converted to prior and posterior odds and the weight factors are the logarithms of the multiplication factors used to express the posterior odds as the prior odds plus the weight factor. We refer for a full derivation of the odds formulation is referred to Bonham-Carters's book. Here we will directly quote from this derivation the posterior odds :

Taking the natural logarithm from both sides of these equations yields:

and

The following section explains how the weights W+ and W- can be calculated by combining the number of pixels from a particular terrain unit (B) and landslide pattern (S). Four combinations of B and S are possible in the map: B S, S ,B and .

The conditional probability of choosing a pixel with a landslide, given that the cell contains pattern B, is:

and the three other conditional probabilities are:

The positive and negative weights (Wi+ and Wi-), which combine these conditional probabilities can be calculated by combining equations 1.21 till 1.28:

and

In GIS the method can be implemented rather easily. Consider the simple crossing of a binary landslide map with a binary variable map (in our example a particular terrain unit). The four possible resulting combinations are given in the table.

Landslides	Variable class represented as binary pattern
1 (present)	0 (absent)
Present 1	Npix1	Npix2
Absent 0	Npix3	Npix4

The weights of evidence can be written in numbers of pixels as follows:

and

Although the method was originally designed for the use of binary maps, it is also possible to calculate in ILWIS the weights for all classes in a multiclass map simultaneously. Let us look at an example again. The use of a classified slope map, with 9 classes.

The cross table will look as follows:

slopec	Active	Npix
1

1

2

2

3

3

4

4

5

5

6

6

7

7

8

8

9

9	0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1	2000

10

1500

20

1200

50

1000

100

900

90

800

60

700

50

100

10

50

3

In order to calculate the weights using equations 1.31 and 1.32 you need to create the following columns:

nmap	= total number of pixels in the map

nslide	= number of pixels with landslides in the map

nclass	= number of pixels in the class

nslclass	= number of pixels with landslides in the class

The values needed for formulas 1.31 and 1.32 are:

npix1	= nslclass

npix2	= nslide-nslclass

npix3	= nclass-nslclass

npix4	= nmap-nslide-nclass+nslclass

The weights are calculated as:

W+	=log((npix1/(npix1+npix2))/(npix3/npix3+npix4))

=log((npix1*(npix3+npix4)/((npix1+npix2)*npix3)

W-	=log((npix2/(npix1+npix2))/(npix4/npix3+npix4))

=log((npix2*(npix3+npix4)/((npix1+npix2)*npix4)

Let us now calculate with ILWIS the weights for the map geom..

We will only use the landslides that are classified as active, since some of the old landslide have occurred in other climatic conditions. If they would also be used, the results would be difficult to interprete.

Because the same calculation has to be performed for the other maps, it would be useful to automate the procedure. This can be done with a SCRIPT, i.e. a text file in which the exact procedure followed in the program and the formulas are stored.

You can find all operations that you have done with ILWIS in the so-called ILWIS.LOG file. The commands from this file can be copied to a script file in ILWIS.

When the script runs without problems, it can also be used for the other maps.

1.1.4	Calculating contrast factors

To quantify the spatial association between a map class and the occurrence of landslides, the contrast factor as mentioned in Bonham-Carter (1996) is defined:

The contrast factor is 0, when the landslide pattern and map class pattern overlap only by the expected amount due to chance, positive when there is a positive association between the two patterns and negative when there is a negative association between the two patterns.

1.1.5	Combining the weights of several patterns: posterior probability

In real world situations, the occurrence of the phenomena of our interest (landslides, mineral or groundwater deposits) will seldomly be associated with one spatial variable. Landslide occurrence for example is not only dependent on slope steepness but also on the geological unit on which it occurs. Ultimately we want make a prediction based on a number of maps that are likely to be useful to predict the occurrence of land slides.

We can use the weights calculated for the individual map patterns for an assessment which maps are important for the prediction in combination with our expert knowledge. In addition we can use the calculation of contrast factors as criteria to convert the input maps to binary patterns, in which the 1's indicate areas favorable for the occurrence of the phenomena of interest and the 0's the areas unfavorable for the occurrence of the same phenomena.

The calculation of the posterior probability for landslides given the presence of several input maps, is given by the formula:

W+ becomes W- if the i-th pattern is absent instead of present. Where data is missing the weights are 0.

Up till now we have only calculated the W+ and W- separately for each class in a multiclass map. If we want to use this map in the combination, we will have to define a final weight, which for each class consists of the positive weight for the record where the class is present, plus the negative weights of all the other classes. Besides of that we have to make sure that the output tables with the weights can be linked to the maps: the tables should have an equal number of records than the units in the maps.

1.1.6	Final analysis and report writing

In the previous exercises you have learned the various steps for executing a weights of evidence analysis. However, the final map was not made by careful consideration of the input data. Nor were all the input layers used.

For the final assignment you are asked to perform weights of evidence modelling on the same data set more accurately.

The following input maps are at your disposal:

geom	the geomorphological map

geol	the geological map

landuse	the landuse map

slope	the slope map

aspect	the slope direction map

roads	a vector map with the roads

fault	a vector map with faults

drainage	a vector map containing drainage lines with different stream orders

sources	a point file containing the starting points for first order drainages.

Two landslide maps are available:

slide60	containing the landslides mapped from airphotos from the 1960's

slide80	containing the landslides mapped from airphotos from the 1980's , and checked by fieldwork

The analysis will be done on all active landslides.

The evaluation of the input data should be done by calculating the weights and contrast factors for the different classes and to calculate the most optimal buffers around linear and point features using the method presented earlier.

The calculation of weights will be done using the active landslides from the 1960's. You can test the overall goodness of prediction by combining the resulting probability map with the landslides from the 1980's.

1.2	Optional exercises: weights of evidence modelling

In the previous exercises you have used the weights of evidence method, as explained in the book from Bonham-Carter, chapter 9. However, we have not calculated the actual probability, nor tested the assumption on conditional independence, since the maps used for landslide hazard zonation, are never conditionally independent. For those interested in the procedure for testing conditional independence, and calculating posterior probabilities and testing the final map, some optional exercises are given. They are not required to be included in the report, which will only be on the topics mentioned before..

1.2.1	Optional exercise: Calculating contrast factors for continuous fields

The contrast factor can not only be applied on nominal categorical maps (maps containing classes) but also useful to evaluate maps with ordinal, interval or ratio measurement levels, such as the distance from a feature with respect to the occurrence of landslides. Bonham-Carter discusses this, for another application, on pages 319-320. Instead of calculating the contrast factor for each map class, the calculation proceeds in a cumulative form by defining classes with increasing upper limits (increasing distance) and keeping the lower limit (distance 0 from the feature) constant. The contrast factor can be plotted with increasing upper limits for the binary pattern. The maximum contrast factor on such a curve would provide the optimal threshold in converting the distance map to a binary map in which the pixels inside the class (value 1) are maximally positively associated with landslide occurrences with respect to pixels outside (value 0) that class.

In this exercise we will apply this concept for the distance to roads, considering that the presence of roads may have triggered landslides in the nearby in the area.

1.2.2	Optional exercise: Conditional probability using two data layers: conditional independence

As was explained in the chapter from Bonham-Carter, the weights of evidence method can only be used when the individual maps in the analysis are conditionally independent.

If two maps are conditionally independent, the following should be true:

This can be expressed as number of pixels:

The left hand side of equation 1.17 is the observed number of deposits occurring in the overlap region where both B1 and B2 are present.

The right hand side is the predicted number or expected number of deposits, in this overlap zone, which should equal the number of landslide pixels in B1 times the number in B2 divided by the total number of landslide pixels, if the two maps are conditionally independent.

Let us look how to solve this in ILWIS with two binary (bit) maps B1 and B2.The various parts of this formula can be calculated by first crossing the two maps B1 and B2 . The result map C is then crossed with the landslide map. From the resulting cross table the values fro B1 and B2 can be obtained by tabel joining. After joining the cross table may look like this:

C	B1	B2	S	Npix
1

2

3

4

5

6

7

8	1

1

1

1

0

0

0

0	0

0

1

1

0

0

1

1	0

1

0

1

0

1

0

1	1000

50

1000

130

800

10

1500

10

The first thing to do is to delete all combinations without landslides. Than the table looks like this:
C	B1	B2	S	Npix
2

4

6

8	1

1

0

0	0

1

0

1	1

1

1

1	50

130

10

10

By using a key on the columns B1 and B2 you can find out the number of pixels with landslides in each class of B1 and B2: NpixB1 and NpixB2.

By using a key on the column S you can find out the number of pixels with landslides: NpixS. The table then looks like this:

C	B1	B2	S	Npix	NpixB1	NpixB2	NpixS
2

4

6

8	1

1

0

0	0

1

0

1	1

1

1

1	50

130

10

10	180

180

20

20	60

140

60

140	200

200

200

200

The expected frequency of landslides can now be calculated:

Npixexp:=(NpixB1*NpixB2)/NpixS

After that check whether the two maps are conditionally independent is done with the formula:

And then compared with tabled values for 2 with one degree of freedom. These are for various probability levels:

80%	90%	95%	97.5%	99%
1.64	2.71	3.84	5.02	6.63

Degrees of freedom may be defined as the number of observations in the sample, minus the number of parameters estimated from the sample. In this case 2-1=1.

In Bonham-Carter all examples are given for binary (=bit) maps, only containing 0 and 1 values. In ILWIS we don't need to convert all our input maps to bit maps first. The procedure can also be done on multiclass maps. Suppose that B1 has two units, and B2 has 3 units.

The various parts of this formula can be calculated by first crossing the two maps B1 and B2 . The result map C is then crossed with the landslide map. From the resulting cross table the values fro B1 and B2 can be obtained by table joining. After joining the cross table may look like this:

C	B1	B2	S	Npix
1

1

2

2

3

3

4

4

5

5

6

6	1

1

1

1

1

1

2

2

2

2

2

2	1

1

2

2

3

3

1

1

2

2

3

3	0

1

0

1

0

1

0

1

0

1

0

1	500

20

1200

55

1500

75

2000

10

3000

100

2500

25

You also have to delete all records without landslides first. Then the table will look like this:

C	B1	B2	S	Npix
1

2

3

4

5

6	1

1

1

2

2

2	1

2

3

1

2

3	1

1

1

1

1

1	130

55

75

10

100

25

By using a key on the columns B1 and B2 you can find out the number of pixels with landslides in each class of B1 and B2. However, you now make two columns for B1: NpixB1cl	containing the number of pixels with landslides in each class.

NpixB1no	containing the number of pixels with landslides in all other classes. You do that by first calculating the total number of pixels with landslides in all classes: Npixtot, and subtracting NpixB1cl from Npixtot will give us NpixB1no.

Also for B2 you have to calculate the same: NpixB2cl and NpixB2no.

By using a key on the column S you can find out the number of pixels with landslides: NpixS.

The predicted number of pixels in each class you can find with the formulas:

NexB1B2:=(NpixB1cl*NpixB2cl)/NpixS

NexB1noB2:=(NpixB1cl*NpixB2no)/NpixS

NexnoB1B2:=(NpixB1no*NpixB2cl)/NpixS

NexnoB1noB2:=(NpixB1no*NpixB2no)/NpixS

In order to do the 2 test you still have to know the number of observed pixels with landslides for each combination.

NobB1B2	is the number of pixels when both classes are present. this is in fact equal to Npix.

NobB1noB2	is the number of pixels in other combinations of the same class of B1 with different combinations of B2. You can find this out by subtracting Npix from NpixB1cl

NobnoB1B2	is the number of pixels in other combinations of the same class of B2 with different combinations of B1. you can find this out by subtracting Npix from NpixB2cl

NobnoB1noB2	is the number of pixels in all other combinations of B1 and B2. You can find this out by subtracting NpixB1cl+NpixB2cl from NpixS.

Now that you know the observed pixels with landslides and the expected ones, you can execute the formula for chi-square testing as given in 1.18, with the difference that you do not sum up the records. The four different cases are now represented in columns.

Although the procedure for testing conditional independence of multiclass maps is more complicated, it works faster than subdividing all multiclass maps into binary maps first and then checking per pairs.

Now that you know the procedure for testing conditional independence, let us try it out with two examples, one using binary maps, and one using multiclass maps.

You can imagine that when you have to check many different maps it will be very boring to type in the same formulae everytime. Because the same calculation has to be performed for the other maps, it would be useful to automate the procedure. This can be done with a script, i.e. a file in which the exact procedure followed in ILWIS is stored.

1.2.3	Combining the weights of several patterns: posterior probability

The calculation of the posterior probability for landslides given the presence of several input maps, is given by the formula:

W+ becomes W- if the i-th pattern is absent instead of present. Where data is missing the weights are 0.

Up till now we have only calculated the W+ and W- separately for each class in a multiclass map. If we want to use this map in the combination, we will have to define a final weight, which for each class consists of the positive weight for the record where the class is present, plus the negative weights of all the other classes. Besides of that we have to make sure that the output tables with the weights can be linked to the maps: the Calculation of the posterior probability from the logarithmic weights factors defined in the "odds formulation" requires some additional algebraic manipulations:

To calculate the final posterior probability we need still to calculate the logarithm of the prior odd for landslide occurrence which is defined as:

Now that we have both arguments from the right hand side of the equation we can calculate the logarithm of the posterior odds with the formula:

Then you have to convert to posterior odds by applying the exponential:

And you can calculate the posterior probability using the formula:

1.2.4	Optional exercise: overall test of the final map

A simple test can be applied to determine the total number of predicted landslide pixels. The predicted number N{S}calc is determined by adding together the product of the area in pixels, N{A}, times the posterior probability , P, for all probability classes in the map, thus:

where there are k=1,2,...m unique combination in the map.

Then a 2 test can be performed:

The observed number can be found by crossing the classified probability map with the landslide map active, and calculating the number of pixels with landslides in each class. The predicted number of pixels with landslides per class can be found by multiplying the average probability value of each class with the number of pixels in each class.

Then the chi-square test is done by summating the results for all classes using formula 1.40.