The Rasterscript


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] =…

where a large number of opertaions are available for the 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:

### input maps ###
mask =; # mask, 1 and missing value
DEM =; # digital elevation model, area must be <= mask
unitmapbase =; # consisting of basically two units: slope clayloam and wetland clay
road =; # road map, 0-1 for road fraction cover, contains all tarred roads
drains =; # primary and secondary drains
out =; # main outflowpoint, needed for a correct flow network
house_cover0 =; # fraction of house cover, from image
grasswid0 =; # grass strips in m
murrumroad =; # fraction murrum road in pixe, from sat img
veg_cover0 =; # fraction veg cover from sat img
barriers =; # anything that obstructs flooding: roads, very large houses
culverts =; # location with main culverts in primary drain
soiltbl = unitbase.tbl; # table with hydrological parameters:
# 01 ksat (mm/h)
# 02 porosity (cm3/cm3)
# 03 psi initial (cm)
# 04 initial moisture content (cm3/cm3)
# 05 RR (cm)
# 06 Manning’s n (-)

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.30 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.30 1.5 0.10

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.

### output maps ###
# basic topography related maps
Ldd =; # Local Drain Direction surface runoff
grad =; # slope, sine!
id =; # pluviograph influence zones (1-n)
outlet =; # location main outlet
landuse =; # land units (1-n)
roadwidth =; # road width, 0-cell size (m)
# vegetation maps
coverc=; # cover fraction (-)
lai=; # leaf area index (m2/m2) for interception storage
cropheight=; # plant height in m, for erosion, not used
grasswid =; # width of grass strips for infiltration
# Green and Ampt infiltration maps
ksat =; # sat hydraulic conductivity (mm/h)
pore =; # porosity (-)
thetai =; # initial moisture content (-)
psi =; # suction unsat zone (cm)
soildep =; # soil depth (mm), assumed constant here
# Surface maps
rr =; # surface roughness (cm)
mann =; # mannings n ()
stone =; # stone fraction on surface (-)
crust =; # crusted soil (-), not present
comp =; # compacted soil (-), murrum roads
hard =; # impermeable surfaces (0 or 1)
# erosion maps
cohsoil =; # cohesion (kPa)
cohplant =; # added root cohesion (kPa)
D50 =; # median of texture (mu)
aggrstab =; # aggregate stability number (-)
chancoh =; # channel cohesion (kPa)
# channel maps
lddchan =; # channel 1D network
chanwidth =; # channel width (m)
changrad =; # channel gradient, sine
chanman =; # channel manning (-)
chanksat =; # Ksat of the channel if it is not lined with concrete (mm/h)
chanside =; # angle channel side walls, 0 = rectangular, 1 = 45o
# channel flooding maps
chandepth =; # channel depth (m), if 0 the channel will not flood (infinitely deep)
chanmaxq =; # maximum discharge (m3/s) in culverts located in channel
# houses
housecov =; # house cover fraction (0-1)
roofstore =; # roof interception (mm)
raindrumsize = scalar(0); # raindrum size (m3)
drumstore =; # locations of rainwater harvesting in drums (0/1)

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 !!!
report grad = max(sin(atan(slope((DEM+barriers)*mask))), 0.005); # sine gradient (-), make sure slope > 0.005
report = 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.

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.

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.

report ksat = lookupscalar(soiltbl, 1, unitmap)*mask; # read column 1 of soil table
report pore = lookupscalar(soiltbl, 2, unitmap)*mask; # read column 2 of soil table
report psi = lookupscalar(soiltbl, 3, unitmap)*mask; # read column 3 of soil table
report thetai = lookupscalar(soiltbl, 4, unitmap)*mask; # read column 4 of soil table
report soildep = 2000.0 * mask; # soil depth 2 m everywhere
report = 2.78 * mask; # ksat of murrum roads, measured at 2.78 mm/h
report = 0 * mask; # ksat for crusted soils (mm/h)
report = 32.0 * mask; # ksat for grass strips (32mm/h)
report crust = mask*0; # crust fraction assumed zero everywhere
report comp = murrumroad*mask; # fraction compacted (murrum roads)