# The Rasterscript

### The PCRASTER GIS

All input and output maps in LISEM are in the format of the PCRaster Geographical Information System. PCRaster is a free GIS is produced by PCRaster environmental software, and can be obtained through their website. There you will find a complete list of functions, research and coding examples and literature. PCRaster is not only a normal GIS but a generic programming language for dynamic landscape modelling.

PCRaster is used in two contexts here:

1. LISEM is largely built in a mixture of PCRaster commands and C++
2. The modelling capabilities of PCRaster are used to create the input maps with one or more macros (scripts).

LISEM is developed in an environment that is a mix of C++ and PCRaster libraries. A large part of the source code is written in the macro language of PCRaster The advantage of this close integration is primarily for the model developer: the PCRaster macro-language includes a powerful ‘calculator’ that is able to parse virtually any mathematical or logical equation. Equations that are commonly found in erosion models can be typed almost literally from a text to include them into the model. An example: the Manning’s equation can be written in PCRaster as a string that is given to a ‘calculator’ (the statements are translated to C++ at compilation):

V=if(n gt 0, R**(2/3)*sqrt(max(slope(DEM), 001))/n,0);

where V, R, DEM and n are raster maps, and slope is the GIS operation to calculate the gradient from the DEM but it cannot be smaller then 0.001. The instruction is carried out for all cells in the raster maps. As a precaution, the logical “if” statement sets the Velocity to 0 where n is 0. The advantage of this approach is that the code is very compact and that several developers can work with the code at the same time, without the need of knowing a higher programming language such as C++ or Pascal.

Note that for processes where iteration is required (such as the kinematic wave and Richard’s type infiltration) special functions are developed. The advantage for the user lies in the fact that the output maps of runoff and erosion/deposition can be viewed directly with PCRaster without the need of conversions.

PCRaster can also be used to vcreate input maps using separate statements for each map. These commands have the general syntax:

pcrcalc [options] out.map = in.map…

where a large number of opertaions are available for the in.map. Note that PCRaster is very strict on the variable types which have to be specified explicitely if the type type of resulting map is ambiguous. The types are:

• boolean = 0/1
• directional = 0-360
• ldd = network direction (1-9)
• nominal = class integer
• ordinal = class integer
• scalar = real value

## Understanding the PCRaster script

### binding section

The top of the script has the binding section, which links input and output variables in the script with actual files and tables. Comments in the script are in grey and start with a “#”. The syntax is: “script var name = file name”. The input variables are as follows:

The unitbase.tbl table (below) with soil parameters is a simple matrix of values. In the script, the table gets the variable name soiltbl. This area has 4 map units to which the soil hydrology is linked: two types of soils – clay loam and heavy clay (units 2 and 4), and the same soil type but with grass and shrubs causing a higher infiltration (units 1 and 3). Unit 5 is a unit that can be used for grass strips, assuming a higher infiltration. The row numbers have to correspond to map unit numbers, they can be any integer as long as the map unit number occurs in the table. Note that the table can have more lines than there are units in the map, so you can add information for different scenarios. The column numbers are the variable numbers and used in the script to make the maps. For instance the command:

report ksat = lookupscalar(soiltbl, 1, unitmap)*mask;

creates the map ksat by taking the values from soiltbl column 1 and row numbers from unitmap and multiply the result with map mask to ensure the proper catchment boundary.

In this script the map is created from crossing the basic soil units with a the vegetation cover map to cerate the 4 basic units. At the top of the script is the option “#! –matrixtable” which assures that this table is read correctly.

 0 1 2 3 4 5 6 1 32 0.51 30 0.3 1.5 0.1 2 6.3 0.55 30 0.37 1 0.05 3 1.4 0.52 40 0.39 1.5 0.1 4 1 0.52 40 0.39 1 0.05 5 40 0.51 30 0.3 1.5 0.1

Below the input file in the binding section is the list of output files. The output filenames are the default filenames used in openLISEM. You can change these but then you have to change these also in the interface and runfile of openLISEM.

### initial section

The following section of the script is the initial section, where all the maps are created and reported. The operation report means that the result of a command is written to disk, using the name as specified in the binding section. This part is one long list of commands, separated for visibility by comments.

The first maps are the DEM derivatives used in the overland flow and kinematic wave.

 ### BASE MAPS ### chanm = if(drains > 0, 1, 0)*mask; # give all channels a value 1 barriersc = if(chanm > 0, 0, barriers); # no barrier when channel, this is a culvert report Ldd = lddcreate(DEM-out*10-chanm*2+barriersc, 1e20,1e20,1e20,1e20); # runoff flow network based on dem, # main outlet, channels and barriers report outlet = pit(Ldd); # should be the same now as mainoutlet.map !!! report grad = max(sin(atan(slope((DEM+barriers)*mask))), 0.005); # sine gradient (-), make sure slope > 0.005 report ups.map = accuflux(Ldd,1); # for display, a color map with the flow network

Next are the maps that define the land cover, infiltration parameters and soil surface characteristics. Note that the LAI is the leaf area index (used for interception) of the plants in the gridcell, not the average of plants and bare soil the gridcell. If there is one mature tree in an otherwise bare gridcell the LAI has a value of for instance 6 m2/m2. Cover and LAI are used in the interception and splash detachment modules.

 ### LAND COVER MAPS ### report coverc = veg_cover*mask; # fraction plant soil cover, rfom satellite img lai = ln(1-min(coverc, 0.99))/-0.4; # LAI of plants inside gridcell (m2/m2) report lai = if(coverc gt 0, lai/coverc, 0);

The soil surface maps are used in the surface storage modules where also the flowwidth of the overland flow is determined. The roads and stony surfaces effect also the splash detachment.

 ### SOIL SURFACE MAPS### report rr = max(lookupscalar(soiltbl, 5, unitmap) * mask, 0.01); # micro relief, random roughness (=std dev in cm) mann = lookupscalar(soiltbl, 6, unitmap)*mask; # read column 6 of soil table report mann = mann * (1+0.5*housecover); # Manning’s n (-), increased with houses report hard = mask*0; # hard surface cells, not used here, included in house cover report roadwidth = road*celllength()*mask; # width of tarred roads in m report stone = mask*0; # stone fraction assumed zero everywhere

The infiltration maps in this example are for a 1 layer Green and Ampt infiltration. Additional surface types (crusted & compacted soils) each need two maps: one with the fraction in the gridcell and one with the ksat value. The ksat maps need not to be homogeneous, they can be derived from a classified map or from kriging for instance.