GRASS SEEDS A BEGINNER'S TUTORIAL Developed by Mitchel Langford Joseph Wood With the assistance of Caroline Crumpton Alistair Duncan Alan Strachan David Unwin ASSIST - Academic Support For Spatial Information Systems Copyright @ 1995 Midlands Regional Research Laboratory Step One: In and out =================================================================== A: Introducing GRASS Introduction This Training Session introduces you to the GRASS GIS and allows you to gain some familarity with its user interface. Many of the commands that are most commonly used for managing and displaying data within GRASS will be covered here. The concept of seperate map layers, each relating to a particular theme is also introduced. Starting GRASS Before you can carry out the exercises contained within this volume, you must first start GRASS and tell it which dataset you are going to work with. When you start GRASS (probably by typing grass4.1), you should be presented with the following screen. If for any reason this screen does not appear, contact your system manager, or refer to installation documentation to check that everything has been set up correctly. grass@congo[1]> grass4.1 Display 1 The usual GRASS opening screen. Note that the Location, Mapset and Database may vary from machine to machine. ----------------------------(clear)-------------------------------- GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: tutor_________ (enter list for a list of locations) MAPSET: PERMANENT_____ (or mapsets within a location) DATABASE: /usr/gis__________________________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) You new need to tell GRASS which data set you will be working with. The opening screen prompts you to identify the LOCATION, MAPSET, and DATABASE. LOCATION identifies the geographical extent of the study region. For these exercises, you will be working with a small part of Leicestershire, England. Press the return key a few times until the cursor is positioned on the LOCATION line. Then type tutor to identify the study region. Don't worry about overtyping the default setting - just make sure that any remaining letters are deleted by presing the space bar. Press return once more to move the cursor onto the MAPSET line. A mapset in GRASS is your own set of files that reside in the current LOCATION. For convenience, you should give your "user name". This will ensure that any files you alter or create will note be inadvertantly changed by anyone else. (Note that the default setting of PERMANENT contains the set of files that can be read by anyone using the system). The DATABASE should not need changing on this occasion. It is simply the directory in which the sample locations and mapsets are stored. A note for system managers - the DATABASE directory can be found by typing setenv (or equivalent shell command) before starting GRASS. The environment variable GISDBASE should indicate the directory holding the data. When you have set the correct options, press to start. ----------------------------(clear)-------------------------------- Welcome to GRASS 4.1 (Spring 1993) Update package 2 Geographic Resources Analysis Support System (GRASS) is a Trademark of U.S. Army Construction Engineering Research Laboratories (USACERL) New releases of GRASS are coordinated and produced by the Office of GRASS Integration (OGI) located at USACERL. This version running thru the C Shell (/bin/csh) Help is available with the command: g.help When ready to quit enter: exit Mapset in Location GRASS 4.1 > exit ----------------------------(clear)-------------------------------- GRASS SESSION WRAPUP You have just finished working on mapset: The following RASTER maps belong to it: contours image plant rail segment spillage urban crash landcov popln roads source topo water There are no VECTOR maps in this mapset There are no SITES maps in this mapset Do you wish to selectively remove data files? y/n [n] n ----------------------------(clear)-------------------------------- GOOD BYE from GRASS grass@congo[2]> pwd /usr/gis grass@congo[3]> ls ./ .show_pro conserv.ung* sgi/ ../ .workspace/ dumpster/ src/ .cshrc* Install* exactowarp/ src.alpha/ .grassrc Mosaic-sgi-2.7b5* jim.prm src.tar.Z .login* Upgrade* man/ system/ .notepad WorkSpace/ niwot-grey.tiff tutor/ .profile* arc/ rastest.asci xv-3.10a/ .seahavensave conserv.dxf robb/ grass@congo[4]> grass@congo[51]> cd tutor grass@congo[5]> ls ./ PERMANENT/ ../ ch_grass_notes.txt grass@congo[6]> pwd /usr/gis/tutor grass@congo[7]> Step Two: Intro to Grass Module A ================================================================================ (repeat above startup until grass prompt as follows) Mapset in Location GRASS 4.1 > ASSIST GRASS Seeds - A Beginner's Tutorial -------------------------------------------------------------------------------- A: Introducing GRASS Aims: Familiarisation with the GRASS GIS. Objectives: On completion of this Training Session you will be able to :- Obtain a list of the current map layers Investigate the contents of any map layer Display the map layers on the computer screen Construct a compostie thematic map Update map layer support files Display 2 GRASS commands to be used g.list d.mon d.rast r.report d.legend d.erase d.colormode r.colors d.colortable d.histogram r.patch r.support g.remove exit Step 1: Listing the available map layers. You first need to obtain a list of all the information that is currently held on the system. Type the command g.list and press the Return/Enter key. GRASS 4.1 > g.list You will then be presented with 8 possible file types to be listed. All the map layers in this exercise are raster files, so select 1 and press the return key. ----------------------------(clear)-------------------------------- LIST FACILITY This program allows you to list files from mapsets in your search path Please select the type of file to be listed 1 raster files 2 binary vector files 3 paint icon files 4 paint label files 5 site list files 6 region definition files 7 imagery group files 8 3D view parameters RETURN to exit You will then be presented with 8 possible file types to be listed. All the map layers in this exercise are raster files, so select 1 and press the return key. > 1 ---------------------------------------------- raster files available in mapset PERMANENT: contours image plant rail segment spillage urban crash landcov popln roads source topo water ---------------------------------------------- hit RETURN to continue --> Discussion: You will notice that the information is held in GRASS as a number of independent 'map layers'. Each map layer holds data for a single 'theme' such as roads, urban areas, or population statistics. This feature of seperate data 'layers' is very common in Geographical Information Systems. You will often need to use the command g.list to remind yourself of the map layers that are available, especially when you start to generate new layers through GIS operations. Check that all the map layers listed in Display 3 appear on your computer screen in mapset PERMANENT; if not please alert your system manager or course tutor. Press return once more when you have finished. Display 3 Rastor map layers available for this tutorial. File name Image title crash Site of motorway crash contours Topographic contour data image Satellite image landcov Land cover / land use plant Sewage treatment plant popln Population data rail Railways roads Road network segment Polluted river segment source Sewage source spillage Chemical spillage topo Degital elevation model urban Major urban areas water Rivers and reservoirs > Step 2: Viewing the study area. To display any graphics in GRASS, you first need to open a 'graphics window'. GRASS 4.1 > d.mon start=x0 Graphics driver [x0] started Discussion: In GRASS you can open up to 7 similargraphics windows simultaneously. These are labeled x0 to x6. Each window can be moved and resized to fill up any portion of the screen. If you wishto move the graphics window, place the pointer over the border of the graphics window, hold down the left mouse button, and move to the desired position. To resize the window, place the pointer over on corner of the window and move the mouse with the left button held down. When you are happy with the window position, you need to tell GRASS to select that window for output. GRASS 4.1 > d.mon select=x0 To help you become a little more familiar with the study area, you will first display a satellite image which provides a synoptic view of the geography of the region. GRASS 4.1 > d.rast image Discussion: This image is not unlike a normal black-and-white photograph taken from space, although it actually records infra-red light intensity. It is intended to provide you with a "real world" view of the area that we will be studying. However, what we say about the "real world" is constrained by the way space is modelled in GRASS; remember this is a rastor GIS, ie one that uses a lattice of grid cells to represent space. Compare this image with the sketch map of the area (see Map 1 on the following page). What land use features can you identify? User Notes:................................................................ Map 1 : The Study Area (missing) Step 3: Describing the contents of a rastor map layer. In this step you will examine the parameters of a rastor data file, in terms of its physical size and the dimensions, or 'resolution', of the cell lattice. GRASS 4.1 > r.report roads units=kilom,percent r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Tue Jul 15 19:00:33 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: Road Network (roads in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | square | % | | #|description |kilometers| cover| |-----------------------------------------------------------------------------| | 0|Background . . . . . . . . . . . . . . . . . . . . . . .|131.817500| 91.54| | 3|Motorway . . . . . . . . . . . . . . . . . . . . . . . .| 0.735000| 0.51| | 6|Proposed . . . . . . . . . . . . . . . . . . . . . . . .| 0.102500| 0.07| | 9|'A' roads. . . . . . . . . . . . . . . . . . . . . . . .| 2.107500| 1.46| |11|'B' roads. . . . . . . . . . . . . . . . . . . . . . . .| 3.447500| 2.39| |12|'C' roads. . . . . . . . . . . . . . . . . . . . . . . .| 5.790000| 4.02| |-----------------------------------------------------------------------------| |TOTAL |144.000000|100.00| +-----------------------------------------------------------------------------+ Discussion: The system will respond with a table of information. Under the banner containing the location and date, the boundaries and resolution of the region are indicated. In this case units are in national grid coordinates and are therefore in meters. From this it should be possible to calculate the number of rows and columns in the image. Below the region information, the title of the rastor layer is given along with the name of each category in that layer. The extent of each catagory is given in the units requested (km2 and % in this case). For each of the map layers in Display 4 below, use r.report to find the most and least frequent category. For which rastor map layers is this a useful process, and for which is it not? Why? Display 4 Rastor map Most Frequent Least Frequent layer catagory catagory landcov ............. ............. popln ............. ............. rail ............. ............. roads ............. ............. topo ............. ............. Discussion: You should have found that the report for the topo map layer scrolled through your window too fast to read. GRASS does not take into account the number of categories in a map layer when displaying textual reports. For a more 'user friendly' display of the data, you can combine the GRASS command r.report with the UNIX command more. more is a text display utility that pauses between successive screens when displaying larger amounts of text. To move on a screen press the space bar, or to quit press q. GRASS 4.1 > r.report topo units=kilom,percent | more Discussion: The ability to mix GRASS and UNIX commands in this way is a feature of GRASS that gives it far greater power and flexibility. This is especially true when interfacing GRASS with other UNIX programs. Step 4: Veiwing the map layers. You will now use the display mechanism, first seen in Step 2, to view all the GIS map layers that are listed in Display 3. GRASS 4.1 > d.rast roads Discussion: The command d.rast is the principal means of displaying raster map layers in GRASS. If you refer back to the map of the study area (Map 1) it will become apparent that the road network seen on the screen is a simplified version of that seen on the map. The degree of detail and accuracy of GIS map layers varies depending on the data source and the intended application, but it is always the case that they are simplifications or models of the real world. To make the map more useful, it is possible to display the legend associated with any map layer. As a first step, display the legend in a new graphics window. GRASS 4.1 > d.mon start=x1 select=x1 Error: no such monitor 'x1' No such monitor as 'x1' Problem selecting x1. Will try once more No such monitor as 'x1' GRASS 4.1 > d.legend roads If necessary, reposition the new window with the mouse. Now try displaying another raster layer with its associated legend. GRASS 4.1 > d.mon select=x0 GRASS 4.1 > d.rast water GRASS 4.1 > d.mon select=x1 No such monitor as 'x1' GRASS 4.1 > d.erase GRASS 4.1 > d.legend water lines=20 Discussion: The first step tells GRASS to direct output to the original graphics window x0. The rastor layer water is then displayed. To display the legend, the second 'legend window' x1 is used. Before drawing a second legend in this window it is necessary to erase the old one. When displaying legends GRASS will, by default, make each legend box and key as large as possible to fit the available space. Since the raster map layer water has only one catagory (labeled as four), a very large legend is produced. By using the option lines=20, GRASS scales the legend as if 20 catagories are to be displayed. Repeat the proceedure for rail, urban, and topo, adjusting the lines option if you feel it necessary. For what type of rastor map layers is d.legend most appropriate? Display the satellite image again. GRASS 4.1 > d.mon select=x0 GRASS 4.1 > d.rast image Discussion: The image should consist of approximately 8 shades of grey - enough to identify most of the features of the study area. As GRASS is capable of displaying up to 240 colors on a SUN workstation, why are there only 8 grey shades making up the image? GRASS has a fixed palette of colors all of which can be displayed simultaneously. Since an image can consist of any range of colors, only a few can be reserved for grey scales. GRASS simply matches each raster cell value to the nearest available color. We can produce an improved image by abandoning GRASS' fixed pallette and dedicating most of the range of colors to one image. Make sure the pointer is resting over the graphics window x0. GRASS 4.1 > d.colormode float GRASS 4.1 > d.rast image Discussion: You should see an improved satellite image of the study area. The majority of the SUN's 256 colors have now been allocated to grey values in this one window. When using this floating color table, GRASS allocates the correct colors to the window on which the pointer is resting. Try moving the pointer between different graphics windows. What happens to the colors in the satellite and legend windows when you do this? When is it most appropriate to use a floating color table and when is a fixed colortable more appropriate? User Notes:................................................................ Finally, try changing the color scheme of the satellite image. GRASS 4.1 > r.colors image color=wave Color table for [image] set to wave GRASS 4.1 > d.rast image Discussion: Every raster file in GRASS is associated with a color file that stores the color associated with each raster category. The color file can be altered in a number of ways. One of the simplest is to use one of GRASS' preset palettes. In the example above, the palette has been changed from a grey scale to a red-green-blue color wave. Display the complete palette associated with a particular file. GRASS 4.1 > d.erase GRASS 4.1 > d.colortable image Display the range of possible preset color palettes that can be called using r.colors. GRASS 4.1 > r.colors help Usage: r.colors [-wq] map=name color=type Flags: -w Don't overwrite existing color table -q Quietly Parameters: map raster map name color type of color table options: aspect,grey,grey.eq,gyr,rainbow,ramp,random,ryg,wave,rules Where color type is one of: aspect (aspect oriented grey colors) grey (linear grey scale) grey.eq (histogram equalized grey scale) gyr (green through yellow to red colors) rainbow (rainbow color table) ramp (color ramp) ryg (red through yellow to green colors) random (random color table) wave (color wave) rules (create new color table by rules) Discussion: All the options available within r.colors should now be presented. The available options for any command can be found by typing the command followed by the word help. Alternatively, if you just type the command by itself, you will be prompted for each of the available options. Appendix A contains a short summary of the uses and all the options of all GRASS commands used in this tutorial. Change the image color palette back to some, or all, of the preset palettes found from r.colors help. Can you think of uses for each of the palettes? Finally, change the color palette back to the original grey scale and reset the colormode back to the fixed palette. GRASS 4.1 > r.colors image color=grey Color table for [image] set to grey GRASS 4.1 > d.colormode fixed GRASS 4.1 > d.erase User Notes:................................................................ Step 5: Displaying simple statistics of a map layer. You will now use GRASS to display some statistical information concerning the map layers. GRASS 4.1 > d.histogram roads r.stats: 100% Discussion: This provides a graphical summary of the cell values held within a map layer. A range of potentially useful information is presented. This includes the number of values within a given class interval and the proportion of total cells in each. This form of output also provides a clear representation of the relationship between cell value and image color. If you wish to see the labels associated with each catagory, use r.report. GRASS 4.1 > r.report roads r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Tue Jul 15 20:47:26 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: Road Network (roads in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | | #|description | |-----------------------------------------------------------------------------| | 0|Background | | 3|Motorway | | 6|Proposed | | 9|'A' roads | |11|'B' roads | |12|'C' roads | +-----------------------------------------------------------------------------+ Find how many cells in the road map layer are coded as 'C' class roads bu using a d.histogram and r.report. What proportion of the total number of cells is this? If you can't remember how to use r.report to do this, refer back to step 3. Now have a go at generating image histograms of any of the map layers that are available to you. Refer to Appendix A or use d.histogram help to find out how to display pie-charts instead of histograms - experiment a little! User Notes:................................................................ Step 6: Combining the map layers to give a composite 'map'. You will now attempt to preform a fundimental GIS operation. This is to take two map layers registered to the same area of the earth's surface and combine them to produce a single map layer. The map layers to be combined are: urban water rail roads Discussion: Overlay can be thought of in two different ways. Firstly, it is possible to preform a 'geographical overlay' in which several raster map layers are displayed together in a graphics window. This process is purely for visualisation as no permanent record of the overlay is created. Secondly, raster layers can be combined to create new raster layers. This is often a more important process, as the new layers can be further manipulated within the GIS. First, create a graphic overlay. GRASS 4.1 > d.rast urban GRASS 4.1 > d.rast -o water GRASS 4.1 > d.rast -o rail GRASS 4.1 > d.rast -o roads Discussion: By using the -o option when calling d.rast, the previous raster layer is not erased as the new one is plotted. Since most of each raster layer is made up of 'background' (category zero), each successive overlay does not completely cover the previous one. Now try the alternative approach and create a new 'composite' map layer. GRASS 4.1 > r.patch input=roads,rail,water,urban output=composite r.patch: percent complete: 100% CREATING SUPPORT FILES FOR composite Discussion: r.patch allows different map layers to be combined whilst preserving the colors and labels associated with each category. Only empty (category zero) cells can be replaced with the contents of a new layer. Display the new map layer composite (refer back to Step 4 if you have forgotten how to do this). What is the significance of the order in which the raster map layers are overlain? Can you account for the differences in order for the graphical and file based overlays? If you are unsure, try changing the order and see! Step 7: Putting the final touches on the composite map. The final map layer generated in Step 6, and held in the file composite, is produced by combining a number of map layers. Unfortunately, this new map still has the title of the first input layer to r.patch (roads). As the final step in this Training Session you will use the GRASS command r.support to change the title. By issuing this command you should get some idea of the support files associated with a raster map layer in GRASS. GRASS 4.1 > r.support Enter name of raster file for which you will create/modify support files Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > composite Edit the header for [composite]? (y/n) [n] y Edit header for [composite] cellhd compression: 1 3.0 compression indicated pre 3.0 compression not indicated hit RETURN to continue --> ----------------------------(clear)-------------------------------- Please enter the following information for [composite] 240 Number of rows 240__ Number of cols 1____ Number of bytes per cell AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- IDENTIFY CELL HEADER ============================= DEFAULT REGION ======== | Default North:322000 | | | | ======= CELL HEADER ======= | | | NORTH EDGE:322000_____ | | | | | | Def. West |WEST EDGE | |EAST EDGE | Def. East 444000 |444000_____| |456000_____| 456000 | | SOUTH EDGE:310000_____ | | | ============================= | | | | Default South:310000 | ===================================================== PROJECTION: 1 (UTM) ZONE: 0 AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- projection: 1 (UTM) zone: 0 north: 322000 south: 310000 east: 456000 west: 444000 e-w res: 50 n-s res: 50 total rows: 240 total cols: 240 total cells: 57,600 Do you accept this header? (y/n) [y] > y header for [composite] updated hit RETURN to continue --> Update the stats (histogram,range) for [composite]? (y/n) [n] n Edit the category file for [composite]? (y/n) [n] y ----------------------------(clear)-------------------------------- The current value for the highest category in [composite] is shown below If you need to change it, enter another value Highest Category: 15___ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- [composite] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: Road Network________________________________________________ CAT NEW CATEGORY NAME NUM 0 Background__________________________________________________ 1 ____________________________________________________________ 2 ____________________________________________________________ 3 Motorway____________________________________________________ 4 Water_______________________________________________________ 5 Railway_____________________________________________________ 6 Proposed____________________________________________________ 7 ____________________________________________________________ 8 ____________________________________________________________ 9 'A' roads___________________________________________________ Next category: 10___ (of 15) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) Category file for [composite] updated hit RETURN to continue --> Create/Update the color table for [composite]? (y/n) [n] n Edit the history file for [composite]? (y/n) [n] y ----------------------------(clear)-------------------------------- ENTER/CORRECT FILE HISTORY INFORMATION Map ID ... Tue Jul 15 21:04:10 1997 Title .... composite_______________________________________________________ Project .. PERMANENT Creator .. grass Maptype .. raster__________________________________________________________ Data source _______________________________________________________________ _______________________________________________________________ Data Description generated by r.patch___________________________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- ENTER/CORRECT FILE HISTORY COMMENTS _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _______________________________________________________________ _________________________________________________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- History file for [composite] updated hit RETURN to continue --> You have now created a completely documented GRASS raster map layer! Consolidation exercises Way back at the start of this Training Session we set out a list of objectives along with the GRASS commands that would be incountered. See if you can now link the commands with these objectives, using the the tables that are provided below. It might be useful to consider the way in which GRASS commands are named. d. commands relate to some aspect of display on the screen; g. commands are general file manipulation commands; i. commands relate to image precessing; and r. commands manipulate raster files in some way. Display 5 GRASS commands to be used g.list d.mon d.rast r.report d.legend d.erase d.colormode r.colors d.colortable d.histogram r.patch r.support g.remove exit Display 6 Objectives: GRASS command used: Investigate the contents of any map layer .......................... Obtain a list of the current map layers .......................... Update map layer support files .......................... Construct a compostie thematic map .......................... Display the map layers on the computer screen .......................... Tidying up: The raster map layers that you have generated during this Training Session will not be needed in subsequent Training Sessions. If you wish, you can conserve disk space by removing these map layers, but do this with care. First, type the command g.list to remind yourself of all the map layers now held on the disk. Remember, you are not deleting any of the files in mapset PERMANENT, but only the ones you have created in you own mapset. You can now delete any map layers you wish by using g.remove followed by a filename. If you do not wish to move on to the next chapter in this session you will have to close down GRASS. To do this, firstly you should close down any of the graphics windows you have opened during this session. GRASS 4.1 > d.mon stop=x0 Monitor 'x0' terminated GRASS 4.1 > d.mon stop=x1 Error - No such monitor as 'x1' Finally, to exit GRASS. GRASS 4.1 > exit ----------------------------(clear)-------------------------------- GRASS SESSION WRAPUP You have just finished working on mapset: The following RASTER maps belong to it: composite image popln segment topo contours landcov rail source urban crash plant roads spillage water There are no VECTOR maps in this mapset There are no SITES maps in this mapset Do you wish to selectively remove data files? y/n [n] y ----------------------------(clear)-------------------------------- REMOVE FACILITY This program allows you to remove files found in your mapset Please select the type of file to be removed 1 raster files 2 region definition files RETURN to exit > 1 enter raster file to be removed Enter 'list' for a list of existing raster files Hit RETURN to cancel request > composite Ok to remove [composite]? (y/n) y REMOVE [composite] raster header category color history misc enter raster file to be removed Enter 'list' for a list of existing raster files Hit RETURN to cancel request > ----------------------------(clear)-------------------------------- REMOVE FACILITY This program allows you to remove files found in your mapset Please select the type of file to be removed 1 raster files 2 region definition files RETURN to exit > ----------------------------(clear)-------------------------------- GOOD BYE from GRASS grass@congo[56]> Step Three: Working with information ================================================================================ B: Manipulating and extracting information Introduction GIS can be employed as data storage and management systems. In this role they are often used to extract summary data or simple statistics concerning thematic map layers. In this Training Session you will generate summary and statistical information for a number of the map layers that you have already become familar with. The Training Session is divided into a series of short exercises. The order in which you complete the exercises, and indeed the final choice of exercises undertaken, we leave to your discretion. However, it is advisable to complete the earlier exercises, and then select from the later exercises when you are more familiar with the techniques being used. Since you are now a little more familiar with the GRASS working environment and style of interaction, the level of discussion and explaination has been reduced slightly from that seen in Training Session A. Starting GRASS grass@congo[56]> grass4.1 ----------------------------(clear)-------------------------------- LOCATION: tutor_________ (enter list for a list of locations) MAPSET: PERMANENT_____ (or mapsets within a location) DATABASE: /usr/gis__________________________________________ ----------------------------(clear)-------------------------------- B: Manipulating and extracting information Aims: The aim of this Training Session is to use a variety of data extraction precesses available within GRASS to obtain summary and stastical information about thematic map layers. Objectives: On completion of this Training Session you will be able to :- Calculate the average elevation of woodland, pasture and scrub land Calculate the total area of woodland, pasture and scrub land Calculate the length of river flowing thtough arable and industrial land Calculate the population living within 1km of the motorway and railway Display 2 GRASS commands to be used r.average d.rast d.what.rast r.report r.reclass r.mapcalc d.zoom d.measure d.erase g.region r.buffer g.remove r.mask Display 3 Rastor map layers available for this tutorial. crash Site of motorway crash contours Topographic contour data landcov Land cover / land use popln Population data rail Railways segment Polluted river segment spillage Chemical spillage topo Degital elevation model water Rivers and reservoirs Exercise 1 Calculating the average elevation of woodland, pasture and scrub land. This exercise will involve using the grass command r.average which has the ability to construct statistical information from a map layer. r.average requires two raster map layers for processing. The 'base file' contains the catagories into which you wish to place the averages. The 'values file' contains the data you wish to average. So for example, suppose you had two layers - one that defined urban/non-urban areas, and another that contained population. To calculate the average urban population, the base file would be 'urban' and the values file 'population'. In this example you wish to obtain statistics from the topo map layer which contains information on height above sea level. However, the specific locations that you are interested in generating these statistics for are contained in the landcov map layer which holds information on land use. First use the command r.report to check the values and landuse categories that are used in the landcov map layer. Enter this information into the relevant columns of the table below along with the extent of each category. Display 4 Category Landuse Average Height Total Area ...... .......... .......... ........ ...... .......... .......... ........ ...... .......... .......... ........ ...... .......... .......... ........ ...... .......... .......... ........ GRASS 4.1 > r.report landcov units=kilom,percent r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Wed Jul 16 19:57:14 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: Land Cover / Land Use (landcov in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | square | % | |#|description |kilometers| cover| |-----------------------------------------------------------------------------| |1|Industry . . . . . . . . . . . . . . . . . . . . . . . . | 4.177500| 2.90| |2|Residential. . . . . . . . . . . . . . . . . . . . . . . | 16.580000| 11.51| |3|Quarry . . . . . . . . . . . . . . . . . . . . . . . . . | 1.387500| 0.96| |4|Woodland . . . . . . . . . . . . . . . . . . . . . . . . | 4.372500| 3.04| |5|Arable . . . . . . . . . . . . . . . . . . . . . . . . . | 57.282500| 39.78| |6|Pasture. . . . . . . . . . . . . . . . . . . . . . . . . | 42.040000| 29.19| |7|Scrub. . . . . . . . . . . . . . . . . . . . . . . . . . | 16.395000| 11.39| |8|Water. . . . . . . . . . . . . . . . . . . . . . . . . . | 1.765000| 1.23| |-----------------------------------------------------------------------------| |TOTAL |144.000000|100.00| +-----------------------------------------------------------------------------+ Now calculate the average height for each land use category. GRASS 4.1 > r.average base=landcov cover=topo output=avheight r.stats: 100% Display the map layer containing average elevation values. GRASS 4.1 > d.rast avheight Discussion: The image should look similar to the landcov map layer - but this timewe have a different color scheme. What does each color represent? We already know how to display the category values of a rastor map by using r.report. Let us now try a different method. GRASS 4.1 > d.what.rast avheight Buttons Left: what's here Right: quit Responses: Move the mouse over the image and click with the left button over any color. You will be presented with two lines similar to 455025(E) 314275(N) avheight in PERMANENT (8)99.269121813 Discussion: The first line refers to the geographical coordinates of the point on the raster map that you have clicked. The number in brackets on the second line refers to the category value (corresponding to the landuse type in this case), and the number afterwards, the category label (average height in this case). We are not limited to querying just one map layer at a time. Try querying avheight and landcov together. GRASS 4.1 > d.what.rast landcov,avheight Buttons Left: what's here Right: quit 444825(E) 317425(N) landcov in PERMANENT (8)Water avheight in PERMANENT (8)99.269121813 Check all eight colors noting down the results in the relevant column of Display 4. Click the right mouse button when you have finished. Discussion: 1 83.61 2 76.30 3 135.25 4 142.94 5 104.02 6 133.48 7 148.25 8 99.27 Check that your results are the same (to 2 decimal places) as those shown above. This table is a record of the average elevation, recorded in meters, of each of the eight land use categories contained in the landcov file. Are these values reasonable/sensible? Consider the possible effects of altitude on both the climate and vegitation. What other factors might be relevant? Exercise 2 Calculating the total length of river flowing through arable and industrial land. Suppose that you are a member of a committee that has been set up to explore the problem of river pollution in our study area. The major sources of pollutants are considered to be industrial waste, and agrochemicals leached from arable land. To gain an initial estimate of the magnitude of the problem it might be useful to know the total length of river sections that flow through land under arable and industrial use. This process involves three stagers - extracting the landuse types you are interested in; finding the rivers that flow through arable and industrial areas; and finding the distance covered by these sections of river. The first step then is to reclass the landuse raster map layer so that only arable and industrial landuses are included. GRASS 4.1 > r.reclass ----------------------------(clear)-------------------------------- RECLASS: Program for reassigning category codes Enter name of data layer to be reclassified Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > landcov Enter name of NEW RECLASSIFIED map Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > arab1_ind ----------------------------(clear)-------------------------------- Please indicate how you would like the reclass table initialized 0_ 0 All values set to zero 1 All values set to the same category number AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- ENTER NEW CATEGORY NUMBERS FOR THESE CATEGORIES OLD CATEGORY NAME OLD NEW NUM NUM Industry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 0____ Residential . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 0____ Quarry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 0____ Woodland . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 0____ Arable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 0____ Pasture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 0____ Scrub . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 0____ Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 0____ Next category: end__ (1 thru 8) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) A table of old category name and numbers appears with the cursor positioned on the Industry line. Type 1 and press Return four more times until the cursor is on the Arable line. Type 1 again and press You now need to give the new map layer a title and legend categories, so fill out the 'form' as follows (remember, you can overtype existing text with new text or spaces): ----------------------------(clear)-------------------------------- [arab1_ind] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: Arable and Industrial Landuse_______________________________ CAT NEW CATEGORY NAME NUM 0 Background__________________________________________________ 1 Industrial/Arable___________________________________________ Next category: end__ (of 1) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- You have just reclassified map [landcov in PERMANENT] into new map [arab1_ind in PERMANENT] A word of warning about r.reclass. All you have done by reclassifying a map layer is to change the category values of the original layer. To save on disk space and processing speed, GRASS does not create a completely new map layer when you do this. It simply creates a special 'reclass' file that stores the changes you have made to the original. As a consequence, if the original raster map layer is removed, so will all its reclassified maps. To make sure your reclassification has been sucessful, display the new file arab1_ind using d.rast. Do the same for the layer water which shows the distribution of rivers and reservoirs. The second step involves finding the rivers that flow through the arable and industrial areas. GRASS 4.1 > r.mapcalc mapcalc> polrivers = arab1_ind * water EXECUTING polrivers = ... 100% CREATING SUPPORT FILES FOR polrivers minimum value 0, maximum value 4 mapcalc> Discussion: r.mapcalc is a program that allows the manipulation of raster map layers using the concept of map algebra. That is, to treat map layers as quantities in mathmatical equations. These can be very simple as in the example above, or more sophisticated using a variety of mathmatical functions and operators. Map algebra provides an extremely powerful and flexible way of manipulating spatial information. In the example above, you have 'multiplied' each cell in the raster map layer arab1_ind by the value of the cell in the same postition in the map layer water. The resultant map indicates only the areas that have non-zero values on both map layers. Consider for a moment why this is so. If you multiply two layers which have cells coded with the values 0 and 1, only those cells which have a value of 1 in both layers will retain a value of 1 in the output file. ie, 1 x 1 = 1 the cell is 'retained' 1 x 0 = 0 the cell is 'rejected' 0 x 1 = 0 the cell is 'rejected' 0 x 0 = 0 the cell is 'rejected' This is a form of logical of boolean overlay (effetively the logical AND operation). You now need to find the length of the rivers that flow through potentially polluted areas. Without measuring the length directly, we can get some idea of the length by using r.report. GRASS 4.1 > r.report polrivers units=cell r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Thu Jul 17 02:00:21 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: arab1_ind*water (polrivers in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | cell| |#|description |count| |-----------------------------------------------------------------------------| |0|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . |57092| |4| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 508| |-----------------------------------------------------------------------------| |TOTAL |57600| +-----------------------------------------------------------------------------+ Discussion: The number of cells in the map polrivers should be displayed. An approximation of the length of these sections of river can be calculated by multiplying this number by 50 (remember each cell is 50m by 50m long). What approximations have been made in calculating river length? Could these be overcome? Do you think the final results in an overestimate or underestimate of the true value? A much better way of measuring linear features involves the use of d.measure. This command allows you to measure features in the active display frame on the graphics monitor. However, it would take quite a long time to measure all the river segments displayed in the polrivers map layer, so one segment has been selected for you. Display (using d.rast) the map layer segment. GRASS 4.1 > d.rast segment GRASS 4.1 > d.measure ----------------------------(clear)-------------------------------- Buttons: Left: where am i Middle: set FIRST vertex Right: quit this Point (with the mouse) to the start of the line to be measured and press the middle button on the mouse. ----------------------------(clear)-------------------------------- Left: where am i Middle: set NEXT vertex Right: FINISH LEN: 428.23 meters LEN: 1939.13 meters LEN: 2879.78 meters LEN: 3153.21 meters LEN: 3671.26 meters Direct the mouse around the line, pressing the middle button at bends in the line. Press the right button when you reach the end of the line. ----------------------------(clear)-------------------------------- Buttons: Left: DO ANOTHER Middle: Right: quit this LEN: 3671.26 meters AREA: 58.19 hectares 0.2247 square miles 581895.13 square meters Press the right button again, and note down the length of the line that you have just measured. You are now going to measure the length of the streatch of river in segment again twice, but this time at different scales, using the command d.zoom to 'zoom in' to the image. GRASS 4.1 > d.zoom Buttons: Left: Establish a corner Middle: Check coordinates Right: Accept region north: 321150 south: 318150 east: 447550 west: 444500 This region now saved as current region. Note: run 'd.erase' for the new region to affect the graphics. Point with the mouse to the left hand edge of the image, about halfway up the edge press the left button. Now point to the halfway point on the top edge, noting that a box is displayed as you move the mouse. When this box covers the top left quarter of the image press the right button. A red box will appear in the bottom left corner of the image, click on the YES option with the mouse to accept the new region. You are now required to enter d.erase in order for the zoomed section you just chose to be set to the graphics window. GRASS 4.1 > d.erase Now redisplay the segment image with d.rast. Note that the image is much enlarged - also notice that the raster cells that make up the image appear much coarser. GRASS 4.1 > d.rast segment Use d.measure again to measure the line representing the river segment. Note down the length measured. GRASS 4.1 > d.measure ----------------------------(clear)-------------------------------- Buttons: Left: where am i Middle: set FIRST vertex Right: quit this ----------------------------(clear)-------------------------------- Left: where am i Middle: set NEXT vertex Right: FINISH LEN: 115.95 meters LEN: 261.57 meters LEN: 339.71 meters LEN: 429.03 meters LEN: 555.88 meters LEN: 719.19 meters LEN: 890.08 meters LEN: 1016.41 meters LEN: 1096.64 meters LEN: 1205.83 meters LEN: 1382.06 meters LEN: 1472.37 meters LEN: 1569.04 meters LEN: 1754.95 meters LEN: 1842.65 meters LEN: 2048.48 meters LEN: 2471.38 meters LEN: 2606.80 meters LEN: 2695.88 meters LEN: 3019.96 meters LEN: 3088.32 meters LEN: 3191.28 meters LEN: 3280.36 meters LEN: 3413.15 meters LEN: 3511.18 meters LEN: 3614.15 meters LEN: 3852.15 meters LEN: 4213.68 meters ----------------------------(clear)-------------------------------- Buttons: Left: DO ANOTHER Middle: Right: quit this LEN: 4213.68 meters AREA: 130.08 hectares 0.5022 square miles 1300800.11 square meters Zoom in further, using d.zoom, so that the river segments' end points lie just inside the bottom left and top right corner of the box. Run d.erase and d.rast to display the river segment. Again use d.measure to measure the line, and note down the length. Work out the length using the r.report, as done before with all of polrivers. Compare this with the three lengths found using d.measure. GRASS 4.1 > r.report segment units=cell r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Thu Aug 7 18:06:14 1997| |-----------------------------------------------------------------------------| | north: 321150 east: 447550 | |REGION south: 318150 west: 444500 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: polrivers_segment (segment in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | cell| |#|description |count| |-----------------------------------------------------------------------------| |0|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 3553| |1| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 107| |-----------------------------------------------------------------------------| |TOTAL | 3660| +-----------------------------------------------------------------------------+ Do you notice any trend in the lengths found using d.measure, why are they different from each other? Is there a 'correct' length for the river segment in the segment map? Finally, you need to 'reset' the view of the data back to the original dimentions. (You will look at this in more detail in Section C.) GRASS 4.1 > g.region -d GRASS 4.1 > d.erase Exercise 3 Calculating the 'at-risk' population in the event of an industrial chemical spillage. GIS are potentially very valuable as decision support systems, and a particularly attractive facility is "what-if" scenerio modelling. Such systems and techniques are already being used in several parts of the country for emergency planning and evacuation procedures. With the limited data set available in the tutorial it is only possible to present a very simplistic example, but this will serve to illustrate the basic principles. Suppose there is a need to transport a dangerous chemical cargo through the study area. It is recognised that any accidental spillage would endanger the publlic and necessitate an immediate evacuation of the population living within 1000m of the accident site. You have the choice of transporting the cargo affecting your decision will be, for each of the modes of transport, the population potentially at risk should an accident occur. Let us now use GRASS to answer this simple query. The problem is essentially one of estimating the total population living within 1000m of the railway or M1 motorway. There are two steps needed to tackle this problem. Firstly, you need to generate new map layers, indicating the land areas within 1000m of these transport systems. Secondly, you need to sum the total population within each of these areas from the popln map layer. Step 1 First generate the map layers showing the land area within 1000m of the railway and the M1. This is achieved using the r.buffer command which generates a file in which each cell can be related to the Euclidean (or "as the crow flies") distance from the nearest non-zero cell in the input file. GRASS 4.1 > r.buffer rail distances=1000 output=railbuf_temp Reading input map (rail) ... 100% Finding buffer zones ... 100% Writing output map (railbuf_temp) ... 100% We need to repeat this procedure for the M1 motorway. Therefore, our next step is to identify the motorway from the roads map layer. GRASS 4.1 > r.reclass ----------------------------(clear)-------------------------------- RECLASS: Program for reassigning category codes Enter name of data layer to be reclassified Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > roads Enter name of NEW RECLASSIFIED map Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > mway ----------------------------(clear)-------------------------------- Please indicate how you would like the reclass table initialized 0_ 0 All values set to zero 1 All values set to the same category number AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) You should be a little more familiar with reclassifying maps by now. Make sure all catagories are zero except for motorways which is 1. Give the map layer a suitable title. ----------------------------(clear)-------------------------------- ENTER NEW CATEGORY NUMBERS FOR THESE CATEGORIES OLD CATEGORY NAME OLD NEW NUM NUM Motorway . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 0____ Proposed . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 0____ 'A' roads . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 0____ 'B' roads . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 0____ 'C' roads . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 0____ Next category: end__ (3 thru 12) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- [mway] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: M1 Motorway_________________________________________________ CAT NEW CATEGORY NAME NUM 0 Background__________________________________________________ 1 Motorway____________________________________________________ Next category: end__ (of 1) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- You have just reclassified map [roads in PERMANENT] into new map [mway in PERMANENT] Next calculate the 'buffer zone' around the motorway. GRASS 4.1 > r.buffer mway distances=1000 output=mwaybuf_temp Reading input map (mway) ... 100% Finding buffer zones ... 100% Writing output map (mwaybuf_temp) ... 100% Display the new map layers railbuf_temp and mwaybuf_temp to check that they look sensible. GRASS 4.1 > d.rast railbuf_temp GRASS 4.1 > d.rast mwaybuf_temp Discussion: Creating a 'buffer' around a geographical object is a common process in GIS where the sphere of influence of the object extends beyond its own boundaries. Note that both buffered map layers distinguish between the original object and the buffer around it. For the purposes of this exercise, we do not need to make this distinction. GRASS 4.1 > r.reclass ----------------------------(clear)-------------------------------- RECLASS: Program for reassigning category codes Enter name of data layer to be reclassified Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > mwaybuf_temp Enter name of NEW RECLASSIFIED map Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > mwaybuf ----------------------------(clear)-------------------------------- Please indicate how you would like the reclass table initialized 0_ 0 All values set to zero 1 All values set to the same category number AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- ENTER NEW CATEGORY NUMBERS FOR THESE CATEGORIES OLD CATEGORY NAME OLD NEW NUM NUM distances calculated from these locations . . . . . . . . . . . . 1 1____ 1000 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1____ Next category: end__ (1 thru 2) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- [mwaybuf] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: M1 Motorway Buffer__________________________________________ CAT NEW CATEGORY NAME NUM 0 Background__________________________________________________ 1 1000m_around_M1_____________________________________________ Next category: end__ (of 1) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- You have just reclassified map [mwaybuf_temp in PERMANENT] into new map [mwaybuf in PERMANENT] Repeat the same procedure for railbuf_temp calling the new layer railbuf. This time reclassify the railway and buffer as 2. GRASS 4.1 > r.reclass ----------------------------(clear)-------------------------------- RECLASS: Program for reassigning category codes Enter name of data layer to be reclassified Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > railbuf_temp Enter name of NEW RECLASSIFIED map Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > railbuf ----------------------------(clear)-------------------------------- Please indicate how you would like the reclass table initialized 0_ 0 All values set to zero 1 All values set to the same category number AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- ENTER NEW CATEGORY NUMBERS FOR THESE CATEGORIES OLD CATEGORY NAME OLD NEW NUM NUM distances calculated from these locations . . . . . . . . . . . . 1 2____ 1000 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2____ Next category: end__ (1 thru 2) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- [railbuf] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: Railway Buffer______________________________________________ CAT NEW CATEGORY NAME NUM 0 Background__________________________________________________ 1 Background__________________________________________________ 2 1000m_around_rail___________________________________________ Next category: end__ (of 2) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- You have just reclassified map [railbuf_temp in PERMANENT] into new map [railbuf in PERMANENT] Before calculating the number of people living within the buffer, it makes the process a little easier if you first combine the two buffers into a single file. This is not essential, and in fact it is often the case that there are many ways to obtain the same objectives in both GRASS and, indeed, in most other GIS systems. The trick (some people call it skill) is in knowing the most efficient methods. GRASS 4.1 > r.mapcalc mapcalc> buffers = mwaybuf + railbuf EXECUTING buffers = ... 100% CREATING SUPPORT FILES FOR buffers minimum value 0, maximum value 2 mapcalc> Display the buffers map layer to sheck you have the expected result before continuing. GRASS 4.1 > d.rast buffers You are now ready to calculate the population contained in each of the buffer zones. GRASS 4.1 > r.mapcalc mapcalc> at_risk = buffers * popln EXECUTING at_risk = ... 100% CREATING SUPPORT FILES FOR at_risk minimum value 0, maximum value 2 mapcalc> Display this final map layer to get a visual impression of those people potentially 'at risk' around each transport link. GRASS 4.1 > d.rast at_risk Finally, produce a quantitative comparison of the numbers of people 'at risk' in each group. GRASS 4.1 > r.report at_risk units=cell r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Thu Aug 7 19:35:03 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: buffers*popln (at_risk in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | cell| |#|description |count| |-----------------------------------------------------------------------------| |0|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . |55367| |1| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 838| |2| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 1395| |-----------------------------------------------------------------------------| |TOTAL |57600| +-----------------------------------------------------------------------------+ Based on this evidence only, which mode of transport would be preferable for moving the dangerous chemical cargo through our study area? Display the rail map layer again. Does anything about this layer indicate that the results found should be treated with caution? GRASS 4.1 > d.rast rail User Notes:............................................................... Step 2 Let us suppose, that for reasons of cost, the motorway route was decided upon. Unfortunately, the lorry transporting the chemicals wa involved in a head on collision, leading to a spillage and fire. You are now required to identify the population at risk from the pollution, and evacuate everyone living within 1000m of the crash. The raster map crash shows the point of the collision. Create a buffer of 1000m around the crash site, by using the r.buffer command. GRASS 4.1 > r.buffer crash distances=1000 output=crashbuf We will now limit our interrogation of the database to this buffer zone, using the r.mask command GRASS 4.1 > r.mask ----------------------------(clear)-------------------------------- MASK: Program for managing current GIS mask current mask: none Options: 1 Remove the current mask 2 Identify a new mask RETURN Exit program > 2 Enter name of data layer to be used for mask Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > crashbuf ----------------------------(clear)-------------------------------- IDENTIFY THOSE CATEGORIES TO BE INCLUDED IN THE MASK OLD CATEGORY NAME CAT NUM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0____ distances calculated from these locations . . . . . . . . . . . . 1 1____ 1000 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1____ Next category: end__ (0 thru 2) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) Now a form appears similar to that used by the r.reclass command. Enter a 1 next to catagories 1 and 2 (distances calculated from these locations and 1000 meters respectively). ----------------------------(clear)-------------------------------- MASK: Program for managing current GIS mask current mask: in mapset masking category(ies): 1-2 Options: 1 Remove the current mask 2 Identify a new mask RETURN Exit program > Discussion: The r.mask command creates a special raster map layer, called MASK, that allows the user to block out certain areas from the current geographic region, by "hiding" the areas not wanted from the sight of other GRASS programs. In this case areas outside the 1000m buffer around the site of the crash are not wanted - these have a value of 0 in MASK and are hidden, whilst areas within the buffer have a value of 1. This can be shown by displaying any map layer using d.rast (eg. try displaying landcov). Now you have your MASK use r.report to find the number of populated cells that need to be evacuated. GRASS 4.1 > r.report popln units=cell r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Thu Aug 7 19:54:14 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:crashbuf in PERMANENT, categories 1-2 | |-----------------------------------------------------------------------------| |MAP: Population Data (popln in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | cell| |#|description |count| |-----------------------------------------------------------------------------| |0|Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 1199| |1|Populated cells. . . . . . . . . . . . . . . . . . . . . . . . . . . | 58| |-----------------------------------------------------------------------------| |TOTAL | 1257| +-----------------------------------------------------------------------------+ The procedure you have just carried out was used to evacuate people from the area around the crash. However, the resulting smoke pollution was spread far to the north east by strong south westerly winds. The extent of the pollution was not known untill the next day, when it became clear that areas had been affected that had not been evacuated. The extent of the pollution is shown in the spillage map layer. (You have to remove the MASK in order to allow the whole image to be shown). GRASS 4.1 > g.remove rast=MASK REMOVE [MASK] raster header category color MISSING history misc GRASS 4.1 > d.rast spillage The effects of the pollution could affect the health of those exposed to over 5 units of the pollutant and therefore they must ahve a medical check up. You need to find the amount of people affected by the pollution. This can be done simply, by using the r.mask command. GRASS 4.1 > r.mask Repeat the same procedure as before, except enter spillage as the data layer to be used for the mask, and enter 1 for categories 5 and over. ----------------------------(clear)-------------------------------- MASK: Program for managing current GIS mask current mask: none Options: 1 Remove the current mask 2 Identify a new mask RETURN Exit program > 2 Enter name of data layer to be used for mask Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > spillage ----------------------------(clear)-------------------------------- IDENTIFY THOSE CATEGORIES TO BE INCLUDED IN THE MASK OLD CATEGORY NAME CAT NUM No risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 0____ Low risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 0____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1____ Next category: 10___ (0 thru 30) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- IDENTIFY THOSE CATEGORIES TO BE INCLUDED IN THE MASK OLD CATEGORY NAME CAT NUM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1____ Moderate risk . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1____ Next category: 20___ (0 thru 30) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- IDENTIFY THOSE CATEGORIES TO BE INCLUDED IN THE MASK OLD CATEGORY NAME CAT NUM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1____ High risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1____ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1____ Very high risk . . . . . . . . . . . . . . . . . . . . . . . . . 29 1____ Next category: 30___ (0 thru 30) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- IDENTIFY THOSE CATEGORIES TO BE INCLUDED IN THE MASK OLD CATEGORY NAME CAT NUM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1____ Next category: end__ (0 thru 30) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- MASK: Program for managing current GIS mask current mask: in mapset masking category(ies): 5-30 Options: 1 Remove the current mask 2 Identify a new mask RETURN Exit program > GRASS 4.1 > r.report popln units=cell r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Thu Aug 7 20:14:27 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:spillage in PERMANENT, categories 5-30 | |-----------------------------------------------------------------------------| |MAP: Population Data (popln in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | cell| |#|description |count| |-----------------------------------------------------------------------------| |0|Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 8288| |1|Populated cells. . . . . . . . . . . . . . . . . . . . . . . . . . . | 461| |-----------------------------------------------------------------------------| |TOTAL | 8749| +-----------------------------------------------------------------------------+ This gives the numberof populated cells put into risk by the pollution, however it doesn't take into account the cells that were evacuated. In order to do this you must first use r.reclass for crashbuf as shown on page B-13 to reclass categories 1 and 2 to 1, as a new raster crashbuf_temp. Then use r.mapcalc to take this area away from the MASK. GRASS 4.1 > r.reclass ----------------------------(clear)-------------------------------- RECLASS: Program for reassigning category codes Enter name of data layer to be reclassified Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > crashbuf Enter name of NEW RECLASSIFIED map Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > crashbuf_temp ----------------------------(clear)-------------------------------- Please indicate how you would like the reclass table initialized 0_ 0 All values set to zero 1 All values set to the same category number AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- ENTER NEW CATEGORY NUMBERS FOR THESE CATEGORIES OLD CATEGORY NAME OLD NEW NUM NUM distances calculated from these locations . . . . . . . . . . . . 1 1____ 1000 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1____ Next category: end__ (1 thru 2) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- [crashbuf_temp] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: Crash Buffer________________________________________________ CAT NEW CATEGORY NAME NUM 0 Background__________________________________________________ 1 100m_around_crash___________________________________________ Next category: end__ (of 1) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- You have just reclassified map [crashbuf in PERMANENT] into new map [crashbuf_temp in PERMANENT] GRASS 4.1 > r.mapcalc mapcalc> MASK = MASK - crashbuf_temp MASK - already exists. ok to overwrite? (y/n) y EXECUTING MASK = ... 100% CREATING SUPPORT FILES FOR MASK minimum value 0, maximum value 1 mapcalc> You have now modified the MASK to exclude those populated cells that have been evacuated. Note that you can manipulate MASK with r.mapcalc and use g.remove to remove it, just like any other raster map layer. Finally work out how many people need a health check following the chemical spillage. Use r.report to find out the number of populated cells affected, then multiply this value by 5.4 (the average population within each populated cell). GRASS 4.1 > r.report popln units=cell r.stats: 100% +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: tutor Thu Aug 7 21:06:58 1997| |-----------------------------------------------------------------------------| | north: 322000 east: 456000 | |REGION south: 310000 west: 444000 | | res: 50 res: 50 | |-----------------------------------------------------------------------------| |MASK:MASK in PERMANENT | |-----------------------------------------------------------------------------| |MAP: Population Data (popln in PERMANENT) | |-----------------------------------------------------------------------------| | Category Information | cell| |#|description |count| |-----------------------------------------------------------------------------| |0|Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 7539| |1|Populated cells. . . . . . . . . . . . . . . . . . . . . . . . . . . | 452| |-----------------------------------------------------------------------------| |TOTAL | 7991| +-----------------------------------------------------------------------------+ Consolidation exercises (1) Determine the highest altitude that is attained by travelling along the M1 motorway within the study area. Hints: This can be answered by using the r.mapcalc command to generate the result, and r.report to view it. (2) Returning to Exercise 2 on page B-7, it could be argued that it would be more relevant to know the area of arable and industrial land within, say, 250m of all rivers and reservoirs. There are several ways in which this information can be obtained using GRASS. It is possible to derive a solution using only the commands that you are currently familiar with. Does this alternative method overcome the problems previously encuntered in finding the river lengths? Tidying up: The raster map layers that you have generated during this Training Session will not be needed in future sessions. If you wish, to conserve disk space, you have the option of removing these map layers, but do this with care. First type the command g.list to remind yourself of all the map layers now held on the disk. Remember, you are not deleting any of the files in mapset PERMANENT, but only the ones you have created in your own mapset. To close down any of the graphics windows you have opened during this session. GRASS 4.1 > d.mon stop=x0 Finally, exit GRASS. GRASS 4.1 > exit ----------------------------(clear)-------------------------------- GRASS SESSION WRAPUP You have just finished working on mapset: The following RASTER maps belong to it: MASK crash mwaybuf railbuf topo arab1_ind crashbuf mwaybuf_temp railbuf_temp urban at_risk crashbuf_temp plant roads water avheight image polrivers segment buffers landcov popln source contours mway rail spillage There are no VECTOR maps in this mapset There are no SITES maps in this mapset Do you wish to selectively remove data files? y/n [n] ----------------------------(clear)-------------------------------- GOOD BYE from GRASS C: Optimal routing ================================================================================ Introduction So far, most of the applications of GIS that you have been exposed to have concentrated on processing parcels of land, of area based festures. In this exercise a simple example is presented which serves to demonstrate the potential of GIS for processing linear features. Suppose that a new sewage pipeline is to be laid between a village and a sewage treatment plant. You are required to find the cheapest route for the pipeline. The 'cost' of the pipeline is to be based, rather simplistically, on the land cover over which it passes. In addition, it is desirable that the pipeline should run 'down-hill' as much as possible to minimise the need for the installation of expensive pumping equipment. Starting GRASS If you have not left GRASS since the previous Training Session, you can ignore this page, otherwise, you must first start GRASS and tell it which dataset you are going to work with. To start GRASS, type grass4.1 and press the return key. You should be presented with a screen similar to the following: grass@congo[2]> grass4.1 ----------------------------(clear)-------------------------------- GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: tutor_________ (enter list for a list of locations) MAPSET: PERMANENT_____ (or mapsets within a location) DATABASE: /usr/gis__________________________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) You now need to tell GRASS which data set you will be working with. The opening screen prompts you to identify the LOCATION, MAPSET, and DATABASE. If you have just completed the previous Training Session, you will probably find that these are already set up correctly. If this is the case, just press to start. LOCATION stores the geographical extent of the study region. For this session you will continue to work with the Leicestershire data. Press the Return key a few times until the cursor is positioned on the LOCATION line. Then type leics to identify the study region. Don't worry about overtyping the default setting - just make sure that any remaining letters are deleted by pressing the space bar. Press Return to move the cursor onto the MAPSET line. Remember, a mapset in GRASS is your own set of files that reside in the current LOCATION. You should be using your user name. This will ensure that any files you alter or create will not be inadvertantly changed by anyone else. The DATABASE should not need changing on this occasion. It is simply the UNIX directory in which the sample locations and mapsets are stored. When you have set the correct options, press to start. ----------------------------(clear)-------------------------------- Welcome to GRASS 4.1 (Spring 1993) Update package 2 Geographic Resources Analysis Support System (GRASS) is a Trademark of U.S. Army Construction Engineering Research Laboratories (USACERL) New releases of GRASS are coordinated and produced by the Office of GRASS Integration (OGI) located at USACERL. This version running thru the C Shell (/bin/csh) Help is available with the command: g.help When ready to quit enter: exit Mapset in Location To use graphics in GRASS, you must open at least one graphics window. Therefore before you start this session, type: GRASS 4.1 > d.mon start=x0 select=x0 C: Optimal routing Aims: This training session illustrates the use of GIS techniques for finding an optimal path solution based on an assessment of the cost incurred in traversing the land surface. Objectives: On completion of this training session you will be able to :- Reduce the bounds and resolution fo the study area. Create an aspect map from a digital elevation model. Assign relative costs based on land use and direction of slope. Create a 'cost-distance' surface. Determine the least-cost pathway. Improve the presentation of the final map. Display the final map using a perspective view. Display 2 GRASS commands to be used g.region d.frame d.rast r.slope.aspect d.his r.rescale r.mapcalc r.reclass i.grey.scale r.cost r.drain r.patch d.colormode d.erase d.colors d.3d Display 3 Map layers to be used landcov Land cover / land use plant Sewage plant source Sewage source topo Digital elevation model Step 1: Selecting a sub-region from an original map layer. In this exercise you will use a small sub-region of the full study area employed in Training Sessions A and B. Only the south-west quadrant of the full study region is required since both the sewage source and the sewage treatment plant are located within this sub-region. Processing smaller files leads to greater speed and generally increased efficiency. You have already used a sub-region, in exercise 2, B-10, when d.zoom was used to "zoom-in" on an area of intrest. The g.region command you will now use is much more flexible, with many options, to give you greater control in setting the geographic region and grid resolution of your study area. GRASS 4.1 > g.region ----------------------------(clear)-------------------------------- REGION FACILITY LOCATION: tutor MAPSET: PERMANENT CURRENT REGION: N=322000 S=310000 RES=50 ROWS=240 E=456000 W=444000 RES=50 COLS=240 PROJECTION: 1 (UTM) ZONE: 0 Please select one of the following options Current Region Region Database 1 Modify current region directly 6 Save current region in the database 2 Set from default region 7 Create a new region 3 Set from a database region 8 Modify an existing region 4 Set from a raster map 9 Set from 3d.view file 5 Set from a vector map RETURN to quit > 1 ----------------------------(clear)-------------------------------- IDENTIFY REGION ============================= DEFAULT REGION ======== | Default North:322000 | | | | ======= YOUR REGION ======= | | | NORTH EDGE:316000_____ | | | | | | Def. West |WEST EDGE | |EAST EDGE | Def. East 444000 |444000_____| |450000_____| 456000 | | SOUTH EDGE:310000_____ | | | ============================= | | | | Default South:310000 | ===================================================== PROJECTION: 1 (UTM) ZONE: 0 Default GRID RESOLUTION Region 50 --- East-West --- 50________ 50 -- North-South -- 50________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- projection: 1 (UTM) zone: 0 north: 316000 south: 310000 east: 450000 west: 444000 e-w res: 50 n-s res: 50 total rows: 120 total cols: 120 total cells: 14,400 Do you accept this region? (y/n) [y] > ----------------------------(clear)-------------------------------- REGION FACILITY LOCATION: tutor MAPSET: PERMANENT CURRENT REGION: N=316000 S=310000 RES=50 ROWS=120 E=450000 W=444000 RES=50 COLS=120 PROJECTION: 1 (UTM) ZONE: 0 Please select one of the following options Current Region Region Database 1 Modify current region directly 6 Save current region in the database 2 Set from default region 7 Create a new region 3 Set from a database region 8 Modify an existing region 4 Set from a raster map 9 Set from 3d.view file 5 Set from a vector map RETURN to quit > GRASS 4.1 > d.frame -e Discussion: This process of 'windowing' in on a sub-region is a fundimental operation in most GIS. In GRASS each location is associated with a default set of bounds and grid resolution. g.region allows us to modify the extent of the region for subsequent operations. In the example above, we have retained the same grid resolution, but given the national grid coordinates of the SW quadrant of the study area. The command d.frame -e clears the graphics window and readjusts output for the new geographic window. Check that you have sucessfully determined the new region by displaying some of the rastor maps you are already familiar with (eg image and landcov). Step 2: Display the sewage source and treatment plant. Before continuing, it would be nice to see the location of the sewage source and treatment plants. These are held in the map layers source and plant respectively. Use d.rast with the overlay option to display the two sites in their geographical context. GRASS 4.1 > d.rast image Did you remember to reset the color table of the image back to grey in Training Session A? If you forgot, use r.colors image color=grey to do this. GRASS 4.1 > d.rast -o source GRASS 4.1 > d.rast -o plant Discussion: You should be able to see two red dots marking the location of the sewage source and plant. Point data of this sort are common in spatial databases. It is obviously rather inefficient to store a whole raster map layer just to represent one point. GRASS has the ability to manipulate special 'sites files' which are specifically designed for the processing such point data. However, for the purposes of this exercise, you will stick to representing our point data as raster map layers. Confirm which point represents the sewage source and which the sewage treatment plant by using d.what.rast. Click on both points and note down their geographical coordinates in the table below. Refer back to B-5 if you can't remember how to use this command. Display 4 Easting Northing source ....... ....... plant ....... ....... User Notes:................................................................ Step 3: Creating a new maplayer of slope direction (aspect). In this example the installation cost of the pipeline will be largely determined by the land uses it crosses. However, in order to minimise the future 'running costs', it is desirable to select a route that avoids up- hill sections whenever possible, since these would necessitate the installation and use of expensive pumping equipment. Since the source of the sewage is approximately to the south of the treatment plant, any slopes facing east, through south, and round to the west can be taken as up-hill sections. Clearly this is both an approximate and simplistic model, but it will serve to illustrate the basic principles of site selection. GRASS is able to compute both slope magnitude and direction (ie aspect) from a digital elevation model. GRASS 4.1 > r.slope.aspect elevation=topo aspect=asp1 percent complete: 100% CREATING SUPPORT FILES ELEVATION PRODUCTS for mapset [PERMANENT] in [tutor] 360 categories of aspect Color table for [asp1] set to aspect ASPECT [asp1] COMPLETE Discussion: The output file asp1 contains the direction of local slope at each cell location. These data are actual degree intervals measured from the East (catagory 1) moving anticlockwise. If a cell has no local slope, the aspect cannot be defined, and these areas are assigned a value of 0. Use d.rast to display the newly created file asp1. GRASS 4.1 > d.rast asp1 Can you interpret the meaning of the grey scale that you see? Try changing the color palette (using r.colors). Do any of the alternative color palettes aid in interpretation? GRASS 4.1 > r.colors asp1 color=grey GRASS 4.1 > r.colors asp1 color=grey.eq GRASS 4.1 > r.colors asp1 color=gyr GRASS 4.1 > r.colors asp1 color=rainbow GRASS 4.1 > r.colors asp1 color=ramp GRASS 4.1 > r.colors asp1 color=ryg GRASS 4.1 > r.colors asp1 color=random GRASS 4.1 > r.colors asp1 color=wave GRASS 4.1 > d.rast asp1 Finally change back to the original color palette buy using GRASS 4.1 > r.colors asp1 color=aspect Aspect maps can be used in combination with other map layers to give quite sophisticated graphical displays, where the aspect map adds a sense of depth representing the topography, to a thematic map. For example try the following: GRASS 4.1 > d.rast landcov GRASS 4.1 > d.his h_map=landcov i_map=asp1 100% The d.his command combines the color (or hue) of landcov, and the shading (or intensity) of the aspect map asp1. User Notes:............................................................... Step 4: Assigning costs due to slope direction. You are now in position to assign a cost for laying the pipeline through any given cell, on the basis of its direction of slope. This can be acheived using a combination of the r.rescale commands and r.mapcalc commands, this essentially reclasses the asp1 file's values from degree intervals to either north or south sloping. r.reclass isn't used because of the large amount of categories that need to be reclassed in asp1, the method below is much more efficient. GRASS 4.1 > r.rescale input=asp1 from=0,180 output=asp_north to=1,1 Rescale asp1[0,180] to asp_north[1,1] This has produced a new raster map layer of northerly slopes and flat areas, with a category value of 1. A map layer consisting of southerly facing slopes, with a category value of 10, can be produced as follows. GRASS 4.1 > r.rescale input=asp1 from=181,360 output=asp_south to=10,10 Rescale asp1[181,360] to asp_south[10,10] Now use r.mapcalc to produce a map layer that contains both the northerly slope and southerly slope categories (values of 1 and 10 respectively). GRASS 4.1 > r.mapcalc mapcalc> asp_cost = asp_north + asp_south EXECUTING asp_cost = ... 100% CREATING SUPPORT FILES FOR asp_cost minimum value 0, maximum value 10 mapcalc> Use d.rast to view the map layer that you have just created, make sure that only northerly slopes or flat areas have a category of 1, and southerly slopes a category of 10, in asp_cost, by exploring the map with d.what.rast asp1,asp_cost. GRASS 4.1 > d.rast asp_cost GRASS 4.1 > d.what.rast asp1,asp_cost Buttons Left: what's here Right: quit 447825(E) 312325(N) asp1 in PERMANENT (18)18 degrees from east asp_cost in PERMANENT (1) 447175(E) 311975(N) asp1 in PERMANENT (0)no aspect asp_cost in PERMANENT (0)no data 447375(E) 313225(N) asp1 in PERMANENT (247)247 degrees from east asp_cost in PERMANENT (10) Discussion: The 'base' cost for passing the pipeline through a cell has been assigned a unit value of 1. Thus those cells that incur no penalty (ie those with zero or down-hill slope) are assigned this cost rating. The cells containing up-hill slopes have been assigned a relative cost of 10 units. This rating is simply an arbitrary choice and furthermore takes no account of the magnitude of slope. In a real application weightings would be based on more sophisticated factual information. User Notes:............................................................... Step 5: Assigning costs due to land use. The cost of the sewage pipeline installation will also depend upon the land cover which it lies. These costings can be generated by re-assigning values in the landcov map layer. GRASS 4.1 > r.reclass ----------------------------(clear)-------------------------------- RECLASS: Program for reassigning category codes Enter name of data layer to be reclassified Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > landcov Enter name of NEW RECLASSIFIED map Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > lan_cost ----------------------------(clear)-------------------------------- Please indicate how you would like the reclass table initialized 0_ 0 All values set to zero 1 All values set to the same category number AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- ENTER NEW CATEGORY NUMBERS FOR THESE CATEGORIES OLD CATEGORY NAME OLD NEW NUM NUM Industry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 8____ Residential . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 8____ Quarry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1000_ Woodland . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 8____ Arable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 4____ Pasture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1____ Scrub . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1____ Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1000_ Next category: end__ (1 thru 8) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) When you have done this press and give the layer a suitable title. To save looking at all the categories (remember the maximum category is now 1000), press until the cursor is on the line 'Next category'. Type end and press . ----------------------------(clear)-------------------------------- [lan_cost] ENTER NEW CATEGORY NAMES FOR THESE CATEGORIES TITLE: Cost rating I think_________________________________________ CAT NEW CATEGORY NAME NUM 0 ____________________________________________________________ 1 Scrub_______________________________________________________ 2 ____________________________________________________________ 3 ____________________________________________________________ 4 Arable______________________________________________________ 5 ____________________________________________________________ 6 ____________________________________________________________ 7 ____________________________________________________________ 8 Woodland____________________________________________________ 9 ____________________________________________________________ Next category: 10___ (of 1000) AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) ----------------------------(clear)-------------------------------- You have just reclassified map [landcov in PERMANENT] into new map [lan_cost in PERMANENT] Discussion: Once again the 'cost ratings' selected for this exercise are arbitrary, but the following arguments are given by way of justification. First, pasture and scrub land represents areas with no special difficulty and are assigned the 'base' cost of 1 unit. Arable land is given a relative cost of 4 units, since it is likely that some compensation would be necessary to land owners for the loss of crops during the instalation process. Woodland is given an even higher cost rating (8 units) mainly because the access for pipe laying machinery would be severely impaired. Felling of trees to gain access would slow progress and may necessitate compensation payments to land owners. Finally, quarries and water bodies must be avoided, since it is clearly impossible for a pipeline to pass through these areas easily and cheaply. By assigning an extremely large cost to these pixels they are effectively rendered as impassable barriers. Step 6: Obtaining a combined cost surface. The final cost associated with any cell location can now be generated by combining (simply adding on a cell-by-cell basis) the seperate slope and land cover costings. We can use the r.mapcalc command to do this. GRASS 4.1 > r.mapcalc mapcalc> cost = lan_cost + asp_cost EXECUTING cost = ... 100% CREATING SUPPORT FILES FOR cost minimum value 0, maximum value 1001 mapcalc> Use d.rast to view the resulting map layer cost. GRASS 4.1 > d.rast cost Discussion: You will see that there are only two colors shown on the image - red and yellow. This is because most of the variation in cost occurs in classes under 20. Since the quarried area has a class of 1000, GRASS is unable to scale its color table equally among classes. To produce a more effective image, it is possible to 'histogram equalise' the map layer. GRASS 4.1 > i.grey.scale which layer needs a grey scale? Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > cost Reading cost ... 100% [cost in PERMANENT] now has a grey scale color table which layer needs a grey scale? Enter 'list' for a list of existing raster files Enter 'list -f' for a list with titles Hit RETURN to cancel request > <> GRASS 4.1 > d.rast cost Discussion: The image should show all the classes in different shades of grey. Histogram equalization involves assigning all the available colors as equally as possible to all the available classes, regardless of the class magnitude. The image itself represents the cost of laying the sewage pipeline through any individual cell within the study area. Explore the relationship between aspect, landcover and this cost surface by using the command GRASS 4.1 > d.what.rast cost,landcov,asp1 Buttons Left: what's here Right: quit 446975(E) 312075(N) cost in PERMANENT (0)no data landcov in PERMANENT (0)Background asp1 in PERMANENT (0)no aspect 448025(E) 313225(N) cost in PERMANENT (14) landcov in PERMANENT (5)Arable asp1 in PERMANENT (239)239 degrees from east 448175(E) 313475(N) cost in PERMANENT (2) landcov in PERMANENT (7)Scrub asp1 in PERMANENT (45)north of east User Notes:................................................................ Step 7: Creating a 'cost-distance' surface You will now produce a map layer which records the cumulative cost incurred by moving out radially, cell-by-cell, from the location of the sewage source. This 'cost-distance' surface is similar to the map layers generated by the r.buffer command. However, the cell values increase not simply due to distance, but also accumulate the specific cost traversing each successive pixel. A 'cost-distance' surface is generated in GRASS by using the r.cost command. To use the r.cost command, you are required to give the geographical coordinates of the point from which you wish to calculate the cost-distance from. Refer back to Display 4 and check that the coordinates of the sewage source are the same as those given below. GRASS 4.1 > r.cost input=cost coord=445775,310875 output=cost_dist Use d.rast to view the map layer that you have just created. GRASS 4.1 > d.rast cost_dist Discussion: The area of the image that are colored using darker legend colors (blues and reds) have accumulated high costs in travelling from the source location and are likely to be avoided when determining the route of the pipeline. Areas that are lighter colored (yellows and greens) are 'cheaper' to reach and therefore should provide the preferred route. You will explore the relationship between the optimal route and this cost_distance surface later on. User Notes:............................................................... Step 8: Determination of the least-cost route. You can now proceed with the next step in which you will use the GRASS command r.drain to trace the lowest cumulative cost route from the source cell to the target cell (the sewage works) on the cost-distance surface cost_dist. Check that the geographical coordinates of the sewage plant are the same as those you recorded in Display 4. GRASS 4.1 > r.drain input=cost_dist coord=444525,313875 output=pipeline Display the pipeline map layer using the d.rast command. GRASS 4.1 > d.rast pipeline Discussion: This shows the optimal route determined by GRASS for the proposed pipeline, based on the costings that were assigned to land cover and slope aspect categories. You may be wondering why the command we used to calculate the optimal route is called r.drain. The same command can be used to calculate the flow of a stream over a landscape. The downward draining of water over a surface is itself following an 'optimal route'. In your case, instead of just using slope as a determinant of flow path, you have added the addtional constraints imposed by land use. The optimal route of a stream differs from that of the pipeline in one important respect. Can you identify this, and what effect will it have on the suggested pipeline? User Notes:................................................................ Step 9: Improving the presentation of the final map. Now that you have identified a route for the sewage pipeline it will be useful to locate it relative to the different land use types. This is easily achieved by overlaying the pipeline onto the landcov map layer. However, the category code that is assigned to the pipeline must first be changed to prevent conflict with those that are already being used in the landcov map layer. GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > pipeline Enter name of NEW RECLASSIFIED map > pipeline2 Reclassify the old value of 1 to 9, and add a suitable catagory label and title. Now preform an overlay operation and display the resulting map layer. GRASS 4.1 > r.patch pipeline2,landcov output=pipe_final GRASS 4.1 > d.rast pipe_final You should be able to see the pipeline in the context of surrounding land uses. However, the pink color of the pipe does not show up very clearly. We can change the color of the pipe interactively. GRASS 4.1 > d.colormode mode=float (devote all colors to map) GRASS 4.1 > d.earse (clear graphics window) GRASS 4.1 > drast pipe_final (display layer to change) GRASS 4.1 > d.colors pipe_final (invoke color editor) Responses: Press d until the arrow at the top left is next to the category 9 - pipeline. Press the r key to reduce the red component of the color. You should be able to see this change on screen as you do it. Do the same for g (green) and b (blue). When you think it is sufficiently dark and visible color, press c to save the new color palette. A message should breifly appear on the screen indicating that the new palette has been saved. Finally press Q to quit, and go back into the fixed colormode by entering GRASS 4.1 > d.colormode mode=fixed Step 10: Display the final map using a perspective view. GRASS is able to generate psuedo-3D orthographic perspective displays. This provides an alternative mechanism for the visualization of surface data (such as a digital elevation model). Furthermore, the ability to "drape" a second map layer onto this surface presents an effective means of illustrating the relationships between the two map layers. GRASS 4.1 > d.3d Responses: Enter raster file to be displayed > pipe_final Enter raster file to be used for elevation > topo Enter name of 3-D veiwing options to be used > As a first attempt, the default image that has been drawn on the graphics screen is there to help you adjust the viewing parameters. You can improve the presentation of the image considerably by filling out the 'form' as follows. Eye Position: Run Y/N y 310000 <-Northing Erase Color black 444000 <-Easting Vertical Exag. 4 2500.0 <-Height Field of view 68.00 Lines only Y/N n Center of view: Line color brown 312000 <-Northing Line frequency 5 445800 <-Easting Resolution 50.00 0.00 <-Height Plot zero elvn n Box color none Average elevs. y Try changing these values to create alternative views of the landscape. Try to emphasise the relationship between the pipline route, topography and landuse. When you have finished, replace the entry 'Run Y/N' (at the top-right of the screen) with n before pressing . If you wish to save the viewing parameters, type in a suitable name and press . Discussion: Many GIS contain sophisticated graphical output capabilities of this sort. The results from an analysis undertaken using a GIS often need to be 'sold' to many other interested parties, and so high quality output that provides simple interpretation to the non-expert is often considered to be a valuable facility. Consolidation exercises (1) Generate a perspective display using d.3d that shows only the pipeline and the cost-distance surface (ie, without the land use information). Can you see why the pipe was routed the way it was? Is it the 'best' route? Would the route be different if calculated from plant to source? (2) The weightings used in Step 4 and Step 5 of the training session are fairly arbitrary. If different weightings were to be chosen you would expect the pipeline to take a different route. Returning to Step 5 (on page C-9) re-calculate the landcost map layer using the following alternative weightings (use the same file names, but prefixed with 'alt_'. Industry assign a value of 8 Residential assign a value of 8 Quarry assign a value of 1000 Woodland assign a value of 4 Arable assign a value of 16 Pasture assign a value of 1 Scrub assign a value of 1 Water assign a value of 1000 Repeat Step 6, Step 7, and Step 8, but remembering to use alt_cost and alt_pipeline. Discuss some ways in which the alternative routes, shown in pipeline and alt_pipeline can be compared. Briefly, not the differences in the routes taken, relating this to the land cover and the cost-distance surfaces over which they pass. Tidying up: The raster map layers that you have generated during this Training Session will not be needed in future sessions. If you wish, to conserve disk space, you have the option of removing these map layers, but do this with care. First type the command g.list to remind yourself of all the map layers now held on the disk. Remember, you are not deleting any of the files in mapset PERMANENT, but only the ones you have created in your own mapset. To close down any of the graphics windows you have opened during this session. GRASS 4.1 > d.mon stop=x0 Finally, exit GRASS. GRASS 4.1 > exit Responses: Shall the mapset be saved? > y (only select n if you do not wish to keep any files in the mapset) Do you wish to selectively remove data files? > y Please select type of file to be removed. 1. Raster files > 1 (give names of any files you wish to delete) > > GOOD BYE from GRASS! D: Site selection ================================================================================ Introduction Since you are now more familiar with GRASS, the instructions in this training session are again less comprehensive than those that have been provided earlier. A common application for GIS is in the identification of land parcels that satisfy multiple selection criteria. In this training session you will explore the ability of GRASS to preform this role, using a hypothetical site selection exercise as described below. A leisure center comprised of indoor facilities, playing fields, and car park is to be located on a site which satisfies the following criteria specified by the development agency. It must be :- 1. Within 500m of Shepshed, to provide easy access for the local community; 2. within 450m of the motorway or A and B class roads, for easy access by car; 3. on slopes of less than 2 degrees, in order to avoid high landscaping costs. 4. on land of Agricultural Grade III, to reduce site aquisition costs; 5. at least 2.5 hectares in area, in order to provide sufficient land for all the facilities of an integrated center; 6. out of view of as many populated areas as possible. Starting GRASS If you have not left GRASS since the previous Training Session, you can ignore this page, otherwise, you must first start GRASS and tell it which dataset you are going to work with. To start GRASS, type grass4.1 and press the return key. You should be presented with a screen similar to the following: grass@congo[2]> grass4.1 ----------------------------(clear)-------------------------------- GRASS 4.1 LOCATION: This is the name of an available geographic location. -spearfish- is the sample data base for which all tutorials are written. MAPSET: Every GRASS session runs under the name of a MAPSET. Associated with each MAPSET is a rectangular COORDINATE REGION and a list of any new maps created. DATABASE: This is the unix directory containing the geographic databases The REGION defaults to the entire area of the chosen LOCATION. You may change it later with the command: g.region - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOCATION: tutor_________ (enter list for a list of locations) MAPSET: PERMANENT_____ (or mapsets within a location) DATABASE: /usr/gis__________________________________________ AFTER COMPLETING ALL ANSWERS, HIT TO CONTINUE (OR TO CANCEL) You now need to tell GRASS which data set you will be working with. The opening screen prompts you to identify the LOCATION, MAPSET, and DATABASE. If you have just completed the previous Training Session, you will probably find that these are already set up correctly. If this is the case, just press to start. LOCATION stores the geographical extent of the study region. For this session you will continue to work with the Leicestershire data. Press the Return key a few times until the cursor is positioned on the LOCATION line. Then type leics to identify the study region. Don't worry about overtyping the default setting - just make sure that any remaining letters are deleted by pressing the space bar. Press Return to move the cursor onto the MAPSET line. Remember, a mapset in GRASS is your own set of files that reside in the current LOCATION. You should be using your user name. This will ensure that any files you alter or create will not be inadvertantly changed by anyone else. The DATABASE should not need changing on this occasion. It is simply the UNIX directory in which the sample locations and mapsets are stored. When you have set the correct options, press to start. ----------------------------(clear)-------------------------------- Welcome to GRASS 4.1 (Spring 1993) Update package 2 Geographic Resources Analysis Support System (GRASS) is a Trademark of U.S. Army Construction Engineering Research Laboratories (USACERL) New releases of GRASS are coordinated and produced by the Office of GRASS Integration (OGI) located at USACERL. This version running thru the C Shell (/bin/csh) Help is available with the command: g.help When ready to quit enter: exit Mapset in Location To use graphics in GRASS, you must open at least one graphics window. Therefore before you start this session, type: GRASS 4.1 > d.mon start=x0 select=x0 D: Site selection Aims: This training session aims to extend you knowledge of several fundimentsl GIS operations, and in particular to illustrate the very common technique of 'sieve mapping' for site selection activities, as well as a 'line of sight' function, and the conversion of raster data to vector data, with an introduction to the use of vector functions and concepts. Objectives: On completion of this training session you will be able to:- Prepare map layers that meet relevant selection criteria Construct a composite map layer by sieve mapping Identify and measure contiguous site areas Transfer from raster to vector polygons, creating a vector map Carryout a 'line of sight' investigation Improve the presentation of the 'results' map Display 2 GRASS commands to be used g.region r.reclass r.buffer d.rast r.slope.aspect r.mapcalc r.clump r.report r.poly v.support d.erase d.vect v.area r.los g.remove Display 3 Map layers to be used landcov Land cover / land use popln Population data roads Roads topo Digitial Elevation Model urban Urban areas Step 1: Preparing map layers to satisfy the location criteria This step involves the creation of a number of additional map layers that are derived from those listed in Display 2, each of which satisfies one specific selection criterion. Before we do this, we need to change the region so that we are looking at the NW quadrant of the study area. GRASS 4.1 > g.region Responses > 1 (modify current region directly) North 322000 West 444000 East 449750 South 316000 grid resolution East-West 50 North-South 50 This should give you 120 rows and 115 columns. GRASS 4.1 > d.frame -e Map layers satisfying the first two criteria can be generated using the buffering operation first used in Training Session B. Just to remind you, this is a common GIS operation that creates zones of 'distance' around point, line, or area features. Criterion 1. Identifying the area within 500m of Shepshed You need to derive a map layer that indicates which cells fall within 500m distance of Shepshed. This is a two stage process. First you will generate a map layer using r.buffer which codes Shepshed and the 500m buffer around it. Secondly, you will use the r.reclass so that only the buffer is coded as 1 (or TRUE) whilst inside and outside this buffer are coded as 0 (FALSE). GRASS 4.1 > r.buffer urban distances=500 output=urbanbuf GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > urbanbuf Enter name of NEW RECLASSIFIED map > urbanbuf2 Reclassify the layers as follows. OLD NEW distances from these locations 1 0 500m 2 1 TITLE: 500 meter buffer around urban areas Criterion 2. Identifying areas within 450m of the motorway and A/B class roads This buffering operation is very similar to that used to generate the urbanbuf2 layer. However, in this case it is necessary first to create a map layer containing just the motorway and A/B class roads before using the r.buffer command. GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > roads Enter name of NEW RECLASSIFIED map > main_roads All categories should be zero except the following. OLD NEW Motorway 3 1 'A' roads 9 1 'B' roads 11 1 TITLE: Main Roads (Motorway, A and B roads) Having extracted the motorway and A/B roads you can now proceed to buffer around the roads included in the new layer. GRASS 4.1 > r.buffer main_roads distance=450 output=main_roadbuf GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > main_roadbuf Enter name of NEW RECLASSIFIED map > main_roadbuf2 Reclassify the layer as follows and add a suitable title and category description. OLD NEW distances from these locations 1 0 450m 2 1 Use d.rast to check the results of this operation. Criterion 3. Identifying slopes less than 2 degrees You will begin this operation by creating a new map layer containing slope measurements based on the relief data held in the topo file. This is the same command that was used on page C-8, but this time slope wll be calculated rather than aspect (slope direction). GRASS 4.1 > r.slope.aspect elevation=topo slope=slope1 Now use r.reclass to create a new map layer in which cells with slopes of 2 degrees or less are coded with the value 1, and the remainder are coded with the value 0. GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > slope1 Enter name of NEW RECLASSIFIED map > flat All categories should be zero except the following. OLD NEW 0 degrees 1 1 1 degree 2 1 2 degrees 3 1 (all other categories should be zero) TITLE: Flat areas (slope <= 2 degrees) Again, display the results of this operation using d.rast. User Notes:................................................................ Criterion 4. Identifying Grade III Land This stage is relatively simple. It involves one step only, namely using r.reclass to select areas of agricultural Grade III land. The only problem is that you do not have a map layer that shows agricultural land grades! Therefore, take Grade III land to be equivalent to pasture and scrub land. GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > landcov Enter name of NEW RECLASSIFIED map > gradeIII All categories should be zero except the following. OLD NEW Pasture 6 1 Scrub 7 1 TITLE: Grade III agricultural land Again, display the results of this operation using d.rast. Discussion: When using GIS, it is often the case that the data wished for are not available. It is common for 'surrogate' data to be used instead. Although this can often provide a good indication of how the actual data would behave, conclusions drawn must be treated with some caution. User Notes:................................................................ Step 2: Constructing a composite map layer by "sieve mapping". You will now combine the four map layers that you have just prepared to produce a single map layer that indicates those cells in which all four selection criteria are met. The process of sieve mapping is most easily achieved by multiplying map layers together. Refer back to Training Session B, on page B-8 if you have forgotten why this is so. GRASS 4.1 > r.mapcalc Responses: mapcalc> sites = urbanbuf2 * main_roadbuf2 * flat * gradeIII mapcalc> User Notes:................................................................ Step 3: Criterion 5. Identify sites that are larger than 2.5 hectares. You will now process your sites layer to further select those sites that also satisfy the fifth criterion, ie that of a minimum areal extent. The first stage involves identifying the cells that form contiguous blocks of equal attribute code value. GRASS 4.1 > r.clump sites output=sites2 GRASS 4.1 > d.rast sites2 Discussion: The cells forming each contiguous group have now been allocated a new individual code value. You are now in a position to be able to calculate the area of each contiguous block. GRASS 4.1 > r.report sites2 units=hect Note down the categories for which the areas are larger than 2.5 hectares (excluding background). At this stage in the analysis you will use r.reclass to select the potential sites for leisure center that satisfy all of the original selection criteria. Each potential site needs to be classified uniquely. GRASS 4.1 > r.reclass Responses: Enter name of data layer to be reclassified > sites2 Enter name of NEW RECLASSIFIED map > sites3 You should now reclassify the new map layer so that those categories that you identified as greater than 2.5 hectares, are classified as uniquely, as 1, 2, 3 and so on. Whilst all those areas less than 2.5 hectares are classified as 0. Give the map layer a suitable title. Give each category a suitable label, for example "Potential site 1", "Potential site 2", and so on. Step 4: Criterion 6. Selecting the least obtrusive site. The finished leisure center is going to consist of car parking facilities, various palying fields, and a 10m high sports hall. This sports hall is going to be quite large, so in order for the development agency to win its planning permission, it must be done as unobtrusive as possible, that is it must be out of the view of as many people as possible. The sports hall will be built at the very center of the site. In order to carry out an evaluation of the visual impact, you need to know the grid coordinates of the center of each potential site. The potential sites are destinct spatial entites, having an area, a perimeter and a location, however this information is not explicitly recorded in the raster data structure of the sites3 map layer. GRASS can support the vector data structure, as well as raster. The vector data structure stores point, line (or arc) and area (or polygon) data explicitly, along with associated attributes. Within GRASS it is possible to convert from raster to the vector data structure (and vice versa), then with GRASS's vector functionality it is possible to find the coordinates of the center of each potential sites' polygon exactly. GRASS 4.1 > r.poly input=sites3 output=vect_sites The r.poly command converts the raster input map sites3, to the vector output map vect_sites. In its current state, vect_sites contains only nodes and lines, not the specific topological information that is needed. Therefore you need to use v.support to build the topology of vect_sites. (The prefix v. signifies a vector function). GRASS 4.1 > v.support map=vect_sites option=build Now display your new vector map layer. GRASS 4.1 > d.erase GRASS 4.1 > d.vect vect_sites You are now in a position to find the coordinates of the center of your potential sites polygons. Your new vect_sites vector map layer can now be queried, in order to find each polygons center. The v.area command enables you to find the polygon centroid coordinates, as well as each polygons area and perimeter. GRASS 4.1 > v.area map=vect_sites Responses: Buttons: Left: get area/perimeter Middle: quit this Right: get area/perimeter Point with the mouse to the inside of one of the polygons, click on the left or right button. Several lines of information should appear, including the perimeter and area. The last line gives you the coordinates of the centroid. Note these down in the Display 5 below, but be warned, the y axis coordinate (N) is given before the x axis coordinate (E). Find the centroid coordinates for all the polygons, then press the middle button on the mouse to quit. You may have forgotten the labels you used for each individual site, you need this information to complete Display 5. If this is the case use d.what.rast as follows: GRASS 4.1 > d.what.rast map=sites3 Click the left button when the mouse is pointing at a polygon. This will give you the site category label. Display 5 Potential x axis (E) y axis (N) Site No. coordinate coordinate ... .......... .......... ... .......... .......... ... .......... .......... Discussion: The use of the d.what.rast command shows that even when no raster map is being displayed you can recover information from the raster maps in the database. This shows how a geographical information system can retrieve information from the database using a geographical location coordinate, the map layer it comes from need not be displayed along with the query. As the final stage of selection you need to find the site that will be the least obtrusive. For this you need to use the command r.los (line of sight), this allows you to produce a new map layer of areas that are in view of the 10m high sports hall, within a 1km radius. The r.los command needs the elevation data contained in the raster map topo in order to carryout its function. You have recorded, in Display 5, the x and y coordinates of each sites center, where the sports hall could potentially be built, substitute the coordinates for potential site 1 into the following command line - instead of x and y. GRASS 4.1 > r.los input=topo output=MASK coordinate=x,y obs_elev=10 max_dist=1000 Discussion: The output map you have just produced is called MASK. You last used MASK on pages B-15 and B-16 where it was created by the r.mask command. However a MASK can be created by raster functions, in just the same way as ordinary rasters. Using the d.rast command now, you can see that the MASK is in operation. GRASS 4.1 > d.rast popln Displayed now are all the populated cells that could view the proposed sports hall. Use r.report to find the number of cells and record this in Display 6. Then remove the MASK. GRASS 4.1 > r.report popln units=cell GRASS 4.1 > g.remove rast=MASK In order to complete Display 6 you must now repeat the above procedure for each potential site, substituting their centroid coordinates from Display 5 in the r.los command line. (Don't forget to remove MASK each time). Display 6 Potential Number of Populated Site No. cells within view. .... ................. .... ................. .... ................. Now finally use r.reclass on the sites3 raster to produce a new raster called site, that has only the site within the lowest number of cells within view, classified as 1. Give the map a suitable title and category label. Step 5: The presentation of the final map. Now that you have identified the best potential site for the leisure center it would be useful to be able to locate it relative to known landmarks, such as the road network and urban area. This can be done by overlaying your final site map onto any of the map layers in your mapset. If you need to, use g.list to remind you what map layers you have available. Then use d.rast -o to produce a graphical overlay of a number of those images. Finally, graphically overlay site onto your composite map. There are many other commands that allow you to carryout alterations to the maps final appearance, for example d.his, d.grid, d.label, d.scale, and so on. These can be run ineractively, or with command lines, try experimenting with these commands. Use g.manual's list to find other display commands. Try using vector layers produced by r.poly (for areas), r.line (for lines), and r.contour (to produce a vector contour map from the topo raster map) - remembering to build their topology with v.support first. In conclusion experiment, GRASS contains many more commands than have been used in this beginner's tutorial. Consolidation exercises (1) How many parcels of land, of any size, were there after the sieve mapping exercise on page D-9? How can you obtain this information? (2) Examining your final map, does the potential site appear appropriate for the leisure center? How might you improve your selection process? (3) It can be instructive to repeat Step 2 on page D-9 using the addition operation in r.mapcalc rather than the multiplication operator. Try this and display the results. What is the significance of the various cell categories? (4) The Boolean logic used in sieve mapping is often criticised as it does not allow any grading of suitability other than acceptable/not acceptable. One consequence of this is that small changes in the selection criteria can lead to remarkably different 'results'. Try the site selection procedure with a 1000m buffer zone around the urban areas. How many potential sites are found now? Tidying up: Now that you have finished all four Training Sessions, you will not need any of the new raster map layers that you have created. If you wish to conserve disk space, you have the option of removing these map layers. First type the command g.list to remind yourself of all the map layers now held on the disk. Remember, you are not deleting any of the files in mapset PERMANENT, but only the ones you have created in your own mapset. To close down any of the graphics windows you have opened during this session. GRASS 4.1 > d.mon stop=x0 Finally, exit GRASS. GRASS 4.1 > exit Responses: Shall the mapset be saved? > y (only select n if you do not wish to keep any files in the mapset) Do you wish to selectively remove data files? > y Please select type of file to be removed. 1. Raster files > 1 (give names of any files you wish to delete) > > GOOD BYE from GRASS!