Introduction
The basin delineation system outlined in this blog only requires a Digital Elevation Model (DEM) as input. The quality of the DEM is, however, critical. This first part of the series of instructions deals with hydrological correction of DEMs using GRASS GIS.
Prerequisites
You must have setup GRASS 7 and imported a DEM as described under Installation & Setup in this blog. The DEM must have null (“no data”) for the water body into which the basins to delineate drain.
DEM errors affecting basin delineation
Basin delineation is grounded in the fact that water tends to flow downhill (from a higher to a lower potential, but strictly that does not necessarily mean downhill). This routing of water across the hillside and within water courses is the key for connecting any geographic point in the landscape to an outlet point.
DEM issues that can cause problems for flow routing include:
- regions lacking elevation data,
- artificial pits,
- artificial barriers,
- flat regions,
- wide river mouths.
Regions lacking data are generally difficult to mend. There are, however, special cases where it is necessary. In this tutorial you will mend small no data regions (a few cells) by first assigning a low elevation (typically 0) and then, optionally, apply a pit filling. It is crucial to remove land locked no data cells as these regions will otherwise gorge all the (virtual) water flow entering.
Artificial (virtual) pits are the most common problem, and most flow routing algorithms, including GRASS r.watershed can handle pits. For hydrological modeling it is often, but not always, better to mend artificial depressions. Also that will be done in the this manual.
Artificial barriers are more difficult to handle. And is not approached in this manual. You could, however, invert the DEM and solve problems with barriers as if they were pits.
Flat regions and wide river mouths cause similar flow routing problems. Different algorithms can be applied (i.e. Single flow Direction, or SFD, versus Multiple Flow Direction, MFD) for handling flat area. But a more ideal solution is to increase the vertical resolution and properly steer the routing using elevation data. This post includes a test of different flow directions algorithms, and in part 3 you will create an artificial slope directing outflow from wide river mouths.
Coastal DEM errors
Narrow bays, estuaries, lagoons etc. can cause problems in near shore regions. Both because they can come out as pits, but a worse problem is if such features become defined as land locked no data cells. That will cause problems both for identifying the basin outlets and later for modeling water flow out of the basin. To get rid of small “no data” regions along the coast you need to identify them and then assign an elevation that allows water to pass. The main objective of this post is to present a GRASS based scripting sequence that removes land locked no data holes in Digital Elevation Models.
GRASS preparations
There are several GRASS commands that can be used for mending and adjusting DEMs. The easiest is to use ordinary gap filling (r.fill.stats or r.fillnulls) but these methods do not consider the flow routing. Instead GRASS offers the r.fill.dir algorithm. r.fill.dir identifies pits along flow paths and then fills them. You can also use r.fill.dir for mending small no data holes, by first filling all the holes with a low elevation (typically DEM = 0).
The command r.flowfill is an alternative, available as an addon. To install GRASS addons you must have prepared the GRASS c-compiler as described in the post Install GDAL, QGIS and GRASS.
If you want to try r.flowfill , the command for installing it is:
> g.extension extension=r.fill.gaps
The rest of this manual however, uses r.fill.dir.
GRASS processing
The text sections below explains how to use GRASS for identifying and filling land locked no data regions. If you only want to run the process, all the commands are summarized further down.
Once you have filled the land locked no data holes you can use r.fill.dir fill them up with more realistic elevation values. Fur running the watershed analysis in parts 2 and 3, it is, however, enough to assign a low elevation (i.e. 0) to the no data holes. The GRASS process [r.watershed] manages pits on its own.
If you want to fill up pits in your DEM in a hydrologically sound manner then you should run r.fill.dir. For large grids that will, however, take a very (very) long time. For a DEM of 10 M pixels it can take 24 hours or longer. Thus you need to tile the DEM, with some overlap, before running r.fill.dir separately on each tile. To write such a command structure would also take a day. But with the help of the Python package basin_extract the GRASS commands will be setup in a second. Below both the principles for the pit filling, and how to use basin_extract, are covered. But first you need to remove the land locked no data regions.
Make the target directory
If you want to follow the standard of Karttur’s GeoImagine Framework, your target folder for the exported layers should be (where amazonia is the region used in this example):
> mkdir -p /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/script
Create a map of the terminal drainage
Create a map of the terminal drainage using the GRASS command r.mapcalc. For this to work, your original DEM (as described in the post on GRASS 7 setup ) must have the terminal drainage defined as “null” in the original DEM.
>r.mapcalc "drain_terminal = if(isnull('DEM@PERMANENT'), 1, null())"
Clump (r.clump) all contiguous cells together with a unique id to allow identifying large and small bodies representing the terminal drainage (null in the original DEM).
> r.clump input=drain_terminal output=terminal_clumps –overwrite
Convert the clumps to polygons with the command r.to.vect.
> r.to.vect input=terminal_clumps output=terminal_clumps type=area –overwrite
Add the area of each polygon to the the vector database with the command v.to.db.
> v.to.db map=terminal_clumps type=centroid option=area columns=area_km2 units=kilometers
Optionally export the vector with the polygons representing “no data” in the original DEM with v.out.ogr.
> v.out.ogr input=terminal_clumps type=area format=ESRI_Shapefile output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazoniax/0/terminal-clumps_SRTM_hydro-amazonia_0_cgiar-250.shp –overwrite
Fill in DEM over land
With the command v.to.rast it is possible to apply an SQL while converting vectors to a raster. Use that function to assign polygons below an area threshold to be the “no data” to fill, while polygons above the threshold should then represent the true terminal drainage. Set the pixel value to 0 - the preliminary elevation you will assign these cells. In the example I have set the threshold at 2.0 square kilometers. You should inspect the result and change the threshold if required.
> v.to.rast input=terminal_clumps type=area where="area_km2< 2.0" output=drain_terminal_small use=val value=0 –overwrite
> v.to.rast input=terminal_clumps type=area where="area_km2 >= 2.0" output=drain_terminal_large use=val value=0 –overwrite
if you can not identify a threshold that works you might have to consider manually editing the polygons. For example by deleting individual polygons in QGIS.
Superimpose the small “no data” clumps over the original DEM to create a DEM that is complete over land.
> r.mapcalc "inland_comp_DEM = if(( isnull(DEM@PERMANENT)), drain_terminal_small,DEM@PERMANENT )" –overwrite
Optionally export the completed DEM for land areas:
> r.out.gdal format=GTiff input=inland_comp_DEM output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazoniax/0/inland-comp-DEM_amazonia_0_cgiar-250.tif –overwrite
The DEM that you just created, inland_comp_DEM is sufficient for use with r.watershed in part 2. The remaining part of this post deals with filling up pits. But that is strictly not needed for running r.watershed and might even have detrimental consequences.
GRASS shell script
The shell script below contains all the GRASS commands explained above. You can run the commands below by copying them to a file with the extension .sh. Remember that you have to make the script file executable before you can run it (chmod u+x /path/to/file.sh or chmod 755 /path/to/file.sh). The run it from the Terminal window where GRASS is running.
# Create a map of the terminal drainage
mkdir -p /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/script
r.mapcalc "drain_terminal = if(isnull('DEM@PERMANENT'), 1, null())"
r.clump input=drain_terminal output=terminal_clumps --overwrite
r.to.vect input=terminal_clumps output=terminal_clumps type=area --overwrite
v.to.db map=terminal_clumps type=centroid option=area columns=area_km2 units=kilometers
# v.out.ogr input=terminal_clumps type=area format=ESRI_Shapefile output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazoniax/0/terminal-clumps_SRTM_hydro-amazonia_0_cgiar-250.shp --overwrite
# Fill in DEM over land
v.to.rast input=terminal_clumps type=area where="area_km2< 2.0" output=drain_terminal_small use=val value=0 --overwrite
v.to.rast input=terminal_clumps type=area where="area_km2 >= 2.0" output=drain_terminal_large use=val value=0 --overwrite
r.mapcalc "inland_comp_DEM = if(( isnull(DEM@PERMANENT)), drain_terminal_small,DEM@PERMANENT )" --overwrite
# r.out.gdal format=GTiff input=inland_comp_DEM output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazoniax/0/inland-comp-DEM_amazonia_0_cgiar-250.tif --overwrite
Restrict DEM filling to coastal strip
If you want to restrict the pit filling to the coastal strip you need to create a mask (r.mask). A mask limits the region for raster operations and will reduce the processing time.
To set a mask restricting the processing to the coastal strip you have to create a new layer for the terminal drainage - this time excluding the cells identified in drain_terminal_small.
> r.mapcalc "drain_terminal_v2 = if(isnull('inland_comp_DEM'), 1, null())" –overwrite
Then do an r.buffer starting from version 2 (v2) of the terminal drainage, set the maximum distance to the width of the coastal strip you want to analyse (e.g. 3 kilometers in the example).
> r.buffer input=drain_terminal_v2 output=coastal_strip distances=3000 units=meters –overwrite
Apply an r.mapcalc calculation to create a boolean layer of the coastal strip.
> r.mapcalc "coastal_mask = if((coastal_strip > 1), 1, null())"
To apply the mask by r.mask you just set it and it will remain in place until explicitly removed.
> r.mask raster=coastal_mask
Filling pits in the DEM
With all land areas having an assigned elevation, the next task is to fill up pits. If you want to fill up pits large than a single cell, you need to set a threshold for the maximum pit size to fill. This is done in one or two steps. First you have to decide whether to only fill up single cell pits (by setting the the parameter filldirsingle to 1 [= True]). Filling up single cell pits is more straight forward and is the recommended practice. If you set filldirsingle 0 [= False] you also need to set the maximum pit are to fill up (the parameter pitfillmaxarea), in the area units set to the vector database (square kilometers in the example below).
As noted above, you can not run r.fill.dir on large files. Here I only list the principal steps. To actually run the pit filling make use of the Python package basin_extract to create the required GRASS commands as explained further down.
The principal steps for filling up single cell pits (set the -f flag) in a DEM in GRASS:
- r.fill.dir -f originaldem -> filldem [the f flag forces filling of single cells only]
- r.mapcalc difference between filled and orignal dem
- r.to.vect for cells (points) with difference
- v.db.addcolumn for storing filled DEM value
- v.what.rast extract fill DEM to vector points
- v.patch to collect all vector points to fill in a single data source
- v.to.rast rasterize the DEM values
- r.mapcalc superimpose the filled DEM values over the original DEM
Overlaying the vector points on top of the original DEM can be done using different approaches. Apart from the GRASS comamnds used above (v.patch, v.to.rast and r.mapcalc) you can also make use of GDAL. This is discussed further down.
The principal steps for filling up also pits larger than a single cell in a DEM in GRASS:
- r.fill.dir originaldem -> filldem2 [all pit cells filled regardless of area]
- r.mapcalc difference between filled and orignal dem
- r.to.vect for regions (areas) with difference
- v.db.addcolumn for storing filled DEM value
- v.what.rast extract fill DEM to polygon centroids
-
- v.patch to collect all vector areas to fill in a single data source
- v.to.rast rasterize the DEM values
- r.mapcalc superimpose the filled DEM values over the original DEM
Overlaying the vector points on top of the original DEM can be done using different approaches. Apart from the GRASS comamnds used above (v.patch, v.to.rast and r.mapcalc) you can also make use of GDAL. This is discussed further down.
The process steps above can not be run for very large grids - it takes too long time. But pits are usually small. Larger pits (hundreds of square kilometers or larger) usually represent endorheic basins and should not be filled. Thus you must, and can, partition large DEMs into smaller tiles and then run each tile separately. With an overlap that allows you to capture pits up to the size you deem relevant.
Create tiled process loop
The example DEM used in this tutorial of central South America is approximately 22700 × 22500 cells. This is far too large to process in one go. The Python package developed for assisting the basin delineation, basin_extract can generate a tiled process loop for GRASS.
Parameterisation
The Python script basin_extractor is parameterised from an XML file that follows the general structure of Karttur´s GeoImagine Framework. The structure below shows the setup for generating the GRASS shell script for a tiled processing of pit filling of the DEM over central South America (“Amazonia”).
<basindelineate>
<userproj userid = 'karttur' projectid = 'karttur' tractid= 'amazonia' siteid = '*' plotid = '*' system = 'MODIS'></userproj>
<!-- Define processing -->
<process processid ='BasinDelineate'>
<overwrite>Y</overwrite> # Y(es) or N(o)
<delete>N</delete> # Y(es) or N(o)
<parameters
stage = '0' # 0, 1, 2 or 3
adjacentdist = '330' # any number in map distance units
outlet = 'SFD' # MOUTH, SFD or MFD
distill = '' # MFD or None (nothing)
clusteroutlet = 'central' # central or max
basinthreshold='2000' # any number is cells
watershed='MFD5' # SFD, MFD or MFD? (? in range 1..10)
verbose = '1' # 0, 1 or 2
proj4CRS = '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs'
>
</parameters>
<srcpath volume = "GRASS2020"></srcpath>
<dstpath volume = "GRASS2020"></dstpath>
<srccomp>
<!-- the name <tag> and prefix must be the same and can not be altered. The script looks for which source data is included and compares that with the parameters, and then sets the processing method if these corresponds. Otherwise the script ends with an error message.
-->
<basin-mouth-outlet-pt source = "SRTM" product = "drainage" folder = "basin" band = "basin-mouth-outlet-pt" prefix = "basin-mouth-outlet-pt" suffix = "cgiar-250">
</basin-mouth-outlet-pt>
<basin-MFD-outlet-pt source = "SRTM" product = "drainage" folder = "basin" band = "basin-MFD-outlet-pt" prefix = "basin-MFD-outlet-pt" suffix = "cgiar-250">
</basin-MFD-outlet-pt>
<basin-SFD-outlet-pt source = "SRTM" product = "drainage" folder = "basin" band = "basin-SFD-outlet-pt" prefix = "basin-SFD-outlet-pt" suffix = "cgiar-250">
</basin-SFD-outlet-pt>
<shorewall-pt source = "SRTM" product = "drainage" folder = "basin" band = "shorewall-pt" prefix = "shorewall-pt" suffix = "cgiar-250">
</shorewall-pt>
</srccomp>
<dstcomp>
<basin-outlet>
</basin-outlet>
<!--
To force a non-default destination path and name, you can enter all components explicitly, e.g.
<basin-outlet band = "basin-mfd+mfd-max" prefix = "basin-mfd+mfd-max">
</basin-outlet>
-->
</dstcomp>
</process>
</basindelineate>
The core parameters of relevance for this stage (0) include:
- stage [0, 1, 2 or 3] (this part is for stage = ‘0’)
- adjacentdist [not relevant for stage 0]
- outlet [not relevant for stage 0]
- distill [not relevant for stage 0]
- clusteroutlet [not relevant for stage 0]
- basinthreshold [not relevant for stage 0]
- watershed [not relevant for stage 0]
- pitfill
- tilesize
- tileoverlap
- single [Boolean]
- verbose [0, 1 or 2]
- proj4CRS [‘proj4CRS definition for output data if input layers lack proper projection information’]
Setup xml file and run basin_extract
Use the xml template file above to setup the processing for your own xml. Set the parameter stage to 0 and run basin_extract pointing to the xml.
Running basin_extract with stage = 0 generates two shell script files, one for GRASS and one for OGR2OGR. You do not need to use the OGR2OGR file, it is just supplied in case you prefer that.
For each tile that is processed the GRASS shell script includes the following set of commands (excluding comments):
# This example shows the script for single cell filling (with the -f flag set):
# region is automatically set from the tiling process
g.region -a n=-3660170.456407 s=-3914992.450317 e=-8826107.250690 w=-9080929.244600
# the -f flag restricts the filling to single, isolated cells
r.fill.dir -f input=inland_comp_DEM output=hydro_cellfill_DEM_0_0 direction=hydro_cellfill_draindir_0_0 areas=hydro_cellfill_problems_0_0 --overwrite
# get the difference between the filled and the original DEM
r.mapcalc "DEM_cellfill_diff_0_0 = hydro_cellfill_DEM_0_0 - inland_comp_DEM" --overwrite
# create boolean mask of changed cells
r.mapcalc "inland_fill_cell_0_0 = if(DEM_cellfill_diff_0_0 != 0, 1, null() )" --overwrite
# convert the cells with altered elevation to a point vector
r.to.vect input=inland_fill_cell_0_0 output=inland_fill_pt_0_0 type=point --overwrite
# add column to vector file with altered elevations
v.db.addcolumn map=inland_fill_pt_0_0 columns="filldem DOUBLE PRECISION"
# read the altered DEM values to the point vector
v.what.rast map=inland_fill_pt_0_0 column=filldem raster=hydro_cellfill_DEM_0_0
# export the point vector to an ESRI shape file
v.out.ogr input=inland_fill_pt_0_0 type=point format=ESRI_Shapefile output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-pt-0-0.shp --overwrite
Superimpose filled DEM over original DEM
All areas identified as pits are captured in the sequence of tiles. You now have several alternatives for superimposing these values over the original DEM. All alternatives can be applied using a tile by tile approach, or applied for the whole region after patching the tiles.
Filling of pits affects flow routing and hydrological modeling. When the flow routing enters a flat area, an algorithm traversing the flat area towards the lowest drainage(s) out of the flat must be identified. Dependent on the parameterization of the flow dividsion (single or multiples and the spread if multiple) the flat area will be wetted differently. This will later also affect the estimation of soil moisture, evapotranspiration and flooding tendencies of the cells in the flat area (in general all tend to increase). The most widely spread recommendation is to allow filling of single cell pits, but be more careful with filling larger and larger pits.
On the other hand, DEMs generated by SAR data have a tendency to underestimate the elevation of water surface (due to radar signal double bouncing on water surface and shoreline structures). Thus no general recommendation can be given regarding how to handle larger depressions in DEMs.
You have different alternatives for superimposing The alternatives for filling pits that I have tested include:
-
- superimpose tile by tile with GDAL
- patch up vectors in GDAL, superimpose
- patch up vectors in GRASS, rasterize, superimpose
GDAL is quicker and can be used for burning vector data on a raster. The burning can be done on a tile by tile basis, or after patching up the tiles to a single vector data source. With GRASS you have to first patch the vector tiles, then rasterize them, and then superimpose. This takes longer and demands more processing capacity.
Patching tiles (optional)
If you use GDAL (GDAL_rasterize) for burning the filled DEM on the original DEM you do not need to patch the tiles. Just continue to the section on Rasterize using GDAL below. The shell script for the burning is already produced by the Python package basin_extract.
You can use two different methods for patching together the tiles. The first alternative is to use the GRASS command v.patch. The commands for patching all the times are given at the end of the GRASS script file:
g.region raster=DEM@PERMANENT
MAPS=$(g.list type=vector separator=comma pat="inland_fill_pt_*")
v.patch -e input=$MAPS output=inland_fill_pt
MAPS=$(g.list type=vector separator=comma pat="inland_fill_area_*")
v.patch -e input=$MAPS output=inland_fill_area
The other alternative is to use the GDAL vector command ogr2ogr. The command line process for patching together all vectors are prepared in the script file “region”_ogr_part0.sh, produced and reported by the Python package basin_extractor. The first lines of the shell script looks like this:
ogr2ogr -skipfailures -nlt /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/flowdir-pt-dem_drainage_amazonia_0_cgiar-250.shp /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-pt-0-0.shp
ogr2ogr -skipfailures -nlt /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/flowdir-pt-dem_drainage_amazonia_0_cgiar-250.shp /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-area-0-0.shp
ogr2ogr -append -skipfailures -nlt /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/flowdir-pt-dem_drainage_amazonia_0_cgiar-250.shp /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-pt-0-1.shp
ogr2ogr -append -skipfailures -nlt /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/flowdir-pt-dem_drainage_amazonia_0_cgiar-250.shp /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-area-0-1.shp
The two first lines coopies the input vectors to a new file, all following lines appends the input vectors to the initial output file.
Patch up vectors in GRASS, rasterize, superimpose
It is more complicated to superimpose the filled elevation values using GRASS compare to GDAL. For GRASS you first have to patch up the vectors, then rasterize and finally superimpose. The required sequence of commands is included in the GRASS shell script file that was produced when running basin_extract for stage0. The very last lines of the that shell script file should be something like:
v.to.rast input=inland_fill_area where=”area_km2<2.0” output=inland_fill_accept use=attr attribute_column=filldem –overwrite</span>
r.mapcalc inland_fill_dem = if( isnull(inland_fill_accept), inland_comp_DEM, inland_fill_accept)
Rasterize using GDAL
With GDAL you can instead use GDAL_rasterize to burn the DEM values from the vectors with the filled DEM value directly over the original DEM. The burning can be done using either the vector tiles or a patched vector. By using vector tiles you bypass the need for patching tiles.
That requires fewer step and is already prepared for running in the two shell script files:
- amazonia_GDAL_pt_part0.sh
- amazonia_GDAL_area_part0.sh
Here are the first lines from the latter:
GDAL_rasterize -a filldem /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-area-0-0.shp /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/inland-comp-demm_drainage_amazonia_0_cgiar-250.shp
GDAL_rasterize -a filldem /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/stage0/inland-fill-area-0-1.shp /Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazonia/0/inland-comp-demm_drainage_amazonia_0_cgiar-250.shp
Visualizing the DEM errors and fixes
If you want to try the hydrological fixing of DEM using r.fill.dir, you need to start by setting a subregion to process with g.region. About 1 million cells is a suitable region.
> g.region -a n=”float” s=”float” e=”float” w=”float”
If you do not set a subregion, the process will take several hours, or even days. With a small subregion set , you can run r.fill.dir.
for single cell pits only:
> r.fill.dir -f input=inland_comp_DEM output=hydro_cellfill_DEM_ direction=hydro_cellfill_draindir areas=hydro_cellfill_problems –overwrite
or for all depressions:
> r.fill.dir input=inland_comp_DEM output=hydro_areafill_DEM_ direction=hydro_areafill_draindir areas=hydro_areafill_problems –overwrite
Calculate the difference between the pit filled and the original DEM using r.mapcalc.
> r.mapcalc “DEM_fill_diff = hydro_areafill_DEM - inland_comp_DEM” –overwrite
Set a color ramp to the difference layer.
> r.colors DEM_fill_diff color=differences
Export the difference layer to a goetiff.
> r.out.gdal format=GTiff input=DEM_fill_diff output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazoniax/0/DEM-fill-diff_amazonia_0_cgiar-250.tif –overwrite
Calculate a boolean map of the difference area.
> r.mapcalc “inland_fill_areas = if(DEM_fill_diff != 0, 1, null() )” –overwrite
Convert the boolean map of identified pits to a vector.
> r.to.vect input=inland_fill_areas output=inland_fill_areas type=area –overwrite
Export the vector with identified pits to an ESRI shape file.
> v.out.ogr input=inland_fill_areas type=area format=ESRI_Shapefile output=/Volumes/GRASS2020/GRASSproject/SRTM/region/basin/amazoniax/0/inland-fill-areas_SRTM_hydro-amazonia_0_cgiar-250.shp –overwrite
You can inspect the exported files in for instance QGIS
You can also open a monitor for GRASS directly from the command line and then add the layers from the command line.
> d.mon wx0
> d.shade shade=inland_comp_DEM
> color=inland_comp_DEM
> d.vect inland_fill_areas type=boundary color=red
