Before we start, the weighted analysis function used by Della Bonna can be found in ArcGIS in the Arc Toolbox under Spatial Analysis, Overlay, and Weighted Overlay.
Spatial Analyst-using the Raster Calculator
"When in doubt, use mapcalc!"
In ArcGIS, Spatial Analyst is the raster processing extension. You can perform many different functions, including writing code with ArcObjects,and using various 'canned' functions. But the primary tool is the Raster Calculator. There are many map algebra operators such as Arithmetic, Boolean, Relational, Logical, and Assignment.
Map algebra uses math-like expressions and return numeric values (integer or floating point) for each cell in the output grid. The detailed usage and explanation of each operator and function are available in Help file, and in the Raster Calculator itself.
Please note that the result of any operator or function must be a new output grid and each operator or function acts on all cells in a grids. This is called cell-based operators or functions.
When working with raster data, we need to keep in mind some differences with vector. Raster data needed to be concerned with the cell size, the extent of the analysis, any masks that are used, etc.
Using the Raster Calculator. This is the most powerful tool in the GIS world. It is an implementation of a concept proposed by Dana Tomlin in 1990 called the map calculator (go here for a good overview of the history of GIS). The Map Calulator uses Map Algegra expressions to process raster GIS (or remote sensing) data. This is a new language for you to learn, so accept some frustration in the beginning, but in the old GRASS days we used to say "When in doubt, use mapcalc!" as you can do about anything once you learn how.
Arithmetic operators ( + , - , * , / ) allow addition, subtraction, multiplication, and division. For example, three rasters measuring three different types of fire risk could be added to create an overall risk analysis raster. Arithmetic operators can also be used to convert values from one measurement to another (e.g., feet x 0.3048 = meters).
Relational operators ( < , > , == , <= , >= , <> ) allow you to build logical tests, returning values of true (1) and false (0). For example, this type of operator can be used to find vegetation "equal to" Sierra-type mixed coniferous forest.
Boolean operators (and, or, not) allow you to chain logical tests. Like relational operators, Boolean operators return values of true and false. For example, you could find all slopes that are "greater-than" 45 degrees "and" that have an elevation that is "greater-than" 5000 meters.
Logical operators (DIFF, IN, and OVER) have different results depending on the operator:
A DIFF B: If a cell value in raster A and raster B are different, the cell value in raster A is returned. If the cell values are the same, the value zero is returned.
A IN {value list): If a cell value in raster A is in the value list, the cell value in raster A is returned. Otherwise, NoData is returned.
A OVER B: If a cell value in raster A is not equal to zero, the cell value is raster A is returned. Otherwise, the cell value in raster B is returned.
For example:Convert a DEM from meters to feet. Say you are working for the COE and they want a DEM that comes in meters converted into feet. Start ArcMap and go to the Chapel Hill or Hillsboro(ugh) datasets. Load the ArcGIS Spatial Analyst extension and make the toolbar visible. Then click on Raster Calulator. This is what you will see.
You can expand or contract the calulator using the << or >> box at bottom right.
Feaked out yet? Don't worry.
Building expressions is difficult and confusing at first. Getting help with building expressions is easy. The "About Building Expression" button gives some basic pointers. You can find specific usage help for various expressions by typing the expression into the window, right clicking, and selecting usage...
You will see this little box A LOT! Don't worry, we all do...
Now lets use the raster calculator. In your Raster Calculator you will notice the layers list shows the raster data available. The keypad area contains the numbers and commonly used math operators. In the layers list, double click on your elevation DEM layer to add it to the expression box. Click the multiplication (*) operator. Now enter the number: 0.3048
Your expression should read:
[Elevation] * 0.3048
click Evaluate
A new temporary layer named calculation is added to your data display and to the TOC.
Next, choose a color ramp color table and right click on the new calculation layer. Right click on the new calculation layer to open the layers properties dialog box. Click on Symbology and from the show list, choose stretched.
Each cell value in the new layer contains the product of a raster calculation, where the original DEM value in meters is now converted to feet. As the multiplier contained decimal points, this is a floating point type dataset.
Now convert this floating point raster to integer. The INT Map Algebra function lets us do this. In the Raster Calculator, enter the following expression;
INT([Calculation]) and click Evaluate
A new layer, named Calculation2 is created where cell values are truncated (not rounded) at the decimal place. You can create permanent output grids easily. Open the Raster Calculator and enter the following expression:
Newelevft = INT([Elevation]* 0.3048)
This creates a new permanent raster file called Newelevft. You could now create a vector contour file in Spatial Analyst that would represent contours in feet.
So: We build raster calculator expressions by typing into the expression box, clicking layers in the layers list, and clicking buttons on the keypad. Easy, AND VERY, VERY POWERFUL.
Let's try some more. Look at your elevation data and see the range of values. To create a new map of all locations above or below a certain elevation you could type:
Goodelev = [Elevation] < 200
then click Evaluate and Look at it.
This makes a new raster map where all cells with values less than 200 are 1 (true) and others are false (0). Try this with different elevation ranges. We can do the same for slope. Look at your overall slope map, and then start the Raster Calculator and type:
Goodslope = [slope] < 5 then click Evaluate and look at it. This shows us all areas with slope of less than 5 (degrees or percent, whatever it is). Easy.
NOW find all areas that are both less than 200 feet elevation and less than 5 degrees in slope. We do this using the Boolean operators. Boolean Operators use Boolean logic to conduct analysis on raster data. Results are true (1 in a cell) or false (0 in a cell) depending on the operators used. Open the raster Calculator again and enter the following Map Algebra expression:
Goodfarm = [goodelev] AND [goodslope] and click Evaluate. The new layer has values for each cell of 1 (true) for cells that have BOTH elev less than 100 meters and slope less than 5 degrees. Easy! All other cells are False (o).
We can do this all in one expression:
Goodfarm = ([slope] < 10) AND ([elevation] < 200)
click Evaluate. Like all mathematical expressions, the Map Calculator works from inside expressions out. We can create expressions with dozens of variables (raster layers) to create complex queries. Where are all places with a certain 4 soil types, above this elevation, facing southwest, on public lands, greater than 300 meters from the private lands, within 100 meters of a permanent drainage, in the line of sight of this hillfort. Boom! or as Emeril would say "BAM!" there it is. Powerful stuff. But this is just the start.
Combinatorial Operators
The Combinatorial Operators (CAND, COR and CXOR) are the raster equivalents of vector overlays. CAND only outputs combined attributes for cells where both inputs are true.
Lets add a buffer from permanent water, so go to the ArcToolbox and under Analysis Tools, go to Proximity, and then click on the hammer for Buffer. Choose the maj_riv_CH file and do a 1,000 meter buffer and call it diststreams1.
Open the Raster Calculator and enter the following expression:Farmsoil = [goodfarm] CAND [diststreams1] and Click Evaluate. What this does is not immediately apparent, so right click on the TOC at right and look at the Attribute table. Right click on Farmsoil and click Open Attribute Table. Left Click on the gray area on the far left of each row of data to see it on the screen. These are the different distance to water zones that are found in the goodfarm sites. Try doing this with a mylar and quad map!
Boolean operators
Boolean operators use Boolean logic (TRUE or FALSE) on input
rasters on a cell-by-cell basis. Output values of TRUE are written as 1 and FALSE as 0.
Boolean operators: And, Or, Xor, Not
And (&): Finds where values are true (nonzero) in the cells of both input rasters.
Or ( | ): Finds where nonzero values are present in the cells of one or both input rasters.
Xor (!): Finds where nonzero values are present in the cells of one input raster or another input raster, but not both.
Not (^): Finds where nonzero values are not present in the cells of a single input raster.
Functions
Functions provide most of the power in the map calculator. There are 186 of these (as opposed to 29 operators). This is how slope, aspect, hillshade, and all remote sensing functions are actually done. To create a slope map, you can use the pulldown functions in Spatial Analyst, but what is ACTULLY happening is the following map calculator function:
SLOPE(<grid>, <z_factor>, {degree | PERCENTRISE})
In a previous lab we have created a hillshade display using the transparency tool, but we can create a permanent hillshade layer using the Raster Calculator
The syntax is:
HILLSHADE<grid>, {azimuth}, {altitude}, {ALL | SHADE | SHADOW} {z_factor}
Open the Raster Calculator and enter the following expression:
HILLSHADE2 = Hillshade ([dem_chphl])
You will have to assign symbology yourself when you use functions, so double click on Hillshade and open the Layer Properties dialog. Click on Symbology tab, and from the show list, choose Stretched, with a black to white color ramp. And click OK
Hillshades have many analytic uses as well. The range of 0 to 255 in Hillshade actually represent solar exposure, very useful in agriculture, siting of ske slopes, and archaeology. We can 'Place' the Sun (no wonder GIS people get swolen heads!) in the true position for any day and year in history. Lets say we have data for a given location that places the Sun in these mean locations for 3 months, lets compute the
Month Azimuth Altitude
Nov 186.6 37.2
Dec 183.1 32.8
Jan 179.3 35.2
Open the Raster Calculator and enter the following expression:Nov = HILLSHADE([Elevation], 186.6,37.2) and click Evaluate. Symbolize it with a stretched black to white color ramp. It will look inverted, but its ok. Now do the same for Dec. and Jan. To see the illumination change over time, arrange them by date in the TOC and turn them off. Now turn them on one by one to see the change. Very cool.
Now lets use the MEAN (not your teacher, the math function) command to compute the mean illumination for each cell for the 3 months. Turn off all layers except Elelvation and close their legends.
Open the Raster Calculator and enter the following expression:
Sunexposure = MEAN([Nov], [Dec], [Jan]) and click Evaluate. Set the same color table in symboloby. This shows the 3 month average solar radiation index for the area.
If we want to find the areas of maximum mean Sun exposure, open the Raster Calculator and enter the following expression:
Goodsun = [Sunexposure] > 150 and click Evaluate.
Now right click on the Goodsun symbol for code 1, go to 'properties for selected colors', and chose a yellow color. Right click on the Goodsun symbol for code 0 and choose No Color
Turn off Sunexposure, and turn on Hillshade and move it below Elevation
The yellow areas on the map have higher avg. Solar exposure.
(You can do the opposite to find good ski runs in the mountains).
The capabilities go on and on. If we want to find good locations for rock shelters, we can run a 3x3 neighborhood analysis on our elevation file to find the local maximum and minimum elevation within a 3 by 3 set of cells. Go to the Spatial Analyst pulldown and pick 'neighborhood statistics'. Choose to do a minimum operation on a 3x3 window of the elevation file, and name it minelev. Then do the same operation and do a maximum operation and call it maxelev.
Then we do a simple expression:
Cliffs = ([maxelev] - [minelev]). This will produce a raster where each cell shows the elevation difference between the local maximum and minimul elevation of cells. Lets see where these are near water.
rockshelters = ([cliffs > 17] and [diststreams1 <= 1000])
Can you understand what this expression will do? Then you could find shelters facing SW, etc. etc. Download the centroids to your GPS and go look.
the following exercise comes from http://web.mit.edu/gis/www/arcgis_dem/
You can combine several DEMs together but you may need to fill gaps between contiguous DEMs, which were caused either by the original construction of the DEMs or the translation of the DEM from other formats format to ArcInfo grids. There are many ways to do this, but no way can perfectly reconstruct this missing data.
The following is a procedure to fill in the missing data with an average of the values surround the missing data:
1. Merge two adjacent DEMs. Using the Raster Calculator, enter:
merge ( [dem1], [dem2])
2. Rename this calculated grid to merge1. Create a grid that consists of the cells that are the means of the 3x3 area around each cell in the original merged grid. Using the "data" option instructs Grid to use all 3x3 rectangles even if one or more cells have no data. You must do this since we need to fill in NODATA values with numbers.
focalmean ( [merge1], rectangle, 3, 3, data )
Result: The cells with no data are dropped and the mean calculated for the remaining cells. Rename this grid merge2.
3. Merge the original merged DEM with the mean DEM.
Order in the merge command is important: the original DEM containing the original data is the first grid and the mean grid is the second.
merge ( [merge1], [merge2] )
Result: The mean grid is used to fill in only where there is missing data in the original. Rename this grid to final_merge.
4. Remove the grid created with the "focalmean" function and the original merged DEMs.
Logical Operators
The OVER operator is useful for updating one raster with another. All non-zero values in the first input are 'pasted' over the cells of the second. You can update landcover maps, etc. with this.
Newland = [goodfarm] OVER [landcover] would work- we don't have a landcover, so don't try it.
The DIFF operator finds the difference between two rasters. This is very useful in time series analysis, like finding places where land use or value has changed between 1980 and 1990. So you could do:
Landdiff = [landcover] DIFF [newland] and this would make all cells that have no difference be 0, and all other cells will be the actual difference.
The IN operator is used to 'select' cells based on their value and write them to a new raster. This operator compares the input values to those provided on the list and, if it finds a match, returns the value of that cell. Otherwise it creates a 0 for that cell. Say we wanted to find where all the soils of a certain type are located in our Goodsoil layer, we can just start the map calulator and type:
Goodsoil2 = [soil] IN [201,204]
Where 201 and 204 are the soil codes for certain soil types. This will find all those cells. There are also CON, SETNULL, and SELECT Operators.
Here are some additional simple examples:
To compute the average of two raster map layers a and b:
ave = (a + b)/2
To form a weighted average:
ave = (5*a + 3*b)/8.0
To produce a binary representation of the raster map layer a so that category 0 remains 0 and all other categories become 1:
mask = a/a
This could also be accomplished by:
mask = if(a)
To mask raster map layer b by raster map layer a:
result = if(a,b)
To change all map values below 5 to NULL:
newmap = if(map<5, null(), 5)
Hydrologic Functions
All of the hydrologic functions in ArcGIS are done using the Raster Calculator. Try some of the Hydrologic Functions to make the stream network for the DEM. Use the help window to see what each of these functions do.
A Flow Direction grid is made from the DEM and a Flow Accumulation grid is made from the Flow Direction grid. Remember to double click on a grid to add it to an equation. Start with this query:
flowdirection([name-of-dem])
Change the name of the new grid (called Calculation) to FlowDir to make it easier to find. Right Click on Calculation and then click on Properties. Click on the General tab. Change the layer name to FlowDir. Now make the flow accumulation grid:
flowaccumulation([FlowDir])
Change the name of the new grid to FlowAcc as you did for the Flow Direction above. Now find all of the streams in the DEM as represented by cells that have more than 400 cells flowing into them:
[FlowAcc ] > 400
Change the name of the new grid to RiverGrid
Create a stream network from this layer. Click on the Spatial Analyst menu then click on Convert then click on Raster to Features. You should see this dialog:
Change the defaults to the appropriate values then click OK. This creates the stream network for the DEM.
And on and on....So, the Raster Calculator is the most powerful tool in GIS, but it is somewhat difficult to get to know. Take the time to read some and try things out. If you cannot find a tool in GIS to do what you need, you can do it using the Map Calculator. 'When in doubt, use MapCalc'!
Here are two major references to Dana Tomlin's work:
Tomlin, C.D. (1990) Geographic Information Systems and Cartographic Modeling, Englewood Cliffs, NJ: Prentice Hall.
Tomlin, C.D. (1991) Cartographic modeling, in Maguire,D. Goodchild, M. and Rhind, D. (eds.) Geographical Information Systems: Principles and Applications, Vol 1: ch.23, pp.361-374.Here is an online article that you might find to be interesting:
Feel free to contact me if you have any questions.
~ Main page ~ Schedule ~ Web Resources ~ Bibliography ~
