Handbook for Geographic Information Systems on Dominica

Note: This handbook was prepared to follow up a GIS training workshop provided for the Government of Dominica. The data files and color plates refered to in this document are not available in the online version of the document.

Contents

• Background
• Workshop on GIS and Utilities
• Getting Started
• Customizing a Palette
• Editing the Map
• Getting a Good Look at the Elevation Map.
• Surface Analysis
• Context and Distance Operators with Targets
• An Engineering Problem
• Unique Solution
• Points to Remember in Evaluating the Models
• SPANS and IDRISI: Special Workshop

• Hints for Using SPANS
• SPANS file maintenance
• Hints for using IDRISI
• Moving Files from SPANS to IDRISI
• Moving Files from IDRISI back to SPANS
• Conclusions from the Special Workshop

• References
• Colour Illustrations

Background

The Commonwealth of Dominica maintains modern infrastructure in difficult conditions. All the permanent settlements of the country are on paved roads, and the majority of the population has access to safe drinking water, electricity, and modern telephones. These services are maintained despite frequent landslides and recurrent hurricanes, but they are expensive. The landscape itself adds to the costs, because the population is scattered among plateaux and valleys separated by sheer cliffs. Given these conditions, careful planning is needed to minimize the costs of new installations and emergency repairs.

The Physical Planning Unit (PPU) of the Ministry of Commerce has a modern Geographic Information System (GIS), recently installed under the auspices of Canadian International Development Assistance (CIDA). PPU has been working to add to this digital archive and to use the system for planning at the country and district level. As it happens, a research center at Springfield, Dominica had a similar system, although it did not have a full database. The Geographic Information System-Environmental Planning project (GIS-EP) is funded by USAID and administered by OAS to support environmental planning and to foster the use of modern decision-support technology in the country. The GIS-EP project has supported both the PPU and Springfield GIS with equipment and technical support.

This handbook is a product of three weeks of workshops organized by the GIS-EP on Dominica in October, 1997

• The first two weeks were a Workshop on GIS and Utilities held at Springfield. Participants were engineers from Dominican utility companies, and they worked through a number of GIS operations including some modeling.
• The third week was a Special Workshop to help PPU integrate the GIS software which they had been using with a new software package, and to optimize their system.

The purpose of these workshops was to aid the Dominicans in making the most of their new systems. This required stimulating interest without arousing false hopes. To stimulate interest, it is necessary to get past the graphics screens and colour printers and advance to the mathematical modeling tools that give GIS its true power. The workshops explored some spatial models, and the handbook that follows develops the models further.

On the other hand, false hopes are a danger for novices. There are dozens of GIS software packages on the market, and some of them are especially designed to support facilities maintenance of public utilities. It is easy to watch a prepared GIS demonstration by a team of traveling salesmen without seeing the work and expense that is needed to keep one of these systems going on a day-to-day basis. There are now two general-purpose GIS on Dominica. If some utility on Dominica decides to get its own GIS, it may get a system that is tailored to its special needs, but no organization should miss the opportunity for collaboration and mutual support.

The workshops had a mixture of difficulties and solutions that were typical for this technology. As a result of this experience, the Dominican participants will be better able to judge and use GIS for their particular needs.

Workshop on GIS and Utilities

Following are a series of exercises based on a workshop at Springfield, Dominica in October 1997. The workshop was sponsored by the GIS-EP, with assistance from PPU. The participants were engineers from Dominican utility companies. During the workshop, we identified the following problem for analysis:

• The northwest of the island (See Figure 1, page 35) is showing the most rapid economic development, so the problem we would choose for the workshop should be situated in that area.
• The large estate at Pointe Ronde, on the south end of Prince Rupert Bay, is in the path of development although it has been fallow for many years. There have been repeated rumours of large capital schemes for the property, and there could be a sudden increase in demand for utility services in that direction.
• The terrain in the area is more level and open than is typical for Dominica, and presents little barrier to the expansion of the electrical or communications grids.
• The water supply in the area raises some interesting questions. The intake for the district is relatively low, at 360 feet of elevation on the Picard River. The water line in the direction of Pointe Ronde is only 10 cm (4 inches) in diameter, and it stops at Petite Baie.
• Current water use from the Picard River is only 1.4 of 12 million gallons per day (MGD) available. The immediate constraint is treatment, which could easily be upgraded.
• Even given that there is sufficient water, there might be head losses due to friction in the long narrow pipeline which would make it difficult to supply Pointe Ronde. This was a worthy problem for GIS modeling.

This handbook started as a memory aid for the workshop participants who might like to return to Springfield and repeat the exercises on their own. As I wrote, I tried to make the instructions readable for others as well, similar to the Tutorial Exercises provided with IDRISI. I assume that the user is familiar with Dominica and with computers using Windows 95 or Windows NT. There should also be some familiarity with the IDRISI software for GIS (Geographic Information Systems), although we will start with elementary operations. A paper map of the island at a scale of 1:50,000 or larger would be helpful. The digital elevation map and the vector files are to be distributed along with this document.

File names, directory names, and file formats will be in bold. A name for an IDRISI map will be both bold and with an initial Capital. As new concepts are introduced, they will get bold ALLCAPS. Commands to be executed inside a software package will be presented as a sequence of menu names with initial capitals: File/Import/Vector.

Getting Started

1. Establish a new directory on your hard disk, someplace where you have at least 90 to 100 Mb of space. An empty Zip disk will be fine. Copy the ten files DomDEM.img/doc, Coast.vec/dvh, Streams.vec/dvh, Roads.vec/dvh, and Houses.vec/dvh to the new directory. You may have to uncompress the elevation map using PKunZip or WinZip. Open IDRISI and set the environment path to your data directory using Environment/Environ.
   
   DISPLAY Click on Display/Display Launcher to see the Digital Elevation Model of Dominica (DEM) entitled DomDEM. Use autoscaling and a palette that gives smooth gradations, such as grey256, quant256, or ndvi256. This image will look a little strange with these colours, but we will be making a custom palette later.

DomDEM is a raster image of elevations on Dominica exported from a SPANS map at quadtree 13. The raster has a nominal resolution of 7.933 meters per cell. The SPANS map, in turn, was the product of a Triangulated Irregular Network with 600,770 nodes, produced from a list of point elevations compiled by Darrell Edgett for CIDA in 1996.

If we figure the area of Dominica as 751 square kilometers, division shows that the original point file averaged eight points per hectare, or one point per 1250 m². So, we know right away that most of the cells on this raster map, averaging one per 63 m², were created by interpolation.

2. WINDOW The window button is the icon of a blue rectangle. Use it to window in on the northwest quarter of the island, including at least the Cabrits Peninsula, Prince Rupert Bay, and the peak of Morne Diablotin. This should center your image over the Picard River valley. Use the Cursor Query Operation to find the coordinates for the left, right, top and bottom of this window. I got columns 17 to 1325 and rows 655 to 2050.
   
   Click on the menu for Reformat/Window and create a new image clipped to the window you just chose from Begging. I used columns 50 to 1249 and rows 700 to 1999 and called the new image Rupert. The result displays automatically at the end of the processing. Examine it to make sure that is includes the whole watershed for the Picard River, Portsmouth, and Pointe Ronde.
   
   DOCUMENT Call up File/Document to look at the parameters of the image Rupert. My choice of window gave me 1200 columns by 1300 rows. The minimum and maximum X and Y were updated automatically for the new image. The values of the image range from 0 (zero, naught) to 4700. That is the same range of values as for the DEM of the whole island, because we included a bit of sea as well as the peak of Morne Diablotin, which is the highest spot on the island. The elevations are in feet, at present, while the reference units are meters, to accord with the UTM projection of the original map.

It is time for a little FILE MAINTENANCE. Shrink IDRISI and open Windows Explorer to look at your data directory. The file pair DomDEM.img/doc totals nearly 42 Mb, and Rupert.img/doc is another 3 Mb. Your directory is getting crowded already. You will not need DomDEM again for these exercises, so make sure that you have a copy safe somewhere else, and erase it.

3. You have two ways to get rid of DomDEM. If you do it in Windows Explorer, you must remember to erase both DomDEM.img and DomDEM.doc. Neither file is usable without the other, and you will only confuse yourself if you leave just one. (You might produce an error message, but you will not crash IDRISI that way.) It is simplest for you to go back to IDRISI, go to File/File Maintenance, and erase the pair with one command there.

You will need to do file maintenance often. If you try to run one of the exercises below without enough storage space, the command will proceed until the disk is full. Then you will have to do file maintenance anyway to remove the temporary files that were not finished.

Customizing a Palette

You probably are not comfortable looking at the elevation map with the palette you have been using. Most people need more intuitive colours, so here is how to get them.

4. PALETTE Click on the button for Palette Workshop. This is an icon on the menu bar showing diagonal stripes of color. A window appears which shows 15 black squares and some choice buttons.
   
   Click on the black square numbered "0" and then move down to the blinking button near the word "Red." Drag that button to the right slightly so that the red value is increased to around 100. Now go to the buttons for Green and Blue and drag them to around 225 and 255, respectively. The box number 0 should change colour to a pale blue-green. This is the colour you will be using to indicate the sea on your maps, so adjust it as you please.
   
   Next, click on the black square numbered "1" and make it into a colour appropriate for low ground near the shore. I favour a pale tan.
   
   Next to the black boxes are two arrowheads, pointing up and down. Click on the down arrow and change the numbers for the black boxes until you reach number 255. Using the colour buttons, change box number 255 to a deep forest green. This is the colour that will be painted on the top of Morne Diablotin.
   
   Now go to the area in this window called "palette blend function." Set the lower anchor at 1 and the upper anchor at 255 and click on the Blend button. Suddenly all the boxes that were black get new colours which are linear interpolations of the two anchor colours.
   
   On the window for Palette Workshop, click on File/Save As and pick a name for your new palette. Save it in your working directory as an *.smp file. Now open the Display Launcher and display Rupert using autoscale and your "user-defined" palette. Overlay the vector Coast to see the exact coastline.

There are three problems:

• The new palette does not bring the landscape into sharp focus. It blends so gradually from sea to mountaintop that you get only a vague impression of the shape of the hills.
• There is an area near the top of Morne Diablotin that is just plain wrong. Window in around column 1000, row 1250 to see it. It seems as if a few points in the original data set were enregistered as 3000 feet too short. The result is a sort of crater in the digital mountainside.
• There are blue areas within the coastline, in areas that should be land. The most obvious area is the swamp behind the Cabrits Peninsula.

First, we will finish editing the palette interactively, using simultaneous windows. Then we will address the next two problems in turn.

5. Open the Palette Workshop again without closing the map you have displayed. Go through File/Open to get the palette you just created. (Remember that when IDRISI gives you a text box to enter a name, you can double click on it to get a list of names available for files that exist already.)
   
   Now you can modify the palette to create sharper contrasts. I suggest that you break the range 1 to 255 into four or five subranges, with a blending from the low to the high values within each of them. For instance, leave colour 1 as it is, modify colour 50 to be a bright, deep green, and then blend from 1 to 50. Does this make sense for the coastal scrubland?
   
   To see the results, save the palette on the hard disk again but do not close the Palette Workshop. Go to the Composer window for Rupert, which should still be on-screen. Click on Properties, double click on the text box to get the list of palettes available, and pick the same palette again. IDRISI reads the new palette file which you just saved, and the colours on-screen change.

Modify your palette until you have ranges of colours that suit you. Then, to accentuate the image even more, you can add contrasting lines. For instance, try making colour 51 bright red. When you redisplay, the new red forms a contour line. Figure 2, page 36 shows a palette that accentuates the terrain.

Editing the Map

The crater in the side of Morne Diablotin is obviously wrong. It is also in the watershed of the Picard River, and we do not want it to affect our analysis. On the other hand, it is in a heavily forested, unsettled area, so that we can guess the approximate shape of the ground there without worrying too much about precision.

6. f) Window in around column 1000, row 1250 again. Make your window small enough that you can see the blockiness of the individual cells, but large enough that you can see the trends of the ground that is well outside of the crater. One hundred cells square is about right.
   
   The crater distorts the values for a short distance around it. Outside of the crater, we see the steep hillside, elevations around 3000 feet, with a well-defined ravine leading down towards the Picard River
   
   DIGITIZING Click on the icon that looks like a red and yellow wheel. This starts the module for On-Screen Digitizing. You are going to create some vector files to patch the hole in the mountain. Name the vector file Patch1, choose to make a Polygon, and choose a feature ID of 1. Do not forget to type in an explanatory title such as, "Outside: 1^st patch for the hole in Diablotin."
   
   When you click OK and move your cursor over the map, you see that you have a circle with crosshairs. Move to the upper-left corner of your window and click the left button on the mouse. A small dot appears on the map. Move to the lower-left corner and click again. The dot stretches to a line between the two points. Click on the lower-right and upper-right corners to make three sides of a square. Now click on the right button and a line shoots over to your first point and closes the square. Save this square by clicking on the icon that looks like a bent arrow.
   
   Now digitize another polygon, inside of the first by about ten cells on each side. Take extra care not to cut across the crater. Call this Patch2, value 1, with the title "Inside: 2^nd patch for the hole."
   
   Look at the elevations around the crater. The ravine appears undistorted up to 2800 feet and it resumes a plausible course above 3200 feet. I suggest that we plug the hole with values at 3000. So, digitize a polygon around the crater, Patch3, with ID equal to 3000. My polygon wound up looking like a cucumber. Save the file by clicking the icon of a bent arrow.
   
   You can check to see if you have saved all three patches by clicking on Composer/Add Layer as if you were going to add a vector to the image. The list of vectors to choose should include your three new creations.

What we are going to do is isolate a part of Rupert using Patch1 and cover the crater with Patch3. Then we are going to smooth over the seams a bit, cut out the part we want, and put it back over the hillside on Rupert.

7. MAKING A RASTER Go to Data Entry/Initial to create a blank map with the same spatial parameters as Rupert. Name it Maskpat1, with an initial value of 0. Title it, "Mask to take the area to be patched."
   
   Go to Reformat/Raster-Vector Conversion/Polyras to update your new raster, Maskpat1, with the area enclosed by Patch1.
   
   This new map does not look like much. It is all zeros except for an area of ones, corresponding to the area we want to edit. Maps like this are called BINARY MASKS, and they are powerful tools.
   
   MASKING OUT Go to Analysis/Mathematical Operators/Overlay and multiply Rupert by Maskpat1. Call the result Patchout, with the title "area of RUPERT to be doctored."

Have a close look at the result. Areas that were zeros on the mask are still zeros; areas that were ones are now the same value that they would be on RUPERT. So, multiplying a map by a binary mask is a sort of Boolean algebra. The mask says no to some areas and yes to others.

8. SMOOTH FILTER Update this image with Patch3. The crater disappears. Now go to the command Analysis/ Context Operators/Filter. Choose a mean-filter, 7 * 7 and call it Smooth. A mean-filter re-assigns each cell to the average value of the cells surrounding it. (You can read a full description if you click on the Help button) The effect is to smooth out the rough transitions, and we use it to smooth the edges of the patch over the crater.

Unfortunately the mean-filter also rounds the edges where the patch of good values comes up against the zeros of the background. So, using the operations you just learned, make a new initial map of zeros and update it with the polygon Patch2, and use the overlay multiplication to trim off the rounded outside edges of the square in Smooth.

9. COVER Now go back to Analysis/Mathematical Operators/Overlay. This time, choose "First covers Second except where zero" and put your leveled, smoothed and trimmed patch on top of the crater on the original Rupert. The repair is complete.

Getting a Good Look at the Elevation Map.

We have improved the palette and patched the hole in the mountain. The problem with the coastline requires more thought.

10. Window in on the south shore of the Cabrits Peninsula. From the cruise ship terminal to Portsmouth, the vector of the coastline looks right, but the colours of the raster seem to show massive flooding. Now do a cursor query, using the button on the menu bar that has the question mark. You will find that points offshore are registered as 0, while points inside the coastline have positive values even if they are blue. So the raster has the right information, but it looks wrong on-screen.

All of this is a result of the "autoscaling" option when you display. There are only 255 colours to allocate to 1400 elevation values in the raster Rupert. The lowest elevations on land get grouped with the sea itself, and the whole group is coloured blue.

After spending time creating a good palette for displaying the elevation data, you probably want to fix this problem in displaying the coastline. We will do that, but in order to make a more presentable picture we will have to distort and destroy some information.

11. RECLASS First, we will make a map that shows only the coastline. Click on Analysis/Database Query/Reclass. In the window that appears, choose the options "Image" and "User-defined reclass." Double click on the text box for "Input file" and select Rupert. For the output name, type "Bayside." For a title, enter "Binary mask of Prince Rupert Bay coastline." Now go to the right side of the window and assign a new value of 1 to all values from 1 to 4701. Notice that 4701 is one unit more than the maximum value in the map Rupert. Now click "OK." The status bar on the bottom of the screen shows the progress of the program, and a few moments later, a green and black map appears.
    
    Window in on a distinctive area, such as Cabrits, and overlay the image with the coastline vector. The vector seems to fit this raster much better than it fit Rupert. It is still not a perfect fit, and we will look into that problem when we evaluate the models later. The point for now is that this raster of the coastline was entirely derived from Rupert, so we know that the elevation map did have the information even if it was difficult to display.
    
    Go to File/Document and look at the parameters for Bayside. Look at the minimum and maximum values. This map has only two values: zero for the ocean and 1 for the land, so it is a binary mask. Obviously, we have destroyed a lot of information to make this simple outline from the elevation map. Look at the data type: "Byte." Close the Document window, shrink IDRISI and use Windows explorer to look at the size of the files in your data space. Bayside, the map of data type "Byte" is only half as big as Rupert.

The elevation map, Rupert, is of data type "integer." Each cell in a raster map of integers requires two bytes of disk space, but the integers can range over tens of thousands of values. When you created a map that only needed to range over two values, 0 and 1, the software chose a more compact data type for you. Rasters of bytes take only half as much disk space, but they can only range from 0 to 255. We will save Bayside and use it many times. Close Explorer and return to IDRISI.

12. MAP ALGEBRA We still want to find a way to display the elevations so that there is one colour for the ocean, and a series of colours for the different elevations on land. We know that that there are 4700 feet of altitude to be accommodated in 255 colours. We need to change the values of Rupert by some number so that the entire range is less than or equal to the number of colours available. Figure that 4700 feet divided by 255 colours equals 18.43, so that if we use one colour for each twenty feet we will be safe.
    
    Click on Analysis/Mathematical Operators/Scalar. This is a program that lets you change every cell in a raster by the same operation. Use Rupert as the input file and "Rupby20" as the output name. Choose to divide the map by 20. Enter a long, descriptive title such as, "Elevation map in feet divided by 20 to fit palette." Click "OK" and wait for the program to finish. When the map displays, click on the "Properties" button of Composer to change to the palette you created earlier.

Darn. It looks just the same. What happened?

13. Window on Cabrits again, overlay the coastline vector, and cursor query the blue area inside the coastline. You get decimal fractions: what was two feet high is now figured as 0.1, and so on. The numbers have changed, but once again the autoscaling feature of Display is grouping the low ground with the sea water.
    
    Open File/Document and look at the parameters for this newest file. Minimum value is zero, because you divided zero by 20. Maximum value is 235, which equals 4700 divided by 20. That is all as expected. The map is a different data type however, called "real." Close Document and shrink IDRISI again to use Windows Explorer to look at your hard drive. Rupby20.img is twice as big as the file it was derived from.

What is going on? Did you produce twice as much information?

No. The new map, with decimal fractions, has almost exactly the same information as the integer map. If you multiply 20 times Rupby20, you will restore the cell values almost as they were, with two differences. First, the cell values will still be in the data type "Real" with decimals. Second, you will have to think about rounding errors. In this case, since you divided by 20 you did not produce any irrational fractions, so there is no rounding error.

The Real data type takes four bytes for each cell value in the raster. It gives 7 significant figures; far more than we can justify for this data. Besides, right now we are intent on the appearance of the image, rather than precision to five decimal places. In that spirit, we are going to simplify the data and destroy information, but we will get a picture we can really study.

14. CONVERTING Click on Reformat/Convert. Choose Rupby20 as the input file. The window tells you that you are starting from data type Real and file type Binary. On the right half of the window, click the button to choose output data type Byte. On the bottom of the window, choose to do the conversion by truncation. Click "OK."
    
    Redisplay Rupby20 and investigate the values again. The truncation made each cell fall into the next whole category downward. That is, a cell that had 35 feet in Rupert had a value of 1.75 as a real number and it truncated to 1 as a byte value. A cell of 40 feet elevation is now classed as 2. That is all well and good, but the low areas are now truncated to zero. What to do?
    
    Go to Analysis/Mathematical Operators/Overlay. You have already seen how the command "Scalar" can add, subtract, multiply or divide an entire map by a number you choose. Overlay includes operations to add, subtract, multiply or divide one map by another map. Add Bayside to Rupby20, with the title "twenty-foot contour."

When you display this with your special palette, you finally see the shoreline. The elevations look just as detailed as ever, because you are using almost the same number of colours as before. So, the visual impact of the image is definitely improved, even though the data has been distorted. All of the above may seem like a lot of work for very little improvement, but there will be a reward.

15. Go to Display/Ortho and click on the help button. As you read the description of this command, you will see that Twentyup is just what you need as a "drape image."
    
    To use Ortho, you must adjust the maps a bit to match the resolution of your video monitor. For instance, Rupert has 1200 columns by 1300 rows, and can not be shown on a monitor with only 1024 by 768 pixels. Use cursor query on either Rupert or Twentyup to see how many columns you can trim off the west edge without cutting into Pointe Ronde, and how many columns you must lose from the mountains to the east. Then, to trim Rupert and Twentyup down to size, use Reformat/Window, choose batch window, and pick a prefix like "ort." Type in the rows and colunms you need. I chose:

• Upper-left column = 63
• Upper-left Row = 0
• Lower-right column = 1086
• Lower-right row = 1299

Now, use Ortho with Ortert as the surface image and Ortntyup as the drape image. Pick a resolution to suit your monitor and a palette to suit yourself. Run the command several times to experiment with different viewing angles, directions, and vertical factors.

Orthogonal views like this appeal to us because they bear some resemblance to what we see in real life. For instance, we see a high oblique view of Prince Rupert Bay when we fly from Guadeloupe to Canefield. You might want to look at orthogonal views to search for obvious data errors like the crater you just finished filling, or you might want to use the functions in Composer to annotate the image and prepare it for a presentation.

Aside from such applications, the orthogonal viewer is a dead end. It pleases our eyes, but it does not lead us to new knowledge. In the next few exercises, we will finally enter the field of GIS analysis. We will use powerful tools which show us the maps in ways we could not see without help.

16. So, save your best orthogonal view if you want to. Also, please save

• Bayside, the binary mask of the coast,
• Twentyup, the elevation map divided into twenty-foot intervals, and
• Rupert, the raw elevation map, with the crater patched up.

Make sure you have good explanatory titles for each of the maps above, and then use File/File Maintenance to erase all the others.

Surface Analysis

We look at the orthogonal image, and we can see that some areas are level, others are steep. How steep? It is impossible to be precise, since we have shifted the angle of view and the vertical exaggeration factor.

Abandoning orthographic views, we could display the image in plan view again. We could window in and use cursor query to get the elevation value of a point of interest and of its immediate neighbors. If we jotted the numbers down on a piece of paper, we could calculate the slope by figuring the difference in elevation between our point of interest and a cell next to it. But which neighboring cell should we use for comparison? We would take the lowest neighbor for going downhill, the highest neighbor for going uphill, but what if there are neighbors both above and below? Which way do we go? This could take a long time.

Fortunately, IDRISI has a command which will figure this out for you. It will search the entire image, compare each cell to its neighbors, and fit an inclined plane to that location. This is a "best fit" approach, because hillsides are not inclined planes: they are bumpy.

You can think of the inclined plane as being tangent to the surface at the center of the cell, slanted so that it passes as close as possible to the neighboring cells. The slant is specified by the slope and the direction of the slope, the "aspect."

17. Go to Analysis/Context Operators/Surface. Choose to calculate both slope and aspect on Rupert with output in degrees, and a conversion factor of 0.3048 meters per foot. (Remember Rupert’s elevations are in feet, but the grid coordinates are in meters) Call the outputs Rupslope and Rupaspct. Watch how quickly IDRISI calculates 1.5 million slopes and aspects for you.
    
    As always, we will use the Cabrits Peninsula to check our output. Display Rupaspct with the grey256 palette and window into Cabrits. Look at the southern half of the Peninsula. The image is a bit strange, but it does have some intuitive appeal as a picture of the surface, lighted from the upper left. The north half of the Peninsula is harder to read. Use Cursor Query and you will find that there are strange values. All the sea is marked "-1" and so is the swamp and the bottom of the valley between the hills. That number is a special value meaning "flat." The rest of the image is in azimuths, degrees measured clockwise from North. That is why the northeast side of the hills is dark (values 0 to 90 degrees), while just a bit around the corner, on the northwest, the hillside is at its brightest (values 270 to 360 degrees, the maximum). You can take the time to design a special palette for viewing the aspect map, if you like.
    
    Check the slope map, Rupslope, for the same area. Let IDRISI pick its own palette, Idrisi256, and overlay the coast line with the Qual16 line colours. The two hills look like two rings of fire. The sides are bright reds and yellows, while the hilltops are dark. The valley between the hills is dark, too, and so is the swamp.
    
    Now window out and look over the whole image. You may want to use the coast-line and the hydrology line files to get oriented, but then switch them off temporarily by clicking the check boxes on Composer. You see that the Picard River and the Espagnol River show as dark valleys with bright borders on either side. Some of the ridgelines have the same effect, and they could be mistaken for valleys.

Most people have never seen a map like this. It is confusing at first, with the hills and low ground looking the same, but it makes sense if you think about it. Sometimes a ridgeline will lead down to a promontory with cliffs on three sides, and sometimes a river will go over a precipice and form a waterfall. Those are the exceptions that prove the rule, and the rule is that valleys and ridge are the gentlest slopes. That is were our footpaths go. Think of it as a mathematical limit: valleys and ridges are formed by the intersection of opposing slopes, and the slope of an intersection is less than or equal to the lesser of the opposing planes.

18. HISTOGRAM Ask for a histogram (Analysis/Statistics/Histo) in graphical output. You get a histogram showing the frequency of each range of slopes. Unfortunately, there are a lot of slopes equal to zero, because our image includes Prince Rupert Bay. The average slope with all that flat water included is 16.9 degrees.
    
    Ask for a histogram again, but this time set the display minimum to 0.01 or some other low number slightly above zero. Now the histogram shows only the values for land. The average slope is 22 degrees. That is low for Dominica, but it is still much higher than most parts of the world. Look at the distribution graph: there is more terrain with slopes between 15 and 25 degrees than there is between 0 and 10 degrees, even with this image which includes the swamps of Cabrits and Indian River.

Context and Distance Operators with Targets

We have seen two commands in IDRISI that do calculations involving the immediate neighborhood of each cell. The first such command was Filter, which we used to smooth the patch on the crater hole. The second command was Surface, which we just used to create slope and aspect maps. These two commands have something in common: they process the whole map, impartially, by operating on each location in turn.

Now we are going to use some commands that we can aim at specific areas. With these commands, we create a "target image" and make the computer organize its calculations around the target.

This next operation may possibly take a long time to run on your computer. You may run it overnight, or skip it altogether.

19. Display Rupert or Rupaspct or Rupslope and overlay the hydrology lines to help you find your way. Window into the central Picard River. You need to include the cell at column 533, row 589 as well as at least a kilometer of river above that. The cell I mentioned is the approximate location of DOWASCO’s water intake.
    
    Go to Analysis/Context Operators/ Watrshed and read the help file to learn what a "target image" is. Now digitize a polygon around a portion of the Picard River, from the water intake to about a half kilometer upstream. Avoid Mary Gutter, the side stream that joins the river from the south just a little bit downstream of the waterworks. Convert the polygon to a binary mask using Reformat/Raster-Vector Conversion/Polyras as you did before. Name it Intake.
    
    Run Watrshed with your new binary mask as a target. The aspect image is Rupaspt, output map will be Picshed. The command will look around the target image, adding any cells facing toward the target in the aspect map. It may take over an hour to run, even at this resolution, and you probably will not like the result.
    
    When Picshed displays, overlay the hydrology lines, coloured by Qual16 so that you can interpret it. Notice that some side branches of the Picard River are included in the watershed, while others are not. There may even be a portion of the Espagnol River included by accident. These mistakes derive from the elevation map. They can be edited out, or you can re-run the command with new targets for the side streams. Watrshed is a powerful command, but it is only as good as the data you have to put in.

To get the best analysis of surface flow, you would have to edit this elevation file carefully. You could also filter the image in a variety of ways to make it more realistic. It might even be worth while to go back to SPANS and export the elevation map at a higher resolution, or take a chance on the slope and aspect generator in SPANS. The better the source, the more you can do with it.

20. t) Display the binary mask of the area, Bayside, and overlay the roads. Notice that the roads all look the same this way, even though some of them are no more than tracks through the mountains. Digitize a polygon around the roads that link Portsmouth, the water intake, and Pointe Ronde and turn the polygon into a binary mask. Include the side streets in the area of Picard estate and Petite Baie, but exclude Glanvillia. Name this mask Axis.
    
    Now convert the roads file into a raster named Allroads using Reformat/Raster-Vector Conversion/Lineras. You can look at Allroads and see that the road lines are now shown as strings of little squares. Check the maximum value of this new raster with File/Document, then Reclass Allroads to a binary mask, Roadall, where all the roads have value equal to 1 and background is zero.
    
    Do an overlay multiplication, Axis times Roadall, and name the result Access. (You may have to redisplay Access with autoscaling in order to see it.) Now the only roads that show are the ones inside Axis. Because most of the water line and local power lines run parallel to the roads, Access gives you a new target image of the utility corridors along the south side of Prince Rupert Bay.
    
    Run the commandAnalysis/Distance Operators/Buffer with Access as the feature image and a buffer distance of 50 meters. Let the road itself stay as value 1 and make the buffer equal to 2. Call the image Buffr50 and overlay the coastline to look at it. Fifty meters is about as far from the road as people can build without major expense for power poles, water lines and driveways, so this map shows the area that has "road access."
    
    Run the command Analysis/Statistics/Histo with buff50 as the input map, class width 1.0 and numeric output.. For my test, there were 2355 cells classed as Road and 26101 classed as Road Access. Since each cell is 7.933 meters on a side, the model has estimated something like 160 to 200 hectares of real estate with road frontage between Portsmouth and Pointe Ronde.

You can easily refine this model by digitizing more carefully, by using buffers of different widths, or by masking out the bay.

21. Rather than modeling one buffer width at a time, you can model all of them at once. Try the command Analysis/Distance Operators/Distance, using Access as the feature image again. Call the output Fromroad.
    
    When Fromroad displays, overlay the coast line and change the palette to Qual256. That palette gives contrasting colours to adjacent values, so you can see the bands of same-distance around the target set of roads.

Obviously, the distance rings going out to sea are irrelevant to utility development, so do a multiplication overlay with the binary mask Bayside. Name the result Fromrod2. The overlay operation cuts out the sea and makes a more tidy image. You could also make a binary mask to cut out the bands to the north and south which wrap around the target set of roads.

Window in on Fromrod2 and look at how the distance rings interact in the places between adjacent roads. There are areas where a new street could connect two roads and open more area for development. It looks easy on this map.

22. Dominica is not easy, however. Two roads separated by 100 meters of horizontal distance may be separated by 100 meters of vertical distance as well. Go to Analysis/Distance Operators/Cost and read the help files. This command allows you to combine distance calculations with "friction." Friction can include cost factors such as terrain, soil type, real estate prices, or whatever you wish.
    
    Terrain is the major factor on Dominica. Call up Analysis/Database Query/Reclass and use the slopes map to create a new map of land suitability classes, Landsuit.

• 0 to 5 degrees, Class 1
• 5 to 10 degrees, Class 2
• 10 to 20 degrees, Class 3
• 20 to 35 degrees, Class 4
• 35 to 90 degrees, Class 5

Use Analysis/Database Query/Overlay and multiply Landsuit by Bayside to cut out the sea, again. Call the result Landsui2, but stop to think: if steep terrain is Class 5 because it is difficult to build on, should the sea be Class zero? That implies that it is no effort to build on at all. Obviously not. Reclass the sea from zero to 100.

Now use Landsui2 with Access as the feature file for Analysis/Distance Operators/Cost. Use the Cost Push option to save time. When the program finishes, use an overlay with Bayside to mask away the high values out to sea. Display with the quant256 palette and overlay the coast and roads. Window into Pointe Ronde.

Notice how the high cost values form a "ridge" around column 325, row 660. The area has roads on three sides, but it is difficult to access because of the steep ravine on the Lamoins River. A similar ridge is visible on the high ground above Pointe Ronde. By contrast, the flat areas along the lower Picard River, the old Picard Estate, are easy access.

23. The friction surface only showed five classes. For a more subtle model, try calculating friction as a trigonometric function of the slope. Look at Analysis/Mathematical Operators/Image Calculator.

This calculator lets you do scalar operations, overlays, and a variety of transformations simultaneously. You can enter the name of an image just as easily as a number. (Note: The Calculator will only work on images with names in the form *.img and *.doc, the defaults for IDRISI.)

With Calculator, you can save many steps in any analysis. Try the formula

TANFRICT = 1 + 2*(TAN(RAD([RUSLOPE])))

That is, calculate the new image Tanfrict by converting Ruslope to radians, figuring the tangent of each slope, doubling it, and adding 1 to it. The friction will never be less than 1, on flat land, and it will range to infinity at sheer cliffs. Luckily, we know that the maximum slope registered is around 88 degrees, so the tangents will range up to 28 or so, and Tanfrict should not go above 57.

If you run Analysis/Distance Operators/Cost on this new friction model, the "ridgelines" between access roads are much sharper, where the slopes are steep. Keep Tanfrict.

This is as far was the workshop progressed before leaving the next day to visit the Picard River intake. The exercises presented below are a continuation of the same line of investigation, for the reader to pursue independently.

An Engineering Problem

Let us return to the problem of distributing water from the intake on the Picard River to the new houses that might be built at Pointe Ronde. Suppose that you work for Dowasco, the Dominican water utility. You find out, late on a Friday, that there are fifty new luxury homes being built at Pointe Ronde. Somehow, you acquire a plat of the house sites and digitize their locations as the vector point file Houses. You want to figure if there is any hope of supplying the new houses with water by simply extending the 4-inch water main you already have at Petite Baie.

We will use English Units, and the Hazen-Williams equation for head loss due to friction of liquid flow:

Where

Here are some relevant numbers gathered from discussions during the workshop, from textbooks, and made up, where necessary:

• The pipe system is ductile iron, with a Hazen-Williams coefficient of 120.
• The INTAKE on the Picard river feeds into a ten-inch pipe which runs down the access road toward Portsmouth.
• At the JUNCTION where that pipeline meets the coastal highway, there is a four-inch pipe which leads off to the south as far as Petite Baie.
• We choose to imagine that 2% of the flow turns south at the junction.
• The pipe network parallels the roads.
• Average use at the intake is 1.4 million gallons per day (MGD).
• Peak daily use may be figured as 180% of average daily use
• Peak hourly use may be figured as 150% of average use for the day.
• One foot of pressure head = 0.433 pounds per square inch (psi).
• Minimum pressure in a drinking water supply >= 2 psi.
• Ordinary service for residential areas is 20 to 40 psi.
• One gallon = 0.13376 cubic feet.
• Luxury homes with gardens and full appliances use up to 90 gallons per person per day, on a long term average.

There are three nodes in the pipe network of concern to us:

  We will also need some dimensions for standard pipes:

24. To start on this problem, look at the help file for Analysis/Distance Operators/Cost, and re-read the help file. We are going to use the Cost Grow algorithm this time, and force the program to measure its way along the road lines leading away from the junction.
    
    You need a friction Surface that forces the program to stay on the roads. To get this, you will need to reclass Access so that the non-road areas change from a value of zero to a value of minus one (-1). This is a special flag value that tells the algorithm that those areas are barriers.
    
    You also need to digitize a point at the junction (column 370, row 484), and make it into a source image by Reformat/Raster-Vector Conversion/Pointras. These preparations are all similar to the exercises in the previous pages.
    
    When you have the source image and the friction surface ready, run the Cost Grow program and name the result Roadist. The output is a raster with -1 in most places and distance values assigned to the cells with roads in them.

The distances on Roadist are expressed in cell units. You can use these cell distance units to model completely different scenarios. Each cell is 7.933 meters across, or 26.03 feet. That stays constant. What changes is the head loss per unit of distance, because friction head loss are a power function of flow rate, as shown in the Hazen-Williams formula above.

WORST CASE First, we will model to see if we can supply the new houses with adequate pressure during the hour of heaviest use, during the day of heaviest use.

• Average use at the intake is 1.4 MGD.
• Use on the peak day of the year would be at the rate of
  1.80 X 1.4 MGD = 2.52 MGD.

• Average hourly use on that day would be
  (2.52 X 10⁶ gal./day) X (1 day/ 24 hours) = 105,000 gal./hour

• Peak hourly use on that day would be
  1.5 X (105,000 gal./hour) = 157,500 gal./hour
  (157,500 gal./hour) X( 1 hour/3600 seconds) = 43.75 gal /second
  (29.17 gal./sec.) X (0.13376 cubic feet/gal.) = 5.852 cubic feet/sec. (cfs)

We will assume that all of that flow comes down the hill as far as the junction, and plug the numbers into the Hazen-Williams equation. Express the result as head loss per foot of run:

• Average velocity = Q/A = (5.852 cfs) / (0.545 ft.²) = 10.74 ft./sec.
• Head loss factor = 0.04326 feet lost / foot of run.

Display the map Roadist and check the value at the intake. It should be 237.1372 cell distances. The number is fractional because of diagonal movements along the line. To figure the head loss at the junction multiply,

(237.13 cells) X (26.06 ft./cell) X ( 0.04326 ft./ft.) = 267 feet.

But wait, the intake is at 366 feet, the junction is at 36 feet, so figure the net head,

(366 - 36) elevation difference - 267 feet, friction head loss = 63 feet.

So the junction would have 63 feet of head, or 27 psi in the worst case scenario. Things do not look good, but let us continue down the line to Petite Baie.

We will do this next stretch graphically. We assume that 2% of the flow turns south at the junction Once again we assume that all the flow goes to the end of the link. That may be unrealistic, but it is a conservative estimate which will over-estimate head losses rather than under-estimating them. The link is a 4-inch pipe all the way. We need one head loss factor for the whole link.

• Flow = 0.02 X 5.852 cfs = 0.11704 cfs
• Average velocity = Q/A = 0.11704 cfs / 0.0873 = 1.34066 ft./sec.
• Head loss factor = 0.002692 feet lost / foot of run.

25. Local head, at any point along the 4-inch pipe, equals the head at the junction, plus any head due to changes in elevation, minus the head loss due to distance times friction. Make sure you have plenty of space on your hard disk and then use Analysis/Mathematical Operations/Image Calculator to combine everything in one step. And, while you are doing that, do a multiplication overlay with Access to mask out the off-road areas. Call the new map Hloss4in:

[Hloss4in] = (63 + (36 - [Rupert]) - (0.002692 * 26.03 * [Roadist])) * [Access]

Display the result and window in on the coast road. Remember that the result is ONLY valid along the road that has this flow and this pipe diameter, the road from the junction to the pipe end at Petite Baie.

You have a variety of tricks to help you see this image clearly. You can switch palettes, you can overlay the line file for the coastline, the rivers, or the roads. You can even digitize around the coastal highway and make a binary mask to isolate the area where the values are valid. I use an overlay of the road lines, window in to a large scale, and do a cursor query cell-by-cell along the road.

It does not matter that this image looks strange. It is full of information just the same. The pressure head, during the worst case scenario, is 60 feet at the pipe end. That is better than we expected when we found that the pressure at the junction was only 63 feet. You could have figured that with a calculator and a ruler. Be careful not to dismiss this problem to quickly, however. Work back along the road from the pipe end. There is a small hill a quarter mile north of the Lamoins River, that rises to 120 feet. Pressure there goes negative, to -25 feet. Not only would those people be without water, there would be danger of contaminants being sucked in.

There is no point going further towards Pointe Ronde with the Worst Case Scenario. We are already discovering places where pressure would fail, on the pipeline that exists already. Some of the houses at the new development would be higher than this, so a simple gravity feed through a 4-inch pipe would not do.

Unique Solution

The problem area near Lamoines River , which was revealed by the model above, is served by storage tanks during peak demand. Our next model will show how to locate and supply a storage tank for Pointe Ronde.

We have shown that a four-inch pipe could not deliver water directly to the new development, by gravity, during peak demand. The rapid flow of water in other parts of the system would deplete the energy too much. On the other hand, remember that the water intake on the Picard River is at 360 feet elevation, and the highest house in the new development is at 200 feet. There is, indeed, a difference in energy level, and that energy gradient could drive some flow.

26. Now digitize the obvious route for a pipeline, from (column 270, row 636) at Petite Baie, along the coast road as far as (column 137, row 781) at Pointe Ronde, and up the track on the hillside to (column 184 row 759). Create a raster of that and create a buffer map of it where the buffer is 50 meters wide, the buffer and the road have a value of 5, and the off-road areas have a value of 1. Do not let any areas have a value of zero. Call it Buffat5.
    
    We are going to elaborate the friction map based on tangents of the slopes, Tanfrict. Reclass Bayside to Bayout so sea water is -1, and use Bayout to Overlay/Multiply Tanfrict and make the sea water a barrier. Call this Frict2. Divide Frict2 by Buffat5 to make Frict3 where friction factors are 1/5 as great near the road. (Now you see why no part of Buffat5 may be zero.) Dispose of Frict2, Frict3, Buffat5.
    
    Frict3 is a friction map that is weighted to give heavy penalties for building across slopes, forbids building across the sea, and factors in that it should cost about 1/5 as much to build along an existing alignment.
    
    Digitize the pipe end at Petite Baie, (column 270, row 636). Make it into a source feature with just the one cell as the target, called Pipend.
    
    Caution: This next command takes a lot of time: run it overnight. If you do not have time, you can find a copy of the result in the data base provided with this handbook.
    
    Use Cost Grow to create a cost surface centered on the pipe head at Petite Baie, with Frict3 as the friction surface. Call it Progress.

So now you have a complex model of cost factors and a map of the cost surface grown outward from the pipe end at Petite Baie. You can clearly see the corridor of low cost along the road to the south. The off-road areas are much higher cost, but only the sea is an absolute barrier.

27. You have a vector point file of the locations for the new houses. Make a raster map of them and create a buffer map with a buffer distance of 100 meters, where the buffer and the houses are all value 1. Reclass Rupert to Min250 so that areas below 250 feet are zero. Overlay the elevation mask with the 100 meter buffer to create an image of areas that can hold a storage tank for Pointe Ronde. Call it Localhi.

That creates a target image of the area on the hill above Pointe Ronde where you could site a storage tank to supply the new houses. The area is at least 50 feet higher than the highest house, but no more than 100 meters distant from the nearest part of the new development. You can see that the road goes around three sides of that hill. But it is hard to tell, just by looking, whether you should run the water line along the road or take a shortcut south across the Cario River.

28. Call up Analysis/Distance Operators/Pathway. Use Progress as the cost surface and Localhi as the target surface. Call the output Shortcut. The result is a raster map of the least-cost path from the pipe end at Petite Baie to the area where you could put a storage tank.
    
    Spend some time looking at this result. The first step might be to overlay the road alignment Obvious and see how the least-cost path cuts across the loops of the road twice before leaving the road corridor altogether and heading straight southwest. Next, you might want to digitize the pathway to produce a vector version of it that you can overlay on raster images like Rupert, Rupslope, or Progress.

There are two directions you can go from this model. First, you can modify the model and run it again, if it does not suit you. For instance, we modeled the roadway as being five times more economical than building cross-country. Perhaps that was not high enough. What would happen to the pathway if we figures the road as ten times easier? Perhaps you have some other engineering rule-of-thumb you would like to apply. Your experience and judgment are keys to a good model.

Second, you can use this pathway to predict the flow. There would be no flow during the Worst Case Scenario, but there is 110 feet of elevation difference to drive a gentle flow at night. Most of the assumptions you need are listed at the beginning of this section.

Points to Remember in Evaluating the Models

The exercises above moved from the more common and trivial GIS operations and reached the edge of true spatial modeling. I hope the reader has developed a feel for just how powerful these systems are. Unfortunately, with power comes responsibility, and a GIS modeler must develop judgment.

First, we must remember where the data came from. For instance, the elevation map Rupert models the surface of the land quite directly. It is a rectangular grid of cells, columns and rows, with an elevation value for each cell. This is the format that you would get for a rigourous and direct sampling such as a radar image from an airplane. Rupert is not that good; it has a long and troubled past. In fact, it is

• derived from contour lines,
• digitized as vectors,
• analyzed for inflection points,
• output as a point file,
• compiled as a Triangulated Irregular Network (TIN),
• enquadded in SPANS, and finally
• exported as a raster to use in IDRISI.

29. bc) Display Bayside and window into the Cabrits Peninsula. You derived Bayside from Rupert in an exercise above, so the zero elevation line of Rupert is the coast of Bayside. Overlay the coast line vector. Window in even closer, a little east of the cruise ship terminal, and you see that the fit is not perfect.

I have found places where the coastline vector is 15 meters from the coast shown by the raster. That is nearly two raster cells difference, and it serves to remind us that this data has come through different mathematical models.

The coast vector was digitized from a line on a plastic sheet that served as a master copy for producing the paper maps we use. This does not necessarily mean that the vector line for the coast is perfect. The shore on the south side of Cabrits is rocky. Waves rise and fall. The vector coastline generalizes all of this to a series of points, called "nodes," and the lines between them. Although the nodes are registered with precision, that does not mean that we can know the coastline with precision along the segment between any two adjacent nodes. Maybe the coastline passes inside of the straight line, maybe it passes outside. Indeed, there is so much generalization that the cruise ship terminal does not show at all.

In these digital files which we are using, we rely on the interpretation done by the original cartographers. The paper maps were produced with high professional standards, but interpretation was necessary, and interpretation means simplification. A common standard of precision is that a line on a paper map should not be misplaced by more than one-half of a millimeter. For a map at the scale of 1/25,000 that comes to 12.5 meters on the ground.

The vector of the coastline is in two dimensions. The three-dimensional equivalent is a TIN, which models the land surface as a set of nodes in X,Y, and Z, with little flat triangular facets stretched between them. Now, we know that our elevation data derives from mechanical digitizing of the contour lines, the same process that produced the coastline vector. But the elevation data passed through a TIN model before being enquadded in SPANS and exported as a raster.

30. Display the map of slopes, Rupslope And window in over the central Picard River. The texture of the slope map is distinctly angular. If you look closely, you can see echoes of the triangle of the TIN. They are blurred, because the process of enquading the elevation data and exporting it as a raster caused the computer to interpolate values.

It is common practice to run a smoothing filter over DEMs to remove the angularities of the TIN. Is that really an improvement? The triangles look artificial, but a smoothed surface is artificial too. After several smoothings, the triangles will be hard to distinguish, but the DEM has been altered.

The second element of judgment is flexibility. Despite the name of the exercise in the previous section, there are no "unique solutions" in modeling. The modeler makes approximations, assumptions and guesses from the very start. It is important to remember that.

We do not model to get a single answer. We model to see how a system behaves under differing conditions. Each model is a trial; the computer makes it possible to do more trials so that we have alternative plans ready when they are needed.

The next section, SPANS and IDRISI: Special Workshop, has been extracted to a separate file.

This document was prepared for the CDMP and the workshop was presented by Ross Wagenseil, [email protected]

References

These two workshops were organized around the SPANS and IDRISI software packages, which are distributed under license. For licenses or manuals, contact:

• The IDRISI Project, Clark Labs for Technology and Geographic Analysis, Clark University, 950 Main St., Worchester, MA 01610-1477 USA
• SPANS, TYDAC Research, Inc., 2 Gurdwara Rd., Suite 210, Nepean, Ontatrio, Canada K2E 1A2, (613) 226-5525, (613) 226-3819, [email protected]

Parameters and rules-of-thumb for water works were garnered from:

• Steel, E.W., and McGhee, T.J., Water Supply and Sewerage, McGraw-Hill, 5^th ed. 1979.
• Linsley, R.K., and Franzini, J.B., Water-Resources Engineering, McGraw-Hill, 3^rd ed., 1979.
• Peavy, H.S., Rowe, R.D., and Tchobanoglous., G., Environmental Engineering, McGraw-Hill, 1985.

Colour Illustrations

Figure 1. The town of Portsmouth on the northwest corner of Dominica has room to expand and a ready source of fresh water. Water supply is by gravity flow from an elevation of 360 feet on the Picard River. Place names on this map are referred to in the exercises.

Figure 2. This elevation map is the same as the one used for the background of Figure 1, but it is rendered with a different palette. As a result, the image comes close to overwhelming our ability to interpret it visually. The elevation data is in the Geographic Information System whether we see it or not; contours and colour keys are not needed except for exhibition. This dense data field is a first requirement of many GIS models.

Figure 3. The colours shown range from green (easy construction) to magenta (nearly impossible). The light blue of the bay represents a special value that forbids the program from routing through those cells. The surface analysis of the elevation data is the primary source of this model of friction factors for civil works. Construction difficulty is modeled as a constant plus a factor figured on the tangent of the slope. Other considerations are sea water, an absolute barrier, and the existing road, which is figured to be one-fifth as difficult as undeveloped ground of the same slope. By visual inspection, it appears plausible that the cheapest route for a pipe line would be parallel to the road, passing around the south side of the hill.

Figure 4. Cool green is low cost; yellow, red, and magenta are progressively expensive. Relative costs are figured by working outward from the end of the existing water pipe, doing a summation of the distance, multiplied by the factors from the friction map. The corridor along the highway shows a relative advantage, but it is not as great as expected. Once the cost surface has been calculated, it is a simple matter for the computer to search backwards. If a water tank were to be sited on the hilltop above the possible house sites, the best path for a water line might be down the back side of the hill.

Figure 5. Returning to the elevation map and overlaying the results from the model, we might not be convinced. It appears that the pipe line would cut across the road twice, pass over two ridges and valleys without turning, and climb quite a steep slope. Perhaps the slopes were not weighted heavily enough, or perhaps the road corridor is cheaper than the model allowed. No matter: it is easy to alter the assumptions and run the model again.