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:
- LISEM is largely built in a mixture of PCRaster commands and C++
- 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:
| ### input maps ### |
|
| mask = mask20m.map; | # mask, 1 and missing value |
| DEM = dem20m.map; | # digital elevation model, area must be <= mask |
| unitmapbase = lu20mbase.map; | # consisting of basically two units: slope clayloam and wetland clay |
| road = roads20m.map; | # road map, 0-1 for road fraction cover, contains all tarred roads |
| drains = chanmask20m.map; | # primary and secondary drains |
| out = mainoutlet20m.map; | # main outflowpoint, needed for a correct flow network |
| house_cover0 = house_cover20m.map; | # fraction of house cover, from image |
| grasswid0 = grasswid0.map; | # grass strips in m |
| murrumroad = murrum20m.map; | # fraction murrum road in pixe, from sat img |
| veg_cover0 = veg_cover20m.map; | # fraction veg cover from sat img |
| barriers = barriers20m.map; | # anything that obstructs flooding: roads, very large houses |
| culverts = culverts20m.map; | # 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 = ldd.map; | # | Local Drain Direction surface runoff |
| grad = grad.map; | # | slope, sine! |
| id = id.map; | # | pluviograph influence zones (1-n) |
| outlet = outlet.map; | # | location main outlet |
| landuse = landunit.map; | # | land units (1-n) |
| roadwidth = roadwidt.map; | # | road width, 0-cell size (m) |
| # vegetation maps | ||
| coverc= per.map; | # | cover fraction (-) |
| lai= lai.map; | # | leaf area index (m2/m2) for interception storage |
| cropheight= ch.map; | # | plant height in m, for erosion, not used |
| grasswid = grasswid.map; | # | width of grass strips for infiltration |
| # Green and Ampt infiltration maps | ||
| ksat = ksat1.map; | # | sat hydraulic conductivity (mm/h) |
| pore = thetas1.map; | # | porosity (-) |
| thetai = thetai1.map; | # | initial moisture content (-) |
| psi = psi1.map; | # | suction unsat zone (cm) |
| soildep = soildep1.map; | # | soil depth (mm), assumed constant here |
| # Surface maps | ||
| rr = rr.map; | # | surface roughness (cm) |
| mann = n.map; | # | mannings n () |
| stone = stonefrc.map; | # | stone fraction on surface (-) |
| crust = crustfrc.map; | # | crusted soil (-), not present |
| comp = compfrc.map; | # | compacted soil (-), murrum roads |
| hard = hardsurf.map; | # | impermeable surfaces (0 or 1) |
| # erosion maps | ||
| cohsoil = coh.map; | # | cohesion (kPa) |
| cohplant = cohadd.map; | # | added root cohesion (kPa) |
| D50 = d50.map; | # | median of texture (mu) |
| aggrstab = aggrstab.map; | # | aggregate stability number (-) |
| chancoh = chancoh.map; | # | channel cohesion (kPa) |
| # channel maps | ||
| lddchan = lddchan.map; | # | channel 1D network |
| chanwidth = chanwidt.map; | # | channel width (m) |
| changrad = changrad.map; | # | channel gradient, sine |
| chanman = chanman.map; | # | channel manning (-) |
| chanksat = chanksat.map; | # | Ksat of the channel if it is not lined with concrete (mm/h) |
| chanside = chanside.map; | # | angle channel side walls, 0 = rectangular, 1 = 45o |
| # channel flooding maps |
||
| chandepth = chandepth.map; | # | channel depth (m), if 0 the channel will not flood (infinitely deep) |
| chanmaxq = chanmaxq.map; | # | maximum discharge (m3/s) in culverts located in channel |
| # houses | ||
| housecov = housecover.map; | # | house cover fraction (0-1) |
| roofstore = roofstore.map; | # | roof interception (mm) |
| raindrumsize = scalar(0); | # | raindrum size (m3) |
| drumstore = drumstore.map; | # | 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 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.
| ### INFILTRATION MAPS for 1 layer GREEN & AMPT### | |
| 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 ksatcomp.map = 2.78 * mask; | # ksat of murrum roads, measured at 2.78 mm/h |
| report ksatcrust.map = 0 * mask; | # ksat for crusted soils (mm/h) |
| report ksatgras.map = 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) |
THIS PAGE IS UNDER CONSTRUCTION !!!!