1.2. Point operations

1.2.1. Introduction

Before we continue, type the following:

exec(open("openmaps.py").read())

It reads all maps and assigns them to corresponding variable names. For instance, the file phreatic.map is assigned to the variable phreaticMap, or topo.map becomes topoMap. If you close Python and restart it again, you need to retype this command to get all your maps loaded!!

New maps can be derived from existing maps with PCRaster functions and operators. These provde a very powerfull calculator for maps. It offers you a wide range of operators (for instance plus, minus, slope, windowaverage ) that act upon input maps. If needed, the operators can be combined, in the same way as it is done in algebraic calculations. For instance the statement:

resultMap = slope(topoMap) * (1 + phreaticMap)  <Enter>

Is perfectly valid. Try it! And display the resultMap with Aguila.

The software contains point operations, area operations, neighbourhood operations and map operations. This section covers point operations. These operate only on the values of the map layers relating to each cell; in other words, for each cell the operation is independent of the cell value(s) of neighbouring cells.

1.2.2. Quantitative computation with scalar maps

Quantitative calculation includes several mathematical operations, such as * (multiply), sin, ln, +. The input maps of a quantitative computation must always be maps of scalar data type.

To create an overlay which contains the depth of the unsaturated zone, you can subtract for every raster cell the phreatic level (i.e., the groundwater level, phreaticMap) from the surface level (topoMap). Both maps give the level in metres above sea level. Type:

unsatMap = topoMap - phreaticMap <Enter>

Display the result.

Calculate a map that gives the depth of the unsaturated zone in mm. Use the operator * and call the result map unsatmmMap. Display the map.

Somebody wants to determine a map with the infiltration rate on basis of the soil type (soilsMap) and the unsaturated zone (unsatMap), by typing the following operation:

infilMap = soilsMap * unsatMap <Enter>

Try the operation and note the error message (particularly the last line).

Question: Why is it not possible to multiply soilsMap with unsatMap?

  1. The unsatMap contains negative cell values, which does not make sense.

  2. The data type of soilsMap is nominal, and nominal data types do not allow multiplication.

  3. The data type of unsatMap is nominal, and nominal data types do not allow multiplication.

  4. The unsatMap contains very high cell values, which does not make sense.

Correct answers: b.

Feedback: The data type of the soils map is nominal, because it contains classes. A class cannot be input to a multiplication, for instance one cannot multiply a peat soil by twenty. This is why the error message is generated.

1.2.3. Boolean algebra operations

Boolean algebra is a combination technique that assumes that cells contain only two different values. Operations for boolean algebra can only be applied on maps with the related Boolean data type and the result will be a map of Boolean data type. On such a map, cells are said to be TRUE (a cell value of 1) if a certain attribute is located in that cell, and they are said to be FALSE (a cell value of 0) if the attribute is not located in that cell. Maps containing these conditions can be combined with boolean algebra operations.

You have maps containing the location of the roads (isroadMap) and the location of the surface water (iswaterMap). Display isroadMap and iswaterMap.

If we are interested in the location of bridges, we can assume that bridges occur where water coincides with roads. This is calculated with the & (and) operator as shown in the table below.

isroadMap

&

False

True

iswaterMap

False

False

False

True

False

True

Table: Boolean values in bold provide result of the & operator with two Boolean inputs, isroadMap and iswaterMap, for all combinations of True and False inputs on these maps.

The operator and can be used to evaluate this boolean relation. Type:

isbridgeMap = isroadMap & iswaterMap <Enter>

Display the inputs and results in one statement and check the location of the bridges:

aguila(isroadMap, iswaterMap, isbridgeMap) <Enter>

Other Boolean operators are ~ (not), | (or) and ^ (xor). Type:

x1Map = isroadMap & iswaterMap <Enter>
x2Map = isroadMap | iswaterMap <Enter>
x3Map = isroadMap ^ iswaterMap <Enter>
x4Map = isroadMap & ~ iswaterMap <Enter>

And display the resulting maps.

Listed below are four cross tables.

Table A

isroadMap

False

True

iswaterMap

False

False

False

True

False

True

Table B

isroadMap

False

True

iswaterMap

False

False

True

True

False

False

Table C

isroadMap

False

True

iswaterMap

False

False

True

True

True

True

Table D

isroadMap

False

True

iswaterMap

False

False

True

True

True

False

Question: Which of the cross tables belongs to x1Map = isroadMap and iswaterMap?

  1. Table A

  2. Table B

  3. Table C

  4. Table D

Correct answers: a.

Feedback: None

Question: Which of the cross tables belongs to x2Map = isroadMap | iswaterMap?

  1. Table A

  2. Table B

  3. Table C

  4. Table D

Correct answers: c.

Feedback: None

Question: Which of the cross tables belongs to x3Map = isroadMap ^ iswaterMap?

  1. Table A

  2. Table B

  3. Table C

  4. Table D

Correct answers: d.

Feedback: None

Question: Which of the cross tables belongs to x4Map = isroadMap & ~ iswaterMap?

  1. Table A

  2. Table B

  3. Table C

  4. Table D

Correct answers: b.

Feedback: None

1.2.4. Comparison operators resulting in boolean maps

For each cell, the comparison operators define a relation between the cell value on a first input map and a second input map. If this relation holds, a boolean TRUE is assigned to the cell. If it does not hold, a Boolean FALSE is assigned to the cell. The resulting map is a Boolean map. The two input maps may be real PCRaster maps or a constant value representing a map that is totally filled with cells of that value.

Try:

highMap = topoMap > 40 <Enter>

Display topoMap`` and highMap``.

The syntax of a PCRaster operator defines how the operator is used (for instance the number of input maps, the order in which the input maps must be typed in, the use of brackets). For instance, the syntax of the > (greater than) operator is: Result = expression1 > expression2. Result is the output map that is generated, expression1 and expression2 are the two input “maps”. The name expression is used here instead of a map because it does not always need to be a map name. The expressions may be:

  • names of variables (e.g. topoMap)

  • constant values (like 40 in the example given above), or;

  • operations resulting in maps.

An example of using an operation resulting in a map is:

testMap = topoMap > (0.234 * 19)

In this case Result is testMap, expression1 is topoMap and expression2 is an operation resulting in a map: (0.234 * 19), 0.234 multiplied with 19.

In addition to the > (greater than) operator, you can use == (equal), >= (greater than or equal), <= (less than or equal), < (less than) and != (not equal). These operators have a syntax that corresponds with the > operator.

Display buildgMap and click on the map to find the cell value used to represent the mine (if you don’t see the legend, it is code 5 in the legend). Make a boolean map, call it isMineMap, that is TRUE at cells with the mine and FALSE at cells without the mine. Use the buildgMap and the == operator. Display isMineMap together with the buildgMap.

Give the command you typed to get isMineMap.

Question: What command did you use to get isMineMap.

  1. isMineMap = buildgMap = 5

  2. isMineMap = buildgMap == mine

  3. isMineMap = buildgMap == 5

  4. isMineMap = buildgMap = mine

Correct answers: c.

Feedback: Note that in PCRaster you cannot use names for legend entries in operations. This is why the constructs with mine do not work. You have to refer to the cell codes associated with classes, always. Mine has a code 5 on the map, and to select mine, you need to select cells with code 5.

1.2.5. Conditional operators with a boolean map

The conditional operators use an input map of Boolean data type to determine whether a first expression, a missing value or a second expression must be assigned to a particular cell.

There are two conditional operators. The first conditional operator is ifthen, with the following syntax: Result = ifthen(condition,expression) where Result is the output map, condition is a Boolean map or expression and expression is the expression that is assigned to Result in those cells where condition has a TRUE value (cell value 1). For those cells where condition is FALSE (cell value 0) the value of Result is not defined and a missing value is assigned.

The elevation at the mine can be determined with isMineMap and topoMap. The isMineMap was created during the previous exercise. Type:

topatminMap = ifthen(isMineMap, topoMap) <Enter>

Display topoMap, isMineMap and topatminMap. The map isMineMap is TRUE (cell value 1) at the mine and the resulting map topatminMap has the value of topoMap for that area. The map isMineMap is FALSE (cell value 0) for all the remaining cells and topatminMap is assigned a missing value for these cells.

The second conditional operator is the if then else operator, with the following syntax: Result = ifthenelse(condition,expression1,expression2) where Result is the output map, condition is a boolean expression, expression1 is the expression that is assigned to Result in those cells where condition has a TRUE value (cell value 1) and expression2 is the expression that is assigned to Result in those cells where condition has a FALSE value (cell value 0)

In the future, the mine on isMineMap will be used for open-cast mining (Dutch: ‘dagbouw’): 20 metres of ground will be removed at the mine. As a result the surface level at the mine will decrease with 20 metres compared to the current surface level (given on topoMap).

Calculate the elevation map of the whole study area (no missing values) that will result after digging down 20 metres of ground at the mine. Call this new elevation map toponewMap. Use the ifthenelse operator.

To answer the following question, calculate the lowest elevation value on toponewMap:

lowestPointOnTopoNew = mapminimum(toponewMap)

Question: What is the lowest elevation value on toponewMap?

  1. 1.0

  2. -1.0

  3. 3.0

  4. 12.0

Correct answers: a.

Feedback:

toponewMap = ifthenelse(isMineMap, topoMap-20,topoMap)
lowestPointOnTopoNew = mapminimum(toponewMap)