User Guide for SAGA (version 2.0.5) Volume 2 By Vern Cimmery November, 2010 i User Guide for SAGA (version 2.0.5) Copyright ©2010 Vern Cimmery Permission is granted to copy, distribute and/or modify this document under the terms of the Creative Commons Attribution–ShareAlike license. A copy of the license can be downloaded from the Creative Commons website at www.creativecommons.org. The license applies to the entire text of this book, plus all the illustrations that are by Vern Cimmery. ii Acknowledgment Most of the System for Automated Geo-Scientific Analysis or SAGA was created and developed by the working-group Geosystem Analysis (formerly associated with Göttingen University and currently with Hamburg University), headed by Prof. Dr. Jürgen Böhner. The current versions of SAGA are mainly due to the creativeness and participation of the core set of developers; namely Rüdiger Köthe, Andre Ringeler, Victor Olaya, Dr. Christian Caro, Dr. Volker Wichmann, Prof. Dr. Jürgen Böhner and, in particular Dr. Olaf Conrad, who shouldered the main programming work. However, SAGA would not have reached this level of sophistication without that multitude of methodical innovations cooperatively worked out by the working-group Geosystem Analysis as a whole in context with national and international environment related research projects. Dr. Volker Wichmann has provided a significant amount of assistance in the production of this guide. His feedback and supportive collaboration have made this User Guide for SAGA (version 2.0.5) feasible and possible. I cannot express enough appreciation for his support. Volume 2 of this User Guide is an expansion of Chapter 7 in the 2007 User Guide for SAGA Version 2. As such, it introduces users to the use of SAGA functions and modules to address hypothetical spatial analysis issues. This volume assumes reader familiarity with Volume 1. Please feel free to e-mail me if you have any questions or suggestions for improvement. My e-mail address is:
[email protected] This User Guide for SAGA (Volume 2) is contributed to the SAGA user community to, hopefully, assist the user in successfully applying the SAGA functions, tools, commands, and procedures in addressing applications specific to spatial analysis. Please feel free to make a copy, reference the document, etc., as you desire. I would appreciate if you gave me credit for the effort I have invested. There is no intention to violate any software copyrights indirectly related to SAGA. You are free: • • • To copy, distribute, display, and use this manual. To make derivative works. To make commercial use of this work. Under the following conditions: • Attribution. You must give the original author credit. iii • • • • Share alike. If you alter, transform, or build upon this work, you may distribute the result. Work only under a license identical to this one. For any reuse or distribution, you must make clear to others the license terms of this work. Any of these conditions can be waived if you get permission from the copyright holder. Your fair use and other rights are in no way affected by the above. For any further information, please contact the author at the following e–mail address
[email protected]. iv Table of Content User Guide for SAGA (version 2.0.5) ................................................... i Acknowledgment ............................................................................................................... iii Table of Content ................................................................................................................. v List of Figures .................................................................................................................. viii Chapter 1 – Introduction to Volume 2 ...................................................................... 1 Description: Chapter 2 – Some SAGA Basics................................................................ 1 Description: Chapter 3 – SAGA Buffer Modules........................................................... 2 Description: Chapter 4 – Creating Data Layers from DEMs in SAGA.......................... 3 Description: Chapter 5 – Generating a Watershed Basins Layer ................................... 4 Description: Chapter 6 – Evaluating Proposed Community Park Sites.......................... 5 Description: Chapter 7 – Using SAGA With Digital Satellite Images........................... 6 Description: Chapter 8 – The SAGA Point Cloud (.spc) Data Type .............................. 7 Description: Volume 1.................................................................................................... 8 Chapter 2 - Some SAGA Basics .................................................................................. 9 List of Data Layers, Tables and Maps In Work Session................................................. 9 Rename a Grid or Shapes Data Layer........................................................................... 13 The SAGA Grid Calculator .......................................................................................... 17 Re-Classifying Grid Data Values.................................................................................. 27 Using the Change Grid Values Module........................................................................ 27 Using the Reclassify Grid Values Module .................................................................... 32 Creating a Vector Layer from a Grid Data Layer ......................................................... 42 Viewshed Analysis........................................................................................................ 47 Chapter 3 - SAGA Buffer Modules .......................................................................... 64 Introduction................................................................................................................... 64 The Shapes – Tools/Shapes Buffer Module .................................................................. 65 Example 1 ..................................................................................................................... 66 Example 2 ..................................................................................................................... 71 The Grid – Tools/Grid Buffer Module.......................................................................... 76 Example 1 ..................................................................................................................... 77 Example 2 ..................................................................................................................... 80 The Grid – Tools/Grid Proximity Buffer Module ......................................................... 86 Example 1 ..................................................................................................................... 87 Example 2 ..................................................................................................................... 91 The Grid – Tools/Threshold Buffer Module ............................................................... 100 Example 1 ................................................................................................................... 102 Example 2 ................................................................................................................... 104 Example 3 ................................................................................................................... 108 Chapter 4 – Creating Data Layers from DEMs in SAGA ................................ 111 Introduction................................................................................................................. 111 Creating Slope and Aspect Maps From a DEM.......................................................... 111 Developing Terrain Form Grid Data Layers............................................................... 117 Creating Contour Shapes Data Layers........................................................................ 123 Chapter 5 - Generating a Watershed Basins Layer .......................................... 128 Preparing a DEM Data Layer for Hydrology Analysis .............................................. 128 v Calculating Catchment Area ....................................................................................... 134 Defining Channel Networks ....................................................................................... 141 Defining Watershed Basins......................................................................................... 148 Chapter 6 - Evaluating Proposed Community Park Sites ............................... 151 Selection Criteria ........................................................................................................ 152 Criterion 1: Identifying the boundary for Grapeview ................................................. 152 Criterion 2: identifying Open Space land use in Mason County and Grapeview....... 156 Criterion 3: Open Space parcels must be a minimum size of 5 contiguous acres ...... 158 Criterion 4: Open Space land parcels with road access .............................................. 162 Criterion 5: The park parcel cannot be adjacent to these land uses: industry, commercial, mineral extraction, and utilities.............................................................. 169 Criterion 6: The park parcel may be adjacent to these land uses: agriculture, commercial recreation, public, rural and residential................................................... 169 Criterion 7: The park is most suitable adjacent to these land uses: open space, forest, rural, and water ........................................................................................................... 170 Criterion 8: The highest amount of ground surface with slopes between 0 and 3 degrees ..................................................................................................................................... 170 Criterion 9: The viewing of these land use classes in the foreground and middle-ground will be minimized: industry, commercial, mineral extraction, utilities. This criterion will be used to compare parcel choices....................................................................... 172 Discussion of Results.................................................................................................. 181 Summary ..................................................................................................................... 185 Chapter 7 – Using SAGA With Digital Satellite Images .................................. 186 Imagery Background Information............................................................................... 187 Contrast Enhancements............................................................................................... 190 Filters .......................................................................................................................... 203 Ratios .......................................................................................................................... 215 Composite Images....................................................................................................... 228 Vegetation Indices ...................................................................................................... 232 The Normalized Difference Vegetation Index............................................................ 233 Ratio............................................................................................................................ 235 The Ratio Vegetation Index (RVI) ............................................................................. 237 Transformed Vegetation Index (TVI)......................................................................... 240 Corrected Transformed Vegetation Index (CTVI) ..................................................... 242 Normalized Ratio Vegetation Index (NRVI).............................................................. 244 Image Classification.................................................................................................... 247 Unsupervised Classification........................................................................................ 248 Supervised Classification............................................................................................ 255 Preparing the Input Polygon Shapes Data Layer ........................................................ 263 Chapter 8 – The SAGA Point Cloud (.spc) Data Type ..................................... 273 Point Cloud Type Background Information ............................................................... 273 Importing a .las File .................................................................................................... 274 Cluster Analysis for Point Clouds............................................................................... 278 How to Delete a Point Cloud Attribute....................................................................... 286 Point Cloud Attribute Calculator ................................................................................ 287 Selecting a Subset of a Point Cloud Layer.................................................................. 290 vi The Shapes – Point Clouds/Point Cloud Reclassifier/Subset Extractor Module ....... 292 Thinning a Point Cloud............................................................................................... 298 Creating a Point Cloud Layer from Other SAGA Data Types ................................... 299 Converting a Point Cloud to a Grid Data Layer ......................................................... 304 Saving a Point Cloud Layer to Shapes Data Layer..................................................... 306 Transforming a Point Cloud Layer ............................................................................. 307 The Shapes – Point Cloud Viewer/Point Cloud Viewer Module................................ 309 APPENDIX I – Volume 2 Example Data Layers ............................................... 314 vii List of Figures Figure 2-1. An example of the ‘Data’ tab area list in the workspace. .............................. 11 Figure 2-2. The ‘Maps’ tab area list in the Workspace window....................................... 13 Figure 2-3. The parameters page for a DEM data layer. .................................................. 14 Figure 2-4. Viewing the storage file name in the ‘Description’ area of the ‘Object Properties’ window. .................................................................................................. 16 Figure 2-5. The Grid Calculator parameter settings page................................................ 18 Figure 2-6. The ‘Grids’ input window for the Grid Calculator module. ......................... 19 Figure 2-7. The two population grid data layers in Map Layout view windows.............. 21 Figure 2-8. The ‘Grid Calculator’ settings page used for the ‘ADD’ example. ............... 22 Figure 2-9. The composite grid data layer for total population. ....................................... 23 Figure 2-10. The initial Grid Calculator parameter settings page.................................... 25 Figure 2-11. The ‘Grids’ input window............................................................................ 26 Figure 2-12. The Grid Calculator parameter settings window used with a Boolean formula...................................................................................................................... 26 Figure 2-13. The slope and aspect grid data layers prior to recoding............................... 27 Figure 2-14. The Change Grid Values module parameter settings page.......................... 29 Figure 2-15. The ‘Lookup Table’ for the slope recode..................................................... 30 Figure 2-16. Comparing the original slope data layer with the recoded data layer. ......... 31 Figure 2-17. The ‘Lookup Table’ for the aspect recoding................................................ 32 Figure 2-18. Comparing the original aspect data layer with the recoded data layer. ....... 32 Figure 2-19. The parameter settings page for the Reclassify Grid Values module........... 33 Figure 2-20. Comparing the original ‘MCschoolDist’ data values to the ‘Reclassified Grid’ values............................................................................................................... 35 Figure 2-21. The Mason County census tract grid data layer. .......................................... 36 Figure 2-22. Using the “range” method to combine several data values.......................... 37 Figure 2-23. A blank ‘Lookup Table’ for entering data class values. .............................. 38 Figure 2-24. The ‘Lookup Table’ with new data classes defined based on elevation criteria. ...................................................................................................................... 39 Figure 2-25. Using the “table” method to create input to a wildlife habitat model.......... 40 Figure 2-26. The color lookup table for the elevation input to a wildlife habitat model.. 40 Figure 2-27. The Vectorising Grid Classes parameter settings page................................ 42 Figure 2-28. The parameter settings page for the ‘MCwaterbasins’ shape data layer...... 44 Figure 2-29. The ‘MCwaterbasins’ shapes data layer default display.............................. 45 Figure 2-30. Choosing a shapes data layer polygon boundary color................................ 46 Figure 2-31. The watershed basin shapes data layer as an overlay on the DEM grid data layer........................................................................................................................... 47 Figure 2-32. The Mason County DEM and watershed basins. ......................................... 49 Figure 2-33. The selected watershed basins. .................................................................... 50 Figure 2-34. The ‘New layer from selected shapes’ parameter window. ......................... 50 Figure 2-35. A new shapes data layer based on selected watershed basins...................... 51 Figure 2-36. The ‘Clip Grid with Polygon’ parameters window...................................... 52 Figure 2-37. The DEM grid data layer for the four watershed basins. ............................. 52 Figure 2-38. The observer location (the black plus symbol). ........................................... 53 Figure 2-39. The Visibility (single point) module parameters window. ........................... 53 viii Figure 2-40. The entries for the Visibility (single point) module parameter settings page. ................................................................................................................................... 54 Figure 2-41. The enlarged DEM displaying the observer location and the ‘Visibility’ grid data layers. ................................................................................................................ 55 Figure 2-42. The ‘Visibility’ grid data layer updated. ...................................................... 56 Figure 2-43. Two perspectives for the Bear Gulch Camp flagpole seen area. ................. 56 Figure 2-44. The ‘3D-View’ properties page. .................................................................. 57 Figure 2-45. The ‘Grid Proximity Buffer’ parameters window........................................ 58 Figure 2-46. The distance buffers for the Bear Gulch Campground. ............................... 60 Figure 2-47. The ‘Reclassify Grid Values’ parameters page............................................ 61 Figure 2-48. The lookup table for reclassifying the buffer zones into foreground, midground and background............................................................................................. 62 Figure 2-49. The visual zone grid data layer. ................................................................... 62 Figure 3-1. Settings page for module Shapes Buffer. ....................................................... 65 Figure 3-2. A portion of the ‘MasonFD3roads’ attribute table......................................... 66 Figure 3-3. The settings used to create ROW zones (buffers) adjacent to roads.............. 67 Figure 3-4. The new ‘MasonFD3ROWs’ polygon shapes data layer............................... 67 Figure 3-5. Settings for the Shapes – Grid/Clip Grid with Polygon module.................... 69 Figure 3-6. The settings for the Geostatistics – Grids/Zonal Grid Statistics module execution. .................................................................................................................. 70 Figure 3-7. The ‘Result Table’ with the ‘Acres’ column included................................... 71 Figure 3-8. The settings for the Shapes – Tools/Shapes Buffer module to create the distance zones. .......................................................................................................... 72 Figure 3-9. The output polygon shapes data layer from the Shapes Buffer module......... 73 Figure 3-10. The attribute table for the ‘PointDistZones’ polygon shapes data layer...... 73 Figure 3-11. The settings used with the Clip Points with Polygons module. ................... 74 Figure 3-12. The new point shapes data layers for each distance zone. Zone 25 is upper left and zone 100 is lower right................................................................................. 75 Figure 3-13. A portion of the attribute table for zone 100................................................ 75 Figure 3-14. Settings page for module Grid Buffer. ......................................................... 76 Figure 3-15. The Grid - Tools/Grid Buffer settings used to create buffers around water bodies. ....................................................................................................................... 77 Figure 3-16. The ‘MClakes’ layer (on the left) and the Grid Buffer output on the right. 78 Figure 3-17. The Grid Calculator parameter settings page.............................................. 78 Figure 3-18. The Grid Calculator parameter settings page for creating a new grid data layer........................................................................................................................... 79 Figure 3-19. The grid layer identifying road areas of interest. ......................................... 80 Figure 3-20. Zoomed in area of the ‘MFDrdsROW’ grid data layer................................ 81 Figure 3-21. Settings page for creating the road ROW grid data layer. ........................... 81 Figure 3-22. The Grapeview Road Maintenance District ROW map............................... 82 Figure 3-23. The Grid Calculator settings for creating a new grid data layer containing the data for the 120’ ROW width road...................................................................... 83 Figure 3-24. The Grid – Tools/Grid Buffer settings for the 60’ wide ROW buffer layer. 84 Figure 3-25. Example of the buffer “overlap” problem.................................................... 85 Figure 3-26. A zoomed in portion of the Grapeview Road Maintenance District ROW map............................................................................................................................ 86 ix Figure 3-27. The stream network for Mason County. ...................................................... 87 Figure 3-28. The Mason County general land use grid data layer.................................... 88 Figure 3-29. The Grid – Tools/Grid Proximity Buffer module settings page used to create the ‘5 acre cell groupings.......................................................................................................... 159 Figure 6-10. Grapeview Open Space areas greater than 5 acres in size. ........................ 160 Figure 6-11. Open Space areas in Grapeview greater than 5 acres in size. .................... 161 Figure 6-12. The attribute table for the ‘GVOpSpPoly’ shapes data layer..................... 161 Figure 6-13. The Table-Calculus/Table calculator for shapes module settings for calculating “ACRES”.............................................................................................. 162 Figure 6-14. Grid Calculator settings for overlaying gravel/paved roads with open space land use in Grapeview............................................................................................. 163 Figure 6-15. Display of ‘GVOpSpRoads’ grid data layer. Site 4 does not have a road contact. .................................................................................................................... 163 Figure 6-16. Settings used for the Shapes-Tools/Shapes Buffer module. ....................... 164 Figure 6-17. Settings for the Shapes-Polygons/Polygon Intersection module execution. ................................................................................................................................. 165 Figure 6-18. Zoomed in portions of the ‘GVOpSpEDGE’ output data layer................. 166 Figure 6-19. The settings used for the Grid-Gridding/Shapes to Grid vector to raster conversion. .............................................................................................................. 166 Figure 6-20. The Grid Calculator settings for creating a grid data layer for land use classes along the site borders. ................................................................................. 167 Figure 6-21. The ‘Zonal Grid Statistics’ module settings............................................... 168 Figure 6-22. The table for land use adjacency for each of the sites................................ 169 Figure 6-23. Grid Calculator settings for creating a grid layer containing less than 3 degree slopes........................................................................................................... 170 Figure 6-24. The Zonal Grid Statistics module settings. ................................................ 171 Figure 6-25. Acres and percent of total site acreage with slopes less than 3 degrees. ... 172 Figure 6-26. The settings for the Shapes-Polygons/Polygon Centroid module.............. 173 Figure 6-27. A zoomed in view of the centroid (red dot) in site #9. .............................. 173 Figure 6-28. The observer point shapes data layer displayed with the DEM grid data layer......................................................................................................................... 174 xii Figure 6-29. The settings for the Terrain Analysis-Lighting, Visibility/Visibility (single point) [interactive] module. .................................................................................... 175 Figure 6-30. An example of the settings for the Grid-Tools/Grid Buffer module.......... 177 Figure 6-31. The Grid Calculator settings for creating the foreground, midground, background grid data layer for observer point 1..................................................... 178 Figure 6-32. The ‘Vzones01’ grid data layer displaying the three distance zones for observer point 1....................................................................................................... 178 Figure 6-33. The viewshed grid data layer for observation point number 1................... 179 Figure 6-34. The Grid Calculator settings used to produce a seen area grid data layer showing foreground, middle-ground, and background areas.................................. 179 Figure 6-35. The output grid data layer showing foreground, middle-ground, and background areas for the viewshed for observation point number 1. ..................... 180 Figure 6-36. The visual impact tables............................................................................. 180 Figure 6-37. Weights assigned to the sites based on applying the criteria. .................... 183 Figure 6-38. The Grid – Gridding/Shapes to Grid module settings page for creating the grid data layer for criteria #4. ................................................................................. 183 Figure 6-39. Using the Grid Calculator to display a composite score for all criteria for each site................................................................................................................... 184 Figure 6-40. The site grid data layer displaying the average weight for each site. ........ 185 Figure 7-1. The full Olympic Peninsula Landsat scene (on left) and Mason County subscene on right. ......................................................................................................... 188 Figure 7-2. Comparing information for the full scene and sub-scene. ........................... 188 Figure 7-3. The full central Arizona Landsat scene (on left) and sub-scene on right..... 189 Figure 7-4. Descriptive information for the full scene and sub-scene............................ 189 Figure 7-5. Basic descriptive information for a grid data layer containing spectral data for band 3 of a Landsat TM sub-scene. ........................................................................ 191 Figure 7-6. The histogram for spectral reflectance values for band 3 (grid data layer ‘MCL720020922_3’).............................................................................................. 192 Figure 7-7. A portion of the tabular version of the histogram in Figure 7-6. ................. 193 Figure 7-8. Functions available in the ‘Classification’ sub-menu. ................................. 194 Figure 7-9. Comparing corresponding zoomed in areas on the input and output following applying the ‘Create Normalised Classification’ option to a grid data layer.......... 195 Figure 7-10. The ‘Create Lookup Table’ window. ......................................................... 195 Figure 7-11. Comparing zoomed in portions of input and output grid data layer after using the ‘Create Lookup Table’ option. ................................................................ 196 Figure 7-12. The ‘[CAP] Colors’ window...................................................................... 196 Figure 7-13. A portion of the lookup table created using an entry for the ‘Count’ parameter of 100. .................................................................................................... 197 Figure 7-14. Comparing zoomed in portions of the before and after grid data layers using a grayscale palette option with 20 classes............................................................... 198 Figure 7-15. The lookup table created using an entry for the ‘Count’ parameter of 20. 198 Figure 7-16. Adding color to the grayscale lookup table................................................ 199 Figure 7-17. Comparing the original grid data layer with the one using the modified lookup table............................................................................................................. 199 Figure 7-18. The “bright-dark” shade choice on the left and “dark-bright” on the right using band 4 of a Landsat TM sub-scene................................................................ 201 xiii Figure 7-19. Example settings for “RGB Overlay”. ....................................................... 202 Figure 7-20. Example composite display using the “RGB Overlay” choice.................. 202 Figure 7-21. The settings page for the User Defined Filter (3 x 3) module. .................. 203 Figure 7-22. The ‘Filter Matrix’ parameter with values to calculate mean for the User Defined Filter (3 x 3) module. ................................................................................ 204 Figure 7-23. Comparing zoomed in areas of the input and output grid data layers........ 204 Figure 7-24. Visual comparison of the input and output grid data layers....................... 205 Figure 7-25. The Grid – Filter/Simple Filter module settings for a 3 x 3 “mean” kernel. ................................................................................................................................. 205 Figure 7-26. Settings used for the example of the Grid – Filter/Gaussian Filter module. ................................................................................................................................. 206 Figure 7-27. Visual comparison of the input and output grid data layers....................... 207 Figure 7-28. The “Standard kernel 1”, “Standard kernel 2” and “Standard kernel 3” options in the Laplacian Filter module. ................................................................. 207 Figure 7-29. Settings for the Grid – Filters/Laplacian Filter module............................ 208 Figure 7-30. Comparing the input grid data layer and two output grid data layers........ 209 Figure 7-31. The “Sobel Detection” directional filters................................................... 210 Figure 7-32. The settings used for the Gaussian Filter module. .................................... 210 Figure 7-33. The values used for the ‘Filter Matrix’ for the first directional kernel. ..... 211 Figure 7-34. The Grid – Calculus/Grid Calculator module settings.............................. 211 Figure 7-35. Viewing the output from applying the Sobel Edge Technique. ................. 212 Figure 7-36. Majority Filter module settings page......................................................... 213 Figure 7-37. Comparing the input and output grid data layers related to the Majority Filter module execution.......................................................................................... 214 Figure 7-38. Comparing data values for the input and output grid data layers related to the Majority Filter module execution........................................................................... 214 Figure 7-39. Two kernels for different illumination directions. ..................................... 215 Figure 7-40. Comparing the outputs using the two kernels. ........................................... 215 Figure 7-41. Spectral Reflectance Curves and TM+ Band Spectral Ranges. ................. 216 Figure 7-42. The settings page for the Grid Calculator module. ................................... 216 Figure 7-43. The Grid Calculator settings for the ratio TM1/TM2 image..................... 217 Figure 7-44. Comparing two TM1/TM2 ratioed images. ............................................... 218 Figure 7-45. Spectral values for bands 1 and 2 and the ratio image TM1/TM2............. 219 Figure 7-46. Comparing two TM3/TM4 ratio images. ................................................... 220 Figure 7-47. Spectral values for bands 3 and 4 and the ratio image TM3/TM4............. 221 Figure 7-48. Comparing two TM5/TM2 ratio images. ................................................... 222 Figure 7-49. Spectral values for bands 5 and 2 and the ratio image TM5/TM2............. 222 Figure 7-50. Comparing two TM3/TM7 ratio images. ................................................... 224 Figure 7-51. Spectral values for bands 3 and 7 and the ratio image TM3/TM7............. 224 Figure 7-52. The Shapes – Tools/Create New Shapes Layer module settings. .............. 225 Figure 7-53. The Shapes – Grid/Get Grid Data for Shapes module settings page. ....... 226 Figure 7-54. Summary of spectral reflectance and ratio image values in the Mason County sub-scene. ................................................................................................... 227 Figure 7-55. Summary of spectral reflectance and ratio image values in the Arizona subscene........................................................................................................................ 227 Figure 7-56. Settings page for module RGB Composite................................................. 228 xiv Figure 7-57. Zoomed in areas on the true color composites for Mason County and Arizona.................................................................................................................... 231 Figure 7-58. Zoomed in areas on the false color infrared composites for Mason County and Arizona............................................................................................................. 231 Figure 7-59. The Grid Calculator setup parameters for creating a grid data layer for the Normalized Difference Vegetation Index (NDVI). ................................................ 233 Figure 7-60. A zoomed in portion of the sample grid data layer outputs for the NDVI. 234 Figure 7-61. Sample point NDVI values. ....................................................................... 235 Figure 7-62. The Grid Calculator setup parameters for creating a grid data layer for the Ratio vegetation index (RATIO). ........................................................................... 236 Figure 7-63. The sample grid data layer output for the RATIO. .................................... 236 Figure 7-64. Sample point RATIO values. ..................................................................... 237 Figure 7-65. The Grid Calculator setup parameters for creating a grid data layer for the Ratio Vegetation Index. .......................................................................................... 238 Figure 7-66. The sample grid data layer output for the RVI. ......................................... 239 Figure 7-67. Sample point RVI values. .......................................................................... 239 Figure 7-68. The Grid Calculator setup parameters for creating a grid data layer for the Transformed Vegetation Index (TVI)..................................................................... 240 Figure 7-69. The sample grid data layer output for the TVI........................................... 241 Figure 7-70. Sample point TVI values............................................................................ 242 Figure 7-71. The Grid Calculator setup parameters for creating a grid data layer for the Corrected Transformed Vegetation Index (CTVI). ................................................ 242 Figure 7-72. The sample grid data layer output for the CTVI. ....................................... 243 Figure 7-73. Sample point CTVI values. ........................................................................ 243 Figure 7-74. The Grid Calculator setup parameters for creating a grid data layer for the Normalized Ratio Vegetation Index (NRVI).......................................................... 244 Figure 7-75. The sample grid data layer output for the NRVI........................................ 245 Figure 7-76. Sample point NRVI values......................................................................... 245 Figure 7-77. Band 3 and 4 spectral and vegetation index values for sample points in the Mason County sub-scene. ....................................................................................... 246 Figure 7-78. Band 3 and 4 spectral and vegetation index values for sample points in the Arizona sub-scene................................................................................................... 247 Figure 7-79. Settings page for the module Cluster Analysis for Grids........................... 249 Figure 7-80. The ‘Grids’ window displaying a list of grid data layers available for the active project selected in the ‘Grid system’ parameter. .......................................... 249 Figure 7-81. Example of statistics output from the Cluster Analysis for Grids module. 250 Figure 7-82. Settings for the Cluster Analysis for Grids module. .................................. 252 Figure 7-83. The “pass” and “change” display............................................................... 253 Figure 7-84. Output from the Cluster Analysis for Grids module.................................. 254 Figure 7-85. The settings page for the module Supervised Classification. .................... 256 Figure 7-86. The ‘Grids’ window displaying a list of grid data layers available for the selected grid system. ............................................................................................... 257 Figure 7-87. Example of ‘Grids’ parameter is for choosing one or more input grid data layers. The default is “No objects”. If the module is to execute successfully, one or more grid data layers must be chosen for input. The ‘Grids’ input. One or more input grid data layers can be chosen as input. The first time you click with the mouse pointer in the value field to the right of the ‘>>Grids’ label, an ellipsis appears in the field. When you click on the ellipsis symbol a list of loaded grid data layers displays in the ‘Grids’ dialog window. The grid data layers in the list are the ones that share the spatial characteristics of the grid system chosen for the ‘Grid system’ parameter. You can choose one or more grid data layers from the list for use with the calculator formula. Figure 2-6 shows an example of how the pop-up list can appear. 18 Figure 2-6. The ‘Grids’ input window for the Grid Calculator module. The list on the left side is all of the grid data layers that are currently loaded for the grid system chosen in the ‘Grid system’ value field. Data layers that are going to be used in a calculation have to be moved to the list on the right side. There are two steps involved. You move a layer from the left column to the right column by first highlighting the layer name, clicking on it with the mouse pointer. You can select more than one layer by holding down the keyboard CTRL or SHIFT keys with additional mouse clicks. Next you click on the button. This button will move the chosen layer or layers into the list in the right column. You can also move data layers back and forth between the two lists by double-clicking with the mouse pointer positioned on the file you want to move. You must pay attention to the order that grid data layers appear in the list on the right. When you refer to a grid data layer in the ‘Formula’ parameter value field, a single text character is used to reference the grid data layers based on their position in the list. For example, if you have four grid map layers in your ‘Input’ list, they will be referred to by “a”, “b”, “c”, and “d”. The letters are assigned to the grid data layers according to their order in the list in the right column, from top to bottom. You can use letters up to “z”. Notice that the bottom two buttons in the middle of the ‘Grids’ window are labeled ‘Up’ and ‘Down’. These buttons allow you to move a highlighted grid data layer in the list, on the right, up or down in the list. When you move a layer up or down in the list, the alphabetic reference letter for it correspondingly changes. The ‘Grids’ label displays “2 objects (MCnonPovPop.dgm, MCpoverty.dgm)”. The objects are the two input grid data layers. The formula entry is “a+b”; adding layer “a” (the population above the poverty level) with “b” (the population below the poverty level). The ‘Grids’ parameter list. The ‘MC1AllTer’ layer is the single entry for the parameter, thus, it will be referred to by the letter “a”. The logical condition being evaluated is equality. In the formula the logical condition is “eq”. The equal condition being checked is for the data value of the grid cell in layer “a” compared to the number “-7”. If the data value for the grid cell on the input layer is equal to the value “-7”, then the condition is true. The condition to be checked is represented by “(a,-7)” in the formula above. The rest of the “ifelse” statement reads as follows: if the condition is true, then put the cell value of grid “a” in the corresponding cell of the output layer. If the condition is false, i.e., the value for the cell does not equal the value “7”, then enter a “0” in the corresponding cell in the output layer. These actions are represented by “,a,0”. Using this formula allows me to easily create a new grid data layer representing channels from a grid data layer containing data values for several different topographic landforms. I start the process by executing the Grid - Calculus/Grid Calculator module. The parameters page for the Grid Calculator is displayed (Figure 2-10). Figure 2-10. The initial Grid Calculator parameter settings page. The first parameter in the list is ‘Grid system’. The grid system that I choose for this parameter must be the one that the grid data layer ‘MC1AllTer’ is a part. There can be more than one grid system loaded for a work session. I choose the ‘Grid system’ by clicking in the value field to the right of the ‘Grid system’ parameter. A small triangle appears in the value field and a pop-up list of the grid systems loaded for this work session displays. I click with the mouse on the appropriate grid system definition. 25 The ‘>>Grids’ parameter is where I choose the input grid data layer. As discussed earlier, more than one grid data layer can be chosen for input. In this example, there is only one data layer needed. When I click in the value field to the right of the ‘>> Grids’ label, an ellipsis symbol appears in the field. When I click on the ellipsis symbol the ‘Grids’ window will appear (Figure 2-11). Figure 2-11. The ‘Grids’ input window. The list on the left side is all of the grid data layers that are currently loaded in the SAGA work session that share the spatial characteristics of grid system identified for the ‘Grid system’ parameter. You move a grid data layer from the list to the right column by clicking on the layer name with the mouse pointer so the name is highlighted. If more than one grid data layer is needed as input, I could make these selections by pressing the keyboard CTRL or SHIFT keys and additional mouse clicks before clicking on the move button. This moves the selected layer into the list in button. Next I click on the the right column. I will only be using one grid data layer (‘MC1AllTer’). It appears first in the list. I choose it and move it to the right column. There is a shortcut for moving layers from one side of the list to the other. If you double-click on the layer name it will move to the other side of the list. I press the ‘Okay’ button and return to the ‘Grid Calculator’ settings window. I will leave the ‘Categorical Grids’ parameter. The output table, ‘Shapes’ parameter value field. In order for the input shapes data layer to be chosen, the layer must be loaded in the current work session. If it is not loaded, you can use the ‘Shapes/Load Shapes’ command in the Menu Bar ‘File’ drop down menu to load it. When you click in the value field to the right of the ‘>>Shapes’ parameter, a list of available shapes data layers displays. Clicking on the shapes data layer name will choose it for the value field. Shapes data layers have a linked table called an attribute table that provides characteristics for the layer objects. One or more of the attribute fields may contain data that can be used as buffer distances or zone widths related to the layer objects. The ‘Buffer Distance (Attribute)’ parameter is used to choose such an attribute to provide this data for analysis. This is an optional parameter. In the absence of a buffer distance attribute or choosing not to use one, other module options will provide the required data. The ‘Shapes’ parameter. An attribute stored in the linked attribute table will be used to provide the module the data values for buffer width. The buffer width varies by surface type. As discussed earlier, this attribute is named ‘BUFFER’ in the attribute table. The module is directed to use the information in the ‘Buffer Distance Attribute’ parameter by the choice of the “attribute field” option for the ‘Buffer Distance’ parameter. The other parameters in the ‘Options’ section will be ignored. The output is named, by default, ‘MasonFD3roads [Buffer]’. I renamed the output ‘MasonFD3ROWs’. Figure 3-4. The new ‘MasonFD3ROWs’ polygon shapes data layer. 67 The ‘MasonFD3ROWs’ polygon shapes data layer contains a single polygon object. The object consists of multiple parts. The full extent of the new layer is displayed on the left in Figure 3-4. A zoomed in area appears on the right. The ROW edges are displayed in black and the ROW is filled with pale yellow. I have overlain the roads in the zoomed in map view. Paved roads are in red, gravel are blue, and brown is used for dirt. The next step is to use the ROW polygons to identify and quantify the amount of forestry land use within the ROW widths. There is not a general land use shapes data layer available for the road maintenance districts. But there is a general land use grid data layer. There is a module that uses a polygon shapes data layer to “clip” a grid data layer creating a grid data layer with data for the geographic area covered by the polygon object(s) on the input shapes data layer. This should work for my objective but I need to do some preprocessing with the general land use grid data layer that is available. The grid data layer available for general land use has a cell size of 90’. The road ROW’s are all less than 90’. Ideally, I would like a spatial resolution that is less than the smallest ROW width. SAGA has a module that is used to resample the larger cell size to a smaller resolution more compatible with my requirements. I can use the Grid – Tools/Resampling module to convert the ‘MasonFD3genLU’ grid data layer from a 90’ cell size resolution to a 30’ cell size resolution. After executing the Resampling module, I rename the output grid data layer ‘MFD3genLU30’ to distinguish it from the larger cell size version. The Shapes – Grid/Clip Grid with Polygon module will be used to create a grid data layer containing only the land use classes within the ROW (the buffer generated from the Shapes Buffer module). I can use the Zonal Grid Statistics module to generate a table showing area of land use within the ROW. Figure 3-5 displays the settings I use for execution of the Shapes – Grid/Clip Grid with Polygon module. 68 Figure 3-5. Settings for the Shapes – Grid/Clip Grid with Polygon module. The grid data layer for land use within the maintenance district is ‘MFD3genLU30’ and is chosen for the ‘>>Input’ parameter. The ‘Grid system’ chosen is the one the land use layer is a part. The polygon shapes data layer representing the ROWs is chosen for the ‘>>Polygons’ parameter. Remember this shapes data layer was generated using the ‘BUFFER’ attribute to create a single polygon for the road ROW. The output is a new grid data layer, in a new grid system, capturing land use class data within the road ROW for all roads in the Grapeview road maintenance district. I rename the output grid data layer ‘ROWlu30’. The Geostatistics – Grids/Zonal Grid Statistics module can be used to tabulate the number of cells for each class or category on a categorical grid data layer. That is one of modules basic functions in addition to several much more sophisticated tabulations. Figure 3-6 displays the settings used for the Zonal Grid Statistics module to produce its default output ‘Result Table’. 69 Figure 3-6. The settings for the Geostatistics – Grids/Zonal Grid Statistics module execution. As noted above, I renamed the output from the Shapes – Grid/Clip Grid with Polygon module ‘ROWlu30’. This grid data layer was chosen for the ‘>>Zone Grid’ parameter in the settings. There are other optional inputs that could be used for more complex statistic summaries. This one is a simple one. The output from this module is a table called ‘Zonal Statistics’. The ‘Zonal Statistics’ table lists number of cells in each land use class (0 through 12) within the ROW polygon. The Table – Calculator/Table Calculus module can be used to divide the cell count by the number of cells making up an acre to determine the number of acres per land use class. It takes 48.4 cells to make up an acre. A new column is generated in the ‘Zonal Statistics’ table that contains the result of dividing the number of cells by 48.4 to create an ‘Acres’ column. Figure 3-7 displays the ‘Zonal Statistics’ table with the ‘Acres’ column created. 70 Figure 3-7. The ‘Result Table’ with the ‘Acres’ column included. LU Class 9 in the table is the land use class for forestry. The table shows that about 45 acres of forestry land use are within the roads ROW in the maintenance district. As noted earlier, the “Acres” column in the ‘Zonal Statistics’ table was generated using the Table – Calculus/Table Calculator module. In the Table Calculator module, the table fields or columns are referred to by an alpha representation of their position going from left to right. The first column is referred to by “a”, the second by “b”, etc. Thus, based on the original table, the equation “b/48.4” was used to calculate acres. The values in the “Count” column represent the number of grid cells for that category. It takes 48.4 grid cells to make an acre. The equation says to divide the value in column “b” (the ‘Count’ column) by the value 48.4. The output from the equation was placed in a new table column named “Acres”. Notice in the title bar containing the name of the table, in Figure 3-7, that the new table output by the Table Calculator is ‘Result Table [b/48.4]’. The equation has been concatenated to the default table name. Example 2 The Reach Harbor Marina has applied for a permit for marina improvements including the addition of new pilings. Pile driver equipment is quite noisy and can be disruptive. The County has requested the Marina to develop a contact list of landowners that could be potentially impacted by the noise. The County wants an address list based on distance zones around the improvement site. They want a list for people within ½ mile, from ½ mile to 1 mile, from 1 mile to 1-1/2 miles, and from 1-1/2 to 2 miles. Each of the distance zones is ½ mile wide. The County has provided a point shapes data layer that contains all the addresses in the county as point objects. The first task is to identify the point that identifies the location of the marina project. 71 The ‘MasonAddresses’ point shapes data layer was displayed in a map view window. Using the ‘Action’ tool from the Tool Bar, the point representing the address of the marina was located, selected, and highlighted. The Shapes – Tools/New layer from selected shapes module was then executed, using the ‘MasonAddresses’ point shapes data layer (with the point selected) as input. The output from the Shapes – Tools/New layer from selected shapes module is a new point shapes data layer that I renamed ‘PointSource’. The Shapes – Tools/Shapes Buffer module is executed to create a polygon shapes data layer with the 4 one-half mile concentric zones centered on the marina project. The settings page for this execution is displayed in Figure 3-8. Figure 3-8. The settings for the Shapes – Tools/Shapes Buffer module to create the distance zones. The input for the ‘>>Shapes’ parameter is the point shapes data layer containing the single point representing the location of the marina project. The intent is to define four one-half mile wide zones going out to two miles from the marina. The “fix value” option is chosen for the ‘Buffer Distance’ parameter since the ‘Buffer Distance (Fixed)’ and the ‘Number of Buffer Zones’ parameters will provide the inputs. Two miles is 10560’. This is the value entered for the ‘Buffer Distance (Fixed)’ parameter. The number of zones is 4 (i.e., 4 one-half mile wide zones within the two miles). The number 4 is entered for the ‘Number of Buffer Zones’ parameter. 72 Figure 3-9 displays the output polygon shapes data layer. The output was renamed to ‘PointDistZones’. Figure 3-9. The output polygon shapes data layer from the Shapes Buffer module. The map view window in Figure 3-9 displays the one-half mile distance zone boundaries in black. The orange circle outlines are address points (the circle size is exaggerated at this display scale). The center is the point representing the marina project. I used the Shapes – Points/Count Points in Polygons module to add the expected number of addresses for each of the four zones into the attribute table linked to the ‘PointDistZones’ shapes data layer. The attribute table is displayed in Figure 3-10. Figure 3-10. The attribute table for the ‘PointDistZones’ polygon shapes data layer. The first ½ mile distance zone is ID 2, Zone 25. The Points field indicates it has 216 addresses in it. Zone 50 is the second ½ mile zone, zone 75 the third, and zone 100 the last. 73 The Shapes – Points/Clip Points with Polygons module is used to create a table for each of the distance zones with the potential impact addresses. Figure 3-11 displays the settings used for this execution. Figure 3-11. The settings used with the Clip Points with Polygons module. The input point shapes data layer is the county layer with the address points (‘MasonAddresses’). The polygon shapes data layer input, ‘PointDistZones’, is the layer containing the four polygons, one for each of the one-half mile wide distance zones. The “ZONE” attribute is chosen for addition to the attribute table for the output point shapes data layer. The “ZONE” attribute is the one identifying each zone: 25, 50, 75, and 100. Figure 3-12 displays the output. The output is a new point shapes data layer for each distance zone. The ‘PointDistZones’ polygon shapes data layer is overlain on each of the four output layers. 74 Figure 3-12. The new point shapes data layers for each distance zone. Zone 25 is upper left and zone 100 is lower right. Figure 3-13 shows a portion of the attribute table linked to the layer for zone 100. Figure 3-13. A portion of the attribute table for zone 100. 75 The attribute table in Figure 3-13 contains 209 records. The records include the address for each point as well as other data related to the property. The main objective is to identify the addresses of properties that can be potentially impacted by noise from the project. This has been accomplished. The Grid – Tools/Grid Buffer Module The Grid – Tools/Grid Buffer module is the first of three buffer creation modules that work with grid data layers. The settings page for the module is displayed in Figure 3-14. Figure 3-14. Settings page for module Grid Buffer. The input ‘>>Features Grid’ parameter is where you choose the grid data layer that contains the features the buffer will be created around. The entry for the ‘Grid system’ parameter must be the grid system the input grid data layer is a part. Depending on the options selected, the input grid data layer could also be used for identifying a variable width buffer using the data values portraying the features. The output for the module is a grid data layer displaying the buffers created for the features on the input grid data layer. The default entry for the value field to the right of the ‘Features Grid’. The need is to delineate a buffer with a fixed width of 105 meters. The width is entered in the ‘Distance’ parameter value field. The “Fixed” option is chosen for the ‘Buffer Distance’ parameter. The output grid data layer is renamed ‘MClakesBuffers’ and displayed on the right in Figure 3-16. 77 Figure 3-16. The ‘MClakes’ layer (on the left) and the Grid Buffer output on the right. The two map view windows in Figure 3-16 are showing zoomed in portions of the input ‘MClakes’ grid data layer (on the left) and the output grid data layer produced by the Grid - Tools/Grid Buffer module. The buffer around the water bodies is displayed in red. I will use the Grid Calculator to create a composite grid for the buffer grid layer for lakes and the roads data layer. Figure 3-17 displays the Grid Calculator parameter settings page settings for creating the composite of the two grid data layers before I click on the ‘Okay’ button. Figure 3-17. The Grid Calculator parameter settings page. The composite grid data layer (named ‘MCbuffersRoads’) created by the Grid Calculator contains six data values. The ‘MCRoads10’ grid data layer has two data classes: 0 for no roads and 10 for road. The ‘MClakesBuffers’ data layer contained three values: 0 for no lakes or buffers, 1 for buffer, and 2 for lakes. The data classes in the composite grid data layer are: 0 = No road, lake or buffer 1 = Lakes buffer 78 2 = Lakes 10 = Roads 11 = Road in Buffer 12 = Road in lakes The data value ‘11’ identifies the cells where roads and lake buffers coincide. These are the cells the research group is interested. You might notice that the grid data value 12 represents roads and lakes but no buffer. Some of these cells are where bridges and causeways exist. Some might also indicate a data error in the roads mapping. I am going to use the Grid Calculator module to create a new grid layer containing only the cells having a data value of 11. Figure 3-18 displays the Grid Calculator parameter settings page that will create this new grid data layer. Figure 3-18. The Grid Calculator parameter settings page for creating a new grid data layer. The formula checks the combined lakes buffer and roads grid layer for the data value 11. All cells having the value 11 are retained in a new grid layer recoded with the value “1”. Figure 3-19, on the left shows this new grid data layer with the shapes data file for roads as an overlay. The roads are in white and the cells where the lakes buffer and a road coincide are in red. The right portion of this graphic is an enlargement of the area around Lake Cushman. 79 Figure 3-19. The grid layer identifying road areas of interest. The buffer/road grid cells identified for the research group are in red. They do not show up very well in the map view window on the left, but you can see in the enlarged map view window on the right, the red cells showing up adjacent to Lake Cushman. Example 2 The Grapeview Road Maintenance District is preparing the budget for next year. Management activities vary by surface type. The three surface types in the district are pavement, gravel, and dirt. Activities take place within ROW management zones. This zone is 120’ either side of paved roads, 90’ for gravel, and 60’ for dirt. Activities include planting vegetation, brush cutting, removal of “danger” trees, etc. They have requested a map of their district that displays the ROW widths for all district roads. A grid data layer exists, named ‘MFDrdsROW’, that contains all roads for the district. The data values identify the width of the zone on either side of the road, by road type. The data values are 60, 90, and 120. Roads with a ROW of 60 are coded with the data value 60, etc. This data layer is displayed in Figure 3-20. This data layer identifies the ROW widths but the actual widths are not depicted. 80 Figure 3-20. Zoomed in area of the ‘MFDrdsROW’ grid data layer. The Grid - Tools/Grid Buffer module will be used to create the road ROW map requested by the maintenance district. Figure 3-21 displays the settings page for creating the road ROW grid data layer. Figure 3-21. Settings page for creating the road ROW grid data layer. 81 The buffer created by this module does have variable width but there is no differentiation within the buffer cells for the ROW class. All buffer cells are coded with 1’s. I can use the Grid – Calculus/Grid Calculator to create an individual grid data layer for each ROW class. Then I can use the Grid – Tools/Grid Buffer module to create a buffer for each ROW class. The Grid Calculator can be used to recode the 1’s to 60’s, 90’s, and 120’s and then add them together into the same grid data layer for a final map. The final map is displayed in Figure 3-22. Figure 3-22. The Grapeview Road Maintenance District ROW map. Here is a summary of the process I used. Three new grid data layers were created, one for each road ROW class. Figure 3-23 displays the Grid Calculator settings used for creating the grid data layer containing the data for the 120’ ROW width road. 82 Figure 3-23. The Grid Calculator settings for creating a new grid data layer containing the data for the 120’ ROW width road. The equation checks for the value 120 on the input grid data layer (‘MFDrdsROW.sgrd’). If the value is found, it is output to the corresponding cell on the output grid data layer. If it is not found, a zero is stored in the cell on the output layer. This module was executed three times to create the three new layers: Rd60, Rd90, and Rd120. Next the Grid Buffer module was used. The Grid – Tools/Grid Buffer module was executed once with each new layer. The settings used for one of the three output layers are displayed in Figure 3-24. 83 Figure 3-24. The Grid – Tools/Grid Buffer settings for the 60’ wide ROW buffer layer. The ‘RD60’ layer is the input for the ‘>>Features Grid’ parameter. The data values for the road width in the input layer will be used to define the ROW width on the output grid data layer. The “Cell value” option for the ‘Buffer Distance’ parameter instructs the module to use those values. This module was executed three times to create three new layers: ROW60, ROW90, and ROW120. Thus, the ‘ROW60’ layer contains the buffer for dirt roads with a 60’ ROW; ‘ROW90 contains the buffer for gravel roads, etc. The output grid data layer contained the feature grid cells coded with 2’s and the buffer with 1’s. The Reclassify module was used to reclassify cells with both data values into the actual ROW width value. See the “Re-Classifying Grid Data Values” section in Chapter 2 for more information related to the Reclassify module. This resulted with three buffer layers, one for each ROW width. The data values on the layers represented the actual ROW width, i.e., 60, 90, or 120. A problem emerged in that the buffers had overlap between the three buffer layers. Figure 3-25 displays an example of this problem. 84 Figure 3-25. Example of the buffer “overlap” problem. This happens because the buffer command was applied on road segments. Where roads intersect or join at junctions, the buffers will overlap. The map view window in Figure 325 illustrates this problem. I used the Grid – Calculus/Grid Calculator to eliminate the overlaps and to add the three grid layers together to make a composite layer containing all three buffers. Figure 3-26 displays a zoomed in area of the map. 85 Figure 3-26. A zoomed in portion of the Grapeview Road Maintenance District ROW map. The zoomed in portion of the map allows you to see more detail. The line shapes data layer for the roads is overlain using red for the roads. The boundary for the district is the thin black line. The Grid – Tools/Grid Proximity Buffer Module The Grid – Tools/Grid Proximity Buffer module delineates a buffer from features using a ‘Buffer distance’ parameter value as the buffer width. Grid cells on the input ‘>>Source Grid’ containing valid data values are considered to be features. The buffer calculations only evaluate grid cells on the ‘>>Source Grid’ containing valid data values. Grid cells containing “no data” values are not features. The Grid Proximity Buffer module delineates a buffer around features on a ‘>>Source Grid’. The value entered for the ‘Buffer distance’ parameter defines the width of the buffer. A buffer will be delineated around isolated single feature grid cells, on either side of linear features, and around contiguous groups of grid cells on the ‘>>Source Grid’. The module outputs three grid data layers. The output layers include a ‘Distance Grid’, an ‘Allocation Grid’, and a ‘Buffer Grid’. 86 Grid cells on the ‘Distance Grid’ output contain values representing the distance, in map units, the cell is from the nearest feature grid cell. The output value type for the ‘Distance Grid’ is floating-point. The ‘Buffer Grid’ is a reclassification of the ‘Distance Grid’ output using a user specified equidistance to create a set of discrete distance buffers from source features. The buffer zones are coded with the maximum distance value of the corresponding buffer interval. The output value type for the ‘Buffer Grid’ is integer. The third output for the module is called the ‘Allocation Grid’. In this case, the buffer grid cells contain data values using the feature data value they are closest. The output value type for the ‘Allocation Grid’ is integer. Example 1 Non-point source pollution is an issue in Mason County as well as in surrounding counties. One of the recently formed study groups has requested a table displaying acres by land use within 200’ either side of Mason County streams. The Grid – Tools/Grid Proximity Buffer module can be used to create this output. Figure 3-27 displays on the left the full extent view of the Mason County stream network. The right side displays a zoomed in area. The name of this grid layer is ‘MasonStrmGT6GR’. Figure 3-27. The stream network for Mason County. Figure 3-28 displays the grid data layer for general land use in Mason County, full extent display on left and a zoomed in area on the right. The name of this grid layer is ‘MasonGenLU’. 87 Figure 3-28. The Mason County general land use grid data layer. The ‘Grids’ parameter. The formula “ifelse(eq(a,9),9,0)” checks cell values on the input layer to see if they equal the forestry land use class value of 9. If they do, a 9 is saved to the corresponding cell on the output grid data layer; if not true, then a 0 is saved. The grid data layer output is displayed in Figure 3-35. 92 Figure 3-35. The forest land use map for Mason County. The red areas in the map displayed in Figure 3-35 are the areas of forest land use. This output layer has been renamed ‘MasonForestLU’. The ‘MasonSlpDeg’ data layer contains continuous data for slope gradient. I will use the Grid – Tools/Reclassify Grid Values module to create an output grid data layer containing four slope classes: 0° to 9°, 9° to 22.5°, 22.5° to 36°, and greater than 36°. I select the “simple table” option and build the table displayed in Figure 3-36. Figure 3-36. The slope gradient reclassify table. 93 Notice the values that will be used for the four classes are 1 through 4. The settings for the Reclassify Grid Values module to create the new slope class grid data layer are displayed in Figure 3-37. Figure 3-37. The Grid – Tools/Reclassify Grid Values module settings for the new slope class grid data layer. The output grid data layer was renamed ‘MasonSlpClasses’. It is displayed in Figure 338. 94 Figure 3-38. The ‘MasonSlpClasses’ grid data layer. The attribute table linked to the ‘MasonTrans’ line shapes data layer identifies road surface in the “SURF” attribute with a text string. The Shapes – Tools/Search in attributes table module supports text searching. It will take two executions of the module, one to search for the text “DIRT” and a second one searching for “GRAVEL”. The results of the second search will be appended to the results of the first search. At the same time the records meeting the search in the attribute table are highlighted, the corresponding objects in the layer are also selected and highlighted in the map view window. Once the two road types are selected I will use the Shapes – Tools/New layer from selected shapes module to create a new line shapes data layer. When the layer exists I can convert it from its vector format to raster or grid. The Shapes – Tools/Search in attributes table module was executed twice. The first time to find all the roads that were “DIRT” and the second time those that are “GRAVEL”. The first time I set the ‘Method’ parameter to “New selection” (see Figure 3-39). 95 Figure 3-39. The ‘Query builder for shapes’ settings for the first execution of the module. For the second execution I changed the ‘Method’ parameter and chose “Add to current selection”. This means the results of the first execution are not lost and the results of the two executions are combined together. Figure 3-40. A portion of the ‘MasonTrans’ attribute table on the left and the selected roads on the right. The grayed records in the attribute table in the figure (Figure 3-40) are for roads meeting the two queries. This is a small portion of the larger table but you can see that the selected records have “GRAVEL” or “DIRT” for the “SURF” field. The display of the roads shapes data layer on the right displays selected roads in red and non-selected roads in green. The Shapes – Tools/New layer from selected shapes module is used to create a new line shapes data layer containing only the selected roads on the ‘MasonTrans’ layer. I rename the new layer output from the Shapes – Tools/New layer from selected shapes module ‘MasonDirtGraRds’. The vector to raster conversion of this layer to a new grid 96 data layer (to be named ‘MasonDirtGraRdsGr’) is done with the Grid – Gridding/Shapes to Grid module. A new grid data layer is needed that contains only those roads that are within the forest land use areas. The Grid – Calculus/Grid Calculator is used to create this new layer. The two input layers were ‘MasonForestLU’ (first in the list) and ‘MasonDirtGraRdsGr’ (second in the list). The equation used for the ‘Formula’ parameter was “ifelse(eq(a,9),b,0)”. The data value used on the ‘MasonForestLU’ grid data layer for forest land use is “9”. I rename the output ‘MasonForRds’. The Grid – Tools/Grid Proximity Buffer module is used to create a buffer 6000’ wide around all dirt and gravel roads that are in forest land use areas. Although the module outputs by default three buffer related grid layers, the only one I am interested in is the ‘Distance Grid’. This is the one where the data values stored in the grid cells are the distances the cells are from the nearest cell on the ‘Source Grid’. Figure 3-41 displays the settings page for the execution of the Grid – Tools/Grid Proximity Buffer module. Figure 3-41. The Grid – Tools/Grid Proximity Buffer module settings page. The buffer area produced by these settings overlaps into land use areas outside of the forestry land use class. Before I use the Grid – Tools/Reclassify Grid Values module to reclassify the ‘Distance Grid’ into the distance zones, I need to create a new grid data layer with buffer areas outside of the forestry land use erased. I used the Grid – Calculus/Grid Calculator module for this task. Now I can apply the Grid – Tools/Reclassify Grid Values module to produce a layer containing the needed distance zones. 97 The “simple” table used to recode the ‘Distance Grid’ is displayed in Figure 3-42. Figure 3-42. The “simple” table used to produce the distance zone grid data layer. The “simple” table method is used to reclassify the continuous data of the ‘Distance Grid’ grid data layer into five distance zones. I rename the output ‘MasonForRdDist’. I use the Geostatistics – Grids/Zonal Grid Statistics module to produce a table showing the area for each slope class within the distance zones. In the ‘Zone Grid Statistics’ settings window, I enter the grid data layer ‘MasonForRdDist’ for the ‘>>Zone Grid’ parameter. The ‘MasonSlpClasses’ grid data layer is chosen for the ‘>Categorical Grids’ parameter. The area for each of the four categories on the ‘MasonSlpClasses’ within each distance zone on the ‘MasonForRdDist’ layer will be calculated and output in a ‘Zonal Statistics’ table. 98 Figure 3-43. The output ‘Zonal Statistics [c/5.3777]’. The output ‘Zonal Statistics [c/5.3777]’ table is displayed in Figure 3-43. The “Acres” column was added to the table using the Table – Calculus/Table Calculator module. The equation used in the Table – Calculus/Table Calculator to create the “Acres” column was “c/5.3777”. The fields or columns are referred to by an alpha representation of their position going from left to right. The first column is referred to by ‘a’, the second by ‘b’, etc. Thus, this formula says to create a new field or column by dividing the value in column ‘c’ (the ‘Count’ column) by the value 5.3777. The ‘Count’ column identifies the number of grid cells for that land use class in the buffer. It takes 5.3777 cells to make an acre. It is more convenient to save the ‘Zonal Statistics’ table as a dBase file, open it with an older version of Excel spreadsheet (a version prior to Office 2007 as newer versions no longer support importing dBase files), and develop the two tables asked for by the study group. Alternatives to Excel exist, for example, Open Office Calc. Figure 3-44 displays the table for acres by distance zone by slope class and Figure 3-45 the acres by relative environmental impact categories. 99 Figure 3-44. Acres by distance zone by slope class. Figure 3-45. Acres by relative environmental impact categories. After reviewing the results, the study group has identified some weaknesses in their impact weighting scheme. They decide to continue using this approach and take advantage of the GIS tools SAGA provides. The Grid – Tools/Threshold Buffer Module The SAGA Grid – Tools/Threshold Buffer module is used to delineate land areas contiguous to point, linear features, or polygon boundaries based on decision criteria applied to descriptive values. Figure 3-46. Settings page for module Threshold Buffer. 100 The ‘>>Features Grid’ cell data values should be categorical data for point, line, or polygon features that will serve as the buffer location source. The grid data layer for module input must be part of the grid system chosen for the ‘Grid system’ parameter. This is the case for all of the mandatory and optional input grid data layer parameters for this module. The grid data layer chosen for the ‘>>Value Grid’ parameter provides numeric values that can be used as a buffer attribute decision criteria depending on the options selected. An optional ‘>Threshold Grid’ grid data layer can be chosen as an alternative to using the ‘Threshold’ settings parameter. This optional layer provides threshold numeric values that will be applied as criteria in combination with the ‘Value Grid’. If you use the ‘>Threshold Grid’ parameter, any value entered in the value field to the right of the ‘Threshold’ option will be ignored. If you use the default entry, “[not set]” for the ‘>Threshold Grid’, this parameter will be ignored and the value in the value field for the option ‘Threshold’ will be used. The Grid – Tools/Threshold Buffer module grid output will contain numeric data values, in contiguous grid cells, meeting the threshold approach selected. There are two inter-related options available on the Threshold Buffer settings page. They are the ‘Threshold’ and ‘Threshold Type’ settings. Threshold and Threshold Type The ‘Threshold’ parameter is for entering a number that represents a threshold criterion. When you are using a ‘>>Value Grid’ as the source for criteria values, and the ‘Threshold Type’ option is “Absolute”, the number you enter for the ‘Threshold’ parameter will be applied as a “less than” criteria value. Starting from a feature cell on the ‘Feature Grid’, the module will check adjacent cells on the ‘>>Value Grid’ to see if they have numeric values less than the ‘Threshold’ parameter. If the value meets the test, the cell is identified as a buffer cell. If there are contiguous cells that meet the test, the buffer expands. Once no contiguous cells meet the test, the buffer expansion stops. Note that contiguous means side-to-side or diagonal contact. When you are using a ‘>>Value Grid’ to provide the criteria values, and the ‘Threshold Type’ option is “Relative to cell value”, the number you enter for the ‘Threshold’ parameter will be used to define a numeric range of acceptable values based on ‘>>Value Grid’ values relative to the value on the ‘>>Value Grid’ for the corresponding ‘Feature Grid’ cell. As an example, a ‘Feature Grid’ cell has a value of 5 and the corresponding cell value on the ‘>>Value Grid’ has a value of 20 and the ‘Threshold’ entry is 6. If an adjacent ‘>>Value Grid’ cell has a value of 19 it will be identified as a buffer cell. This is because the acceptance range, based on the ‘Threshold’ entry of 6, is for ‘>>Value Grid’ cell values from 14 through 26. The ‘>>Feature Grid’ cell corresponding to the ‘Value Grid’ cell value plus or minus 6 defines the acceptance range. The range is not based on the 101 numeric value stored on the ‘>>Feature Grid’. If the ‘Value Grid’ cell had a value of 27 or 13, the test would fail and the cell would not be a buffer cell. The initial identified buffer cell will in turn spawn the same test to be applied for grid cells adjacent or contiguous to it. The buffer will expand until no contiguous cells meet the test. Example 1 This example involves grid data layers for the Lake Cushman, Washington area. It is a relatively small area in the Olympic Mountains on the Olympic Peninsula. The goal is to delineate a buffer area adjacent to the streams. The buffer will be based on slopes less than 10 degrees. The numeric value 10 will be entered for the ‘Threshold’ parameter and the ‘Threshold Type’ option of “Absolute” will be selected. Figure 3-47. Settings for the Threshold Buffer module for the first example. The settings page in Figure 3-47 shows that the streams data layer (‘LCchannelnet’) is chosen as the ‘>>Features Grid’ and the slopes layer (‘LCslopeDEG’) as the ‘>>Values Grid’. The ‘>Threshold Grid’ is not used. The ‘Threshold’ parameter value “10” will serve as a “less than” criteria and is applied to the ‘LCslopeDEG’ grid data layer chosen as the ‘>>Value Grid’. 102 Figure 3-48. Module output on the left and less than 10 degree slope map on the right. The threshold buffer output, on the left in Figure 3-48, displays the area meeting the buffer criteria in red. The streams show up in black. The map on the right displays less than 10 degree slopes in red. The narrow black line is the line shapes data layer for streams overlain. Comparing the two maps, you can see that areas of less than 10 degree slopes that were not adjacent to streams, do not show up on the buffer output. Also, it would appear that the stream area not adjacent to less than 10 degree slopes might be in a slightly steep sided valley. The DEM for this area verifies that this is the physical situation. Figure 3-49. Threshold Buffer module output for Example 1. Figure 3-49 displays, on the left, the output from execution of the Threshold Buffer. In addition, the slopes grid data layer, ‘LCslopeDEG’, is displayed on the right. The same zoomed in area is showing for each so that the actual grid cell data values can be viewed. 103 Discussion The ‘Buffer Grid’ output on the left (Figure 3-49) displays grid cells with values of 0 (gray), 1 (red), or 2 (black). Cells with a value of 2 identify features from the ‘Features Grid’ data layer. The data value 2 is used regardless of what data value the features are coded with on the grid data layer being used for the ‘Features Grid’. Cells with a value of 1 are cells meeting the threshold criterion for being a buffer. The cells with 0’s did not meet the criteria. Slopes, in degrees, are the data values for the data layer on the right in Figure 3-49. You can see on the ‘Buffer Grid’ that five cells have been identified as meeting the buffer criteria. Looking at the slope values for the five cells, you see that the numeric values for slope are less than 10. The block of four buffer cells in the center in the graphic on the left, have slope values of 7.81, 9.46, .61, and 1.60. The cell touching the upper left corner of the block of four, has a slope value of 6.31. The buffer could not expand beyond these five cells because all of the contiguous cells contained slope values of greater than 10. Example 2 This example uses the streams data layer (‘LCchannelnet’) as the ‘>>Features Grid’ and the elevation layer (‘LCdem30’) as the ‘>>Values Grid’. The ‘>Threshold Grid’ is not used. The objective is to produce a map showing areas contiguous to streams that differ by less than 5 feet in elevation from adjacent grid cells. Thus, the buffer will include contiguous cells with elevation differences less than 5 feet. The buffer ends when an elevation difference of greater than 5 feet is contiguous. The numeric value 5 will be entered for the ‘Threshold’ parameter; the ‘Threshold Type’ option of “Relative to cell value” will be selected. A second execution of the module will use the same parameters except the ‘Threshold’ parameter will be changed to 10. This will allow a comparison of the two outputs. 104 Figure 3-50. Settings for the Threshold Buffer module using the “Relative from cell value” option. The ‘LCchannelnet’ grid data layer contains stream features and is selected for the ‘>>Features Grid’ input parameter. Elevations are contained on the ‘LCdem30’ grid. It is selected for the ‘>>Value Grid’. The output ‘>Features Grid’. The DEM for the area will provide the elevations. The ‘LCchannelnet’ grid data layer identifies both features (by their location, not by data value) and their flood stage values. Because a grid data layer has been chosen for the ‘>Threshold Grid’ parameter, the ‘Threshold’ parameter will not be used for this execution. 108 Figure 3-54. A zoomed in areas of the ‘LCchannelnet’ layer (upper left), the output ‘Buffer’ layer (upper right), and the ‘Value grid’ DEM layer. Output from the execution of the Threshold Buffer using a ‘Threshold grid’ and a ‘Threshold type’ of “Relative from cell value” is displayed in the upper right in Figure 354. We can see the values for the zoomed in portions of the three grid data layers in Figure 354. The threshold grid displays the data value ‘3’ for a flood stage potential for the seven cells of the stream visible. The buffer grid shows that three cells adjacent to the stream will flood. The “Relative from cell value” option sets an elevation range that serves as the “flood” cell criteria. In this case, we are dealing with a single data value, 3. The elevations for the seven stream cells are 1546, 1554, 1559, 1536, 1514, 1508, and 1481 going from left to right. Based on a plus or minus 3 range for each of these elevations, we see that the range criteria for each cell, from left to right, is 1543-1549, 1551-1557, 1556-1562, 1533-1540, 1511-1517, 1505-1511, and 1378-1484. This means that any cell, adjacent to any of the stream cells, that have an elevation within the range criteria for the stream cell, will be identified as a “flood” cell. 109 The third cell from the left, having an elevation value of 1559, is adjacent to cells that meet the range criteria, having elevation values of 1559 and 1561. The third river cell from the right, having an elevation value of 1514, is adjacent to a cell within the range criteria. The cell has an elevation value of 1513. The range criteria are also applied to cells adjacent to the ones meeting the river cell elevation range criteria. There are no contiguous cells to the “flood” cells that meet the range criteria, therefore, the “flood” area does not expand. 110 Chapter 4 – Creating Data Layers from DEMs in SAGA Introduction This chapter is made up of three sections. The sections are listed below along with the names of any modules used in that section. Creating Slope and Aspect Maps from a DEM Terrain Analysis – Morphometry/Local Morphometry Grid – Calculus/Grid Calculator Developing Terrain From Grid Data Layers Terrain Analysis – Morphometry/Surface Specific Points Grid – Calculus/Grid Calculator Creating Contour Shapes Data Layers Shapes – Grid/Contour Lines from Grid Creating Slope and Aspect Maps From a DEM The SAGA Terrain Analysis - Morphometry/Local Morphometry module can be used for creating slope and aspect grid data layers using an input digital elevation model grid data layer. When I choose the Local Morphometry module, the parameter settings page in Figure 4-1 appears. Figure 4-1. The Local Morphometry module parameters page. The entry for the ‘Grid system’ parameter must be the grid system the input grid data layer for the ‘>>Elevation’ parameter is a part. In order to be chosen, the grid system must be loaded in the current work session. When you click with the mouse pointer in the 111 value field to the right of the ‘Grid system’ parameter, a list of loaded grid systems for the current work session displays. Moving the mouse pointer over a list entry highlights the entry. Clicking the mouse button while the entry is highlighted chooses the grid system. The default for this parameter is “[not set]”. As you can tell by the familiar ‘>>’ symbols, only one input is needed and it should be an elevation grid data layer. When I click in the value field to the right of the ‘>>Elevation’ parameter, a small triangle appears followed by a pop-up list of loaded grid data layers for the grid system chosen for the ‘Grid system’ parameter. I choose the DEM grid data layer ‘MCdem30’ that contains elevations. If the DEM grid data layer had not been loaded, the Grid: Load Grid command in the Menu Bar ‘File’ drop down menu can be used to load it. The rest of this parameter section lists a set of output grid data layers that can be created from the DEM grid data layer. The first two in the list are ‘Slope’ and ‘Aspect’. They are mandatory outputs. You can see that in the value fields to the right of these labels the default entry is “[create]”. I will leave the defaults as the entries. The last three outputs in the list are set to the default “[not set]”. Whether these are chosen for output depends on whether further terrain analysis is planned that requires them as input. I will require them later so I am going to change the default entry from “[not set]” to “[create]”. In the ‘Options’ section, the ‘Method’ parameter allows me to choose a method from a pop-up list of seven for performing the morphometric analysis. Note that the method chosen will affect slope and aspect calculations as well as the three curvature outputs: curvature, plan curvature, and profile curvature. The seven methods are: Maximum Slope (Travis et al. 1975) Maximum Triangle Slope (Tarboton 1997) Least Squares Fit Plane (Costa-Cabral & Burgess 1996) Fit 2.Degree Polynom (Bauer, Rohdenburg, Bork 1985) Fit 2.Degree Polynom (Heerdegen & Beran 1982) Fit 2.Degree Polynom (Zevenbergen & Thorne 1987) Fit 3.Degree Polynom (Haralick 1983) SAGA defaults to “Fit 2.Degree Polynom (Zevenbergen & Thorne 1987)”. The first two approaches relate to flow routing algorithms as well as morphometrical analysis. Since their emphasis is on flow routing, they do not generate local morphometry using a bi-dimensional function and may be less suitable for calculating curvature outputs (curvature, plan curvature, and profile curvature). The “Maximum Slope” method, for example, calculates slope and aspect to the lowest neighboring cell to which the steepest 112 gradient is detected. This equals the single-flow direction method D8 and should be used to calculate channel slope. I will use the default method for calculating slope and aspect. Figure 4-2 shows the parameter settings I will use. Figure 4-2. The parameter settings page for creating Slope and Aspect grid data layers with the Local Morphometry module. I execute the module by clicking the ‘Okay’ button. The output grid data layers will be named with default data layer names: Slope, Aspect, Curvature, Plan Curvature, and Profile Curvature. When I click on the ‘Data’ tab at the bottom of the Workspace window, I see the new data layers at the bottom of the list of grid data layers for the grid system. The new data layers will display as thumbnails when I click on the ‘Thumbnails’ tab beneath the ‘Data’ tab. I need to take a close look at the data. I open the map view windows for the slope and aspect grid data layers by double-clicking on their names in the list. I follow the SAGA process for creating map view windows for the layers. Figure 4-3 displays these two windows. 113 Figure 4-3. Slope and Aspect grid data layers created with the Local Morphometry module. I highlight and make active the ‘Slope’ grid data by clicking on the layer name in the ‘Data’ tab area list in the Workspace window. When I move the mouse pointer across the ‘Slope’ map view window, the grid cell slope values are displayed along the bottom of the application window in the “Z” field. I do the same cursory data value check for the ‘Aspect’ grid data layer; first highlighting the ‘Aspect’ grid data layer name in the ‘Data’ list. The range of values for each data layer appears correct. It is important to make sure that the grid data layer name you want to view data values for is the one you make active because, regardless of which grid data layer is displayed in a map view, the one active in the ‘Data’ tab area list is the one whose data values will be displayed in the “Z” field. Reviewing the histograms for each grid data layer is another check. How to create a histogram is explained in Volume 1, Chapter 4. Figure 4-4. Reviewing the data histograms for the ‘Slope’ and ‘Aspect’ grid data layers. The histograms look okay. The value range for slope is from 0 to 67.055 and aspect from 0 to 360. When I check the tabular values for the histograms (by converting the graphic histogram to a table using the Menu Bar Histogram ‘Convert To Table’ command) for the range of data values for each of the layers, however, the values are not in the same range 114 as the values displayed in the grid data layer view windows or in their graphic histograms (see Figure 4-5). Figure 4-5. Reviewing the data tables for the ‘Slope’ and ‘Aspect’ grid data layers. The data value range for the tabular display of slope data is from 0 to 1.170338 and for aspect from 0 to 6.285185. I now check the parameter settings pages for the two layers. Making sure the ‘Object Properties’ window is displayed and the ‘Settings’ tab area showing, I click on the grid data layer name for one of the layers in the ‘Data’ tab area list of grid data layers. The parameters for the layer are displayed in the ‘Object Properties’ window. I do the same for the second data layer. The upper portion of the parameter settings pages for the two grid data layers is displayed in Figure 4-6. Toward the top of the ‘Parameters’ page, in the ‘General’ section, is a parameter called ‘Z-Factor’. I can see, for both data layers, the value field to the right of the ‘Z-Factor’ contains “57.2958” as an entry. Figure 4-6. The Z-Factor parameters for the new Slope and Aspect grid data layers. The actual data stored for these two data layers is radians and not degrees as displayed in the view windows. The values I see displayed when I move the mouse pointer across the grid cells are the data values in the grid cells (radians) multiplied by the Z-Factor 57.2958; the conversion factor for converting radians to angular degree value. This means that the actual grid cell data values are in radians and not in degrees. 115 I believe there will be times when some of the spatial analysis tasks I am interested will require radian input for slope and/or aspect and other times when slope and/or aspect will be required in degrees. Therefore, I am going to save versions of both, i.e., radians and degrees. First, I will save the new data layers with their radian data. I edit the ‘Name’ parameter in the ‘Object Properties’ for each of the two grid data layers. I name the slope grid data layer ‘MC1SlopeRAD’ and the ‘Aspect’ grid data layer ‘MC1AspectRAD’. Now I save the two grid data layers. Up until now, they have been temporary, un-saved working files. I right-click with the mouse on the ‘MC1SlopeRAD’ data layer name in the ‘Data’ tab area of the Workspace window. When the pop-up list of options appears, I click on ‘Save Grid As …’. When the ‘Save Grid’ dialog window appears I browse to the folder where I am saving data layers. I enter the same name I used in the ‘Object Properties’ window (‘MC1SlopeRAD’) and click the ‘Save’ button. I use the same naming approach for the ‘Aspect’ grid data layer. I am going to use the Grid-Calculus/Grid Calculator module to convert the slope and aspect values stored in the grid data layers from radians to degrees. I choose the Grid Calculator module. The ‘Grid Calculator’ settings page (Figure 4-7) displays. Figure 4-7. The initial Grid Calculator settings page. When I click in the value field to the right of the ‘Grid system’ label, a small triangle appears in the field and a pop-up list of grid systems loaded in this work session displays. I choose the grid system that the slope and aspect grid data layers are a part. Clicking in the value field to the right of the ‘>>Grids’ label causes an ellipsis symbol to appear in the field. When I click on the ellipsis symbol, a pop-up list of grid data layers available for the chosen grid system displays. On the list, I choose the ‘MC1slopeRAD’ grid data layer. In the formula value field, this data layer will be referred to as “a” since it is the first grid data layer in the list (as well as the only one). In the formula data field I 116 enter the equation “a*57.2958” (the value could also be highlighted and copied with CTRL-C from the settings tab area of the ‘Object Properties’ for the layer and pasted here with CTRL-V). This means that each of the grid cell values in the grid data layer will be multiplied by the value “57.2958” with the result being stored in the corresponding cell of the new output grid data layer. Figure 4-8. The Grid Calculator settings page with a formula for converting Radians to Degrees. Figure 4-8 shows the settings page before I click the ‘Okay’ button. SAGA names the new data layer using the equation entered for the ‘Formula’ parameter ‘a*57.2958’ by default. I use the re-naming process described earlier to enter a name that includes the text “DEG”. I now have two slope maps; one with data in radians and one with values in degrees (which, of course, has to be saved to disk storage). I follow the same procedures for creating an aspect data layer in degrees. Before exiting this work session, I go through a renaming process for the other three grid data layers created by the Local Morphometry module. Developing Terrain Form Grid Data Layers The module Terrain Analysis - Morphometry /Surface Specific Points can be used to produce a landform grid data layer. The landform of a specific cell is based on its elevation and/or slope compared to the elevation and/or slope of the 8 adjacent cells, and depending on the computation method chosen. When I execute this module, the ‘Surface Specific Points’ settings page displays (Figure 4-9). 117 Figure 4-9. The parameter settings page for the Surface Specific Points module. I choose an entry for ‘Grid system’ by clicking with the mouse pointer in the value field to the right of the ‘Grid system’ parameter. A list of grid systems loaded for the work session displays. I choose the grid system that the DEM grid data layer (‘MCdem30’) is a part. This is the DEM that will be chosen as input for the ‘>>Elevation’ parameter. A single input is required for elevation. I click in the value field to the right of the ‘>>Elevation’ parameter and a small triangle appears in the field along with a pop-up list of loaded grid data layers for the grid system. I choose the ‘MCdem30’ DEM grid data layer from the list. By default, the module output, the ‘Grids’ label causes an ellipsis symbol to appear in the field. When I click on the ellipsis symbol, a pop-up list of grid data layers available for the chosen grid system displays. I move the mouse pointer over the ‘MC1AllTer’ grid data layer name on the list and click the left mouse button. This selects the layer. I move the layer to the right side of the window by clicking on the button that displays the ‘>>’ symbols in the center of the window. I can move layers from the right side back to the left side by selecting and using the ‘DEM’ parameter a small triangle displays along with a pop-up list of grid data layers for the grid system chosen for the ‘Grid system’ parameter. I choose the ‘MCdem30’ data layer on the pop-up list of grid cell data layers. The default “[not set]” in the value field to the right of the ‘>Grids’ parameter. An ellipsis appears on the right in the field. When I click on the ellipsis with the mouse pointer, a pop-up list of current grid data layers appears. I choose the two DEM grid data layers from the list. It does not make any difference what order I select them in. But remember that when you develop an equation, the data layers are referred to by “a“, “b“, “c“, etc., in the order they appear in the list. In this example, the first one in the list is ‘MCdem30-0-sinks’ and the second one is ‘MCdem30’. So ‘MCdem30-0-sinks’ will be referred to in the formula as “a“ and ‘MCdem30’ as “b“. I am going to leave “[create]“ as the entry for the value field to the right of the ‘Elevation’ input is mandatory. The second one, ‘>Sink Routes’, would be used if I had pre-processed the DEM using the “Sink Drainage Route Detection” option rather than the “Sink Removal” approach. In the “Preparing a DEM Data Layer for Hydrologic Analysis” section I used the “Sink Removal” approach. Therefore, I do not need to use the ‘>Sink Routes’ input parameter here. The third input parameter, ‘>Weight’, is used to adjust the contribution each cell makes when calculating catchment areas. You would use the ‘Weight’ grid data layer if you 136 desire to reduce the area of each cell to the portion actually contributing to runoff. This approach could be used to account for evaporation and deep drainage. The weights are in the range 0-1, with 0 contributing no runoff and 1 contributing all precipitation falling on this cell. Runoff from a specific precipitation event can then be calculated by multiplying the output catchment area grid with the amount of precipitation. I am going to leave it as “[not set]“. Using “[not set]“, each cell will contribute to the catchment area with its’ full spatial area which equals the square of its’ cell size. Here are the parameters I enter. Clicking in the value field to the right of the ‘Grid system’ label, a pop-up list of grid systems available in this work session displays. I choose the grid system that the grid data layer I am going to use as input is a part. The DEM grid data layer is chosen by clicking in the value field to the right of the ‘>>Elevation’ label. I choose from the pop-up list of grid cell data layers the one that has been pre-processed for filling sinks, ‘MCdem30-0-Sinks’. I am going to leave the outputs set at the default. The ‘Catchment Area’ value field has a default of “[create]”. There are several other outputs that can be created. For now, I am going to leave them in their default setting, “[not set]”. This means they will not be created. Note that the default setting for the ‘Method’ parameter is “Multiple Flow Direction”. At a later time, when I have a better understanding of the different approaches, I am going to experiment with the various methods and compare the outputs. The settings I am entering for this module are displayed in Figure 5-9. 137 Figure 5-9. The parameter settings used to generate the ‘Catchment Area’ grid data layer. Last, I click on the ‘Okay’ button. Figure 5-10. The Catchment Area grid data layer. The output catchment area grid data layer in Figure 5-10 does not look, visually, very complex. I can check the data value range by opening the ‘Object Properties’ and displaying the ‘Description’ tab area for the new layer. I highlight and make active the 138 catchment area data layer name in the ‘Data’ tab list of the Workspace window by clicking on it once. The description window for the catchment area grid data layer is displayed in Figure 511. Figure 5-11. The catchment area grid data description window. In this example, the range of data values goes from a minimum of 10,889.966 to a maximum of 594,853,824; a substantial range. The huge range of the data is the reason for the lack of detail in the view window for the grid data layer. The default color classification palette does not provide enough variation. Using the “logarithmic scaling” mode will improve the level of detail. The ‘Settings’ tab area of the ‘Object Properties’ window in Figure 5-12 has been scrolled down part way. 139 Figure 5-12. The catchment area grid data description window. I can see the data value range displayed to the right of the ‘Value Range’ label. Underneath this parameter is another one called ‘Mode’. There are three options supported: Linear (default), Logarithmic (up), and Logarithmic (down). I choose the “Logarithmic (up)” option. After some experimentation, I find that entering 100000 in the ‘Logarithmic Stretch Factor’ value field yields better visual interpretation. Now I am able to get a better impression of the results of generating catchment areas (Figure 5-13). 140 Figure 5-13. The Catchment Area grid data layer using the ‘Logarithmic’ option. The ‘Linear Flow Threshold’ parameter (see Figure 5-9) was not discussed. This threshold can be used to set a catchment area above which the algorithm used switches from multiple flow direction to single flow direction. This threshold would usually be set to a value, above which a channel network is developing (the “blue lines” on a topographic map). Using the D8 also on the hill-slope would usually result in unrealistic parallel flow paths. But along channels, it allows the flow to concentrate. So the combination of both can be used to make use of the advantages of both algorithms. I rename this data layer ‘MCcatchmentArea’. Defining Channel Networks Channel or stream network data layers, grid and vector, are generated using the Terrain Analysis - Channels/Channel Network module. The module settings page is displayed in Figure 5-14. 141 Figure 5-14. The Channel Network parameter settings page. The grid system is chosen by clicking with the mouse pointer in the value field to the right of the ‘Grid system’ label. A list of loaded grid systems for the work session displays. I choose the grid system the input grid layers are a part. You can see from the ‘>>’ symbols on the ‘Channel Network’ settings page that two grid data layers are mandatory input: ‘>>Elevation’ and ‘>>Initiation Grid’. The elevation grid is used to trace the channels. The digital elevation model (DEM) with the sinks filled (‘MCdem30-0-sinks’) will be used for this parameter. The initiation grid supplies information about the starting points of the channels. These cells are referred to as “initiation” cells. The ‘Catchment Area’ grid data layer can usually be used as the ‘Initiation Grid’. I will use it in this example. The ‘Catchment Area’ grid data layer (‘MCcatchmentArea’) was generated and described in the “Calculating Catchment Area” earlier in this chapter. When I click in the value field to the right of the ‘>>Elevation’ parameter, a small triangle appears in the field and a pop-up list of grid data layers available for the chosen grid system displays. I choose the DEM data layer with the sinks filled, ‘MCdem30-0sinks’. The same list of grid data layers displays when I click in the value field to the right of the ‘>>Initiation Grid’ parameter. I choose the catchment area data layer, ‘MCcatchmentArea’, to use for the ‘>>Initiation Grid’ parameter. Note that immediately below the ‘>>Initiation Grid’ label, in Figure 5-14, are two options. The first one is for ‘Initiation Type’. This option defines the condition under 142 which a channel is started. There are three choices: Less Than, Equal To, Greater Than. These choices relate to the threshold value set in the next row labeled ‘Initiation Threshold’. In this example, I will set the threshold to 1000000 and the ‘Initiation Type’ to “Greater Than”. The higher you set the threshold, the greater will be the distance from the divide to the initiation point, i.e., the more area the catchment area needs to be to supply enough water to form a channel. The parameters page shows that there are three output data layers. ‘Channel Network’ and ‘Channel Direction’ are both grid data layers. The third one is a line shapes data layer for the channel network. I will use the default “[create]” that appears in the value fields. The ‘Options’ section has one variable. It is used to control the length (in cells) of the shortest channel segments. I will limit the creation of channel segments to those greater than 10 cells in length. I enter ‘10’ in the value field to the right of the label ‘Min. Segment Length’. The parameter settings I will use for this example are displayed in Figure 5-15. Figure 5-15. The parameter settings used for the Channel Network module. Now I execute the module by clicking on the ‘Okay’ button. 143 Figure 5-16. The Channel network and Channel Direction grid data layers. The channel network grid data layer (on the left in Figure 5-16) contains cell values identifying the stream order after Shreve (1966) of the segment that flows across the cell. The zoomed in portion of the channel network grid data layer (Figure 5-17) shows cells having a stream order value of 33. Figure 5-17. A zoomed in portion of the Channel Network grid data layer. The channel direction grid data layer (on the right in Figure 5-16) shows the flow direction of each cell. The directions are coded from one through eight, counting clockwise from northeast (1) around to north (8). The grid data value in the cell (see Figure 5-18) identifies the direction out of the cell for the flow. A value of 3.00 means the flow is from the SE corner. A 4.00 value identifies the flow moving from the bottom of the cell. 144 Figure 5-18. A zoomed in portion of the Channel Direction grid data layer. The third output element created by this module is a shapes data layer. It has been added as a line shapes data layer to the list of data layers in the ‘Data’ tab area of the Workspace. I can see on the list of shape data layers that there is one called ‘Channel Network’. It captures the same channel network as the grid data layer channel network except in a vector format. If the shapes data layer is going to be used as an overlay for grid data layers, the color used for the stream segments may have to be changed depending on the color palette being used for the grid data layer it would overlay. I can adjust the color. First, I must make sure that the shapes data layer is selected in the ‘Data’ list in the Workspace window. Next, the ‘Object Properties’ window must be displayed in the work area. If it is not already displayed, I can activate it by selecting ‘Show Object Properties’ from the ‘Window’ title on the Menu Bar. The parameter page in the ‘Object Properties’ window will be populated with the settings for the data layer. You can see that the ‘Object Properties’ parameters page in Figure 519 is for the shapes data layer for the ‘Channel Network’. 145 Figure 5-19. The parameter settings page for the ‘Channel Network’ shapes data layer. There are several ways to adjust the color of the stream segments. Since my intent is not to portray attributes by colors but only to choose a suitable color for displaying the segments as a map overlay, I will use a single color for all. In the ‘Display: Color Classification’ section, there is a parameter called ‘Unique Symbol’. The first parameter under it is ‘Color’ and that is where the default color for the stream segments in the ‘Channel Network’ shapes data layer has been specified. The default is lime green. This is the parameter I will change. When I click in the value field 146 to the right of the ‘Color’ label, a list of colors, like the one in Figure 5-20, is displayed. I click on the color I want to use to display the channel network stream segments. Figure 5-20. A color table used with shapes data layers. I choose white. Since the default background for the shapes view window is also white, the segments will seem to disappear. However, this is just because I am using white for the segments against a white background. When I overlay the white segments on a DEM grid data layer the channel network will stand out nicely. Figure 5-21 shows an example of the channel network shapes data layer as an overlay using white on the ‘MC1dem30sinks’ grid data layer. 147 Figure 5-21. The DEM grid data layer with a shape data layer channel network overlain. I rename both the grid and line shapes data layers for the channel network to ‘MCchannelNet’. Defining Watershed Basins Each segment of a channel network grid data layer is associated with a watershed basin. This watershed basin corresponds to the catchment area of the lower point of the segment minus the basins associated to other segments located upslope. The watershed basin is delineated from the junction point at the end of the segment. This point is coded with a –1 while the segment is coded with the stream order. Knowing that junction points are coded with –1’s suggests another approach for calculating watershed basins. An input grid data layer consisting of only the –1 cells with all other cells coded as “NoData” cells could be used as a pour point grid instead of a channel network grid. The SAGA module for creating a watershed basin grid data layer is Terrain Analysis Channels/Watershed Basins. Figure 5-22 displays the ‘ Watershed Basins’ settings page. 148 Figure 5-22. The Watershed Basins module parameter settings page. Two grid data layers are required for input. These are a DEM for elevation with the sinks filled (unless you specify a ‘Sink Route’ grid) and one for channel network. The ‘Grid system’ is chosen by first clicking with the mouse pointer in the value field to the right of the ‘Grid system’ parameter. A small triangle appears in the field and a list of the grid systems loaded for the work session displays. There are a couple of grid systems loaded in this work session. I select the one that the grid data layers that will be chosen for the ‘>>Elevation’ and ‘>>Channel Network’ parameters are a part. The grid data layers to be used for input are chosen by clicking in the value fields to the right of the ‘>>Elevation’ and ‘>>Channel Network’ parameters. A list of grid data layers that are part of the chosen grid system displays. For the ‘Elevation’ layer I choose the ‘MCdem30-0-sinks’ DEM that has been pre-processed by filling sinks. I choose the ‘MCchannelNet’ grid data layer for the ‘Channel Network’ input. I want to create a grid data layer for watershed basins so the default “[create]“ for the ‘Layer A’ and ‘>>Layer B’ parameters are for choosing two shape data layers for input. In the ‘Options’ section for the ‘Method’ parameter, I chose the “Complete Intersection” method. The output shapes data layer from these settings is displayed in Figure 6-4. The output consists of a number of polygon objects. The ShapesPolygon/Polygon Union module was then used to combine the polygon objects into a single object in a data layer named ‘MasonGV’. Figure 6-4. The polygons representing the hamlet of Grapeview (‘MasonGV’). The boundaries for the two districts are overlain on the output in Figure 6-4. The last step is to create a grid data layer version of the shapes data layer. The Grid-Gridding/Shapes to Grid module is used for this operation. Figure 6-5. The Grid-Gridding/Shapes to Grid module settings, vector to raster conversion. 155 The polygon shapes data layer (‘MasonGV’) defines Grapeview with a single polygon object. The ‘GVID’ attribute contains a ‘1’ data value for the polygon. I chose this attribute to provide the data value for the output on the grid data layer. This means the grid cells on the output grid data layer falling within the polygon on the input shapes data layer will be coded with the value of ‘1’. The output grid data layer from this module execution was renamed ‘MasonGVgrid’. Criterion 2: identifying Open Space land use in Mason County and Grapeview A county land use grid data layer exists (‘MasonGenLU’). The open space land use class on the county grid data layer that exists in the Grapeview hamlet will be copied into a grid data layer for the hamlet. The Grid-Calculus/Grid Calculator module is used to build this open space land use layer for the hamlet. The Mason County General Land Use map grid data layer is displayed in Figure 6-6. Figure 6-6. The Mason County General Land Use class map. The majority of the land area is in the forestry land use class. At this scale, it is difficult to see many areas in the open space class (displayed in pink). The data value for open space is “7”. The large white area in the northwest quadrant is federally owned land and includes portions of the Olympic National Park and Olympic National Forest. The county does not determine land use for federally managed lands. Figure 6-7 displays the Grid Calculator settings for creating a grid data layer containing the open space land areas for the Grapeview hamlet. 156 Figure 6-7. The Grid Calculator module settings for creating the open space land use grid data layer for the Grapeview hamlet. The input grid data layers are ‘MasonGVgrid’ (referred to in the formula as variable “a”) and ‘MasonGenLU’ (referred to in the formula as variable “b”). The formula applied is: ifelse(eq(a,1),ifelse(eq(b,7),b,0),0) The ‘MasonGVgrid’ layer uses a “1” for all the grid cells in the Grapeview hamlet. The first part of the formula checks layer “a” to see if the cell data value is “1”. If the value is a “1” (i.e. the condition is true), then the second part of the formula is checked. The data value for the open space land use on the ‘MasonGenLU’ layer is “7”. If this condition is true then the data value is saved to the corresponding cell on the output grid data layer. If the value is not found meaning the condition is false, a “0” is saved to the output layer grid cell. The output grid data layer is displayed in Figure 6-8. The grid data layer is renamed ‘GVOpSpgr’. 157 Figure 6-8. Open Space and Rural land uses in Grapeview (‘GVOpSpgr)’. The map display in Figure 6-8 represents the “population” of open space parcels in Grapeview that meet the first level of suitability for a community park. The parcels showing in the pink color have the open space land use classification. Criterion 3: Open Space parcels must be a minimum size of 5 contiguous acres The Grid-Filter/Filter Clumps module (see Figure 6-9) is used to identify groupings (larger than a specified minimum size) of grid cells having the same cell value. A value is entered to specify a minimum number of cells in a group. The value entered for specifying a minimum group size is the number of cells representing the minimum group size. In this situation, this means I must calculate how many cells make up 5 acres of land area. Each cell in the grid system is 90’ square or 8100 square feet. Five acres is about 217,800 square feet. Thus, 5 acres is approximately 27 grid cells. 158 Figure 6-9. The settings for the Grid-Filter/Filter Clumps module for identifying >5 acre cell groupings. The input grid data layer is ‘GVOpSpgr’. As noted earlier, this layer has open space land areas identified using the data value “7”. The minimum number of cells in a 5-acre group is calculated to be “27”. The output from executing these module settings is displayed in Figure 6-10. 159 Figure 6-10. Grapeview Open Space areas greater than 5 acres in size. Grid cells containing the open space data value “7” but are not part of a contiguous group of cells meeting the 5 acre threshold, are recoded by the Filter Clumps module. The module places ‘no data’ values in these cells. Because these ‘no data’ grid cells are really valid data cells containing zeroes, I used the Grid-Tools/Reclassify Grid Values module to convert the values to zeroes. The name for the grid data layer stayed the same (‘GVOpSpgr’). You can see that small areas have dropped out when you compare Figures 6-8 and 6-10. The open space areas displayed in Figure 6-10 represent the total population of suitable land areas that are land use class open space and meet the minimum size criteria of 5 acres. I used the Shapes-Grid/Vectorising Grid Classes module to create a polygon shapes data layer that can be used in later steps of this evaluation process. I have named this polygon shapes data layer ‘GVOpSpPoly’. Figure 6-11 displays the shapes data layer. 160 Figure 6-11. Open Space areas in Grapeview greater than 5 acres in size. Figure 6-12 displays the attribute table linked to the ‘GVOpSpPoly’ shapes data layer. Figure 6-12. The attribute table for the ‘GVOpSpPoly’ shapes data layer. I will add columns/fields to this table as the evaluation progresses. The labels displayed for the polygons in Figure 6-11 are from the ‘Site’ column in the attribute table. The site numbers are the numbers I have assigned to each of the open space groups. SAGA does not have good cartographic output support; thus, the labels are displayed but not cartographically pleasing. 161 On the far right in the table is a column labeled “ACRES”. This is a generated column using the Table-Calculus/Table calculator for shapes module. The “AREA” field displays the square feet for the polygon object. The Table-Calculus/Table calculator for shapes module settings for generating the “ACRES” column are displayed in Figure 613. Figure 6-13. The Table-Calculus/Table calculator for shapes module settings for calculating “ACRES”. The polygon shapes data layer is chosen for the ‘>>Shapes’ parameter. The sixth column from the left is referred to in the equation as variable ‘f’. That column displays the square feet for each polygon. An acre is 43,560 square feet. The formula divides the polygon square footage by the square feet in an acre to calculate acres. Criterion 4: Open Space land parcels with road access There are several ways to approach this criterion. Since there are only 10 groupings of open space land use greater than 5 acres in size, you could display the county roads (‘MasonRoadsGraPav’) line shapes data layer overlain on the ‘GVOpSpgr’ grid data layer and visually determine the overlap. Another way is to use the Grid Calculator module to add the grid data layer version of the county roads to the ‘GVOpSpgr’ grid data layer. On the grid data layer version of the county roads, a data value of ‘2’ is for gravel surface roads and a value of ‘3’ is for pavement surface roads. Open space on the ‘GVOpSpgr’ data layer is coded as ‘7’. A mathematical overlay (adding the two together into a single layer) would mean that data values of ‘9’ and ‘10’ on the output grid data layer would represent areas where roads overlap open space land use. Figure 6-14 displays the Grid Calculator settings to produce this new grid data layer. I rename the new grid layer ‘GVOpSpRoads’. 162 Figure 6-14. Grid Calculator settings for overlaying gravel/paved roads with open space land use in Grapeview. Figure 6-15 displays the Grid Calculator output. Site 4 is the only site out of the 10 that does not have a road contact. Figure 6-15. Display of ‘GVOpSpRoads’ grid data layer. Site 4 does not have a road contact. It is relatively easy to zoom in on the 10 sites to assess whether any contain roads or have road contact. Criterion 5: The park parcel cannot be adjacent to these land uses: industry, commercial, mineral extraction, and utilities. Criterion 6: The park parcel may be adjacent to these land uses: agriculture, commercial recreation, public, rural and residential. 163 Criterion 7: The park is most suitable adjacent to these land uses: open space, forest, and water. These three criteria (5, 6, and 7) address compatibility with adjacent land uses. I am going to create a table showing the linear distance for adjacency, in feet, for land uses bordering each of the ten sites. This table can be used to evaluate these three criteria. The first step toward developing this table is to create a grid data layer containing a one cell wide zone around each of the sites. This border zone will be one cell width or 90’ (since a cell is 90’ on a side). The data values for the border zone cells for site #1 will have a data value of 100; site #2 a data value of 200; and so on. This data value numbering system retains a numeric link to the sites and is initiated with a new attribute called ‘BUFFNUM’. The Table-Calculus/Table calculator for shapes module was used to create the new attribute ‘BUFFNUM’ for the ‘GVOpSpPoly’ shapes data layer. The existing attribute, ‘SITE’ (numeric values 1 through 10), values were multiplied by 100 to create the values for the ‘BUFFNUM’ attribute field. The Shapes-Tools/Shapes Buffer module was used to develop a polygon shapes data layer using ‘GVOpSpPoly’ as input. Each polygon object in the ‘GVOpSpPoly’ layer will increase in size by the addition of a 90’ buffer added to it in the output layer. Figure 6-16 displays the settings for the execution of this module. Figure 6-16. Settings used for the Shapes-Tools/Shapes Buffer module. 164 The ‘>>Shapes’ parameter shows ‘GVOpSpPoly’ polygon shapes data layer as input. This layer contains polygon objects representing the ten open space sites that are greater than 5 acres in size. The ‘ID’ attribute identifies the site ID for each object. In this case, the attribute selected does not have any influence on the output as the attribute values are not used. The grid cell size is 90’. This is the desired width of the output buffer. A “90” is entered for the ‘Buffer Distance (Fixed)’ value. The objects on the output layer will be the original object dimensions plus the 90’ buffer added around each polygon object. This is an interim output. I can use the Shapes-Polygons/Polygon Intersection module to create a polygon shapes data layer containing only the buffer portion of the output polygons from the Shapes Buffer module. This buffer portion is the single-cell width border zone around each polygon. Figure 6-17 displays the settings for the Shapes-Polygons/Polygon Intersection module for this execution. Figure 6-17. Settings for the Shapes-Polygons/Polygon Intersection module execution. The buffer polygons are created by subtracting the original open space site polygon objects from the output objects of the execution of the Shapes Buffer module. The result is the difference between the two layer objects. This is the buffer that was added by the Shapes Buffer module around each open space site; a 90’ border zone. The output layer is called ‘GVOpSpEDGE’. Figure 6-18 displays an example of the objects on the ‘GVOpSpEDGE’ layer. 165 Figure 6-18. Zoomed in portions of the ‘GVOpSpEDGE’ output data layer. The next step is to convert the ‘GVOpSpEDGE’ output to a grid data layer so it can be used with the ‘MasonGenLU’ grid data layer to identify adjacent land uses for the sites. Unfortunately, the values for the ‘BUFFNUM’ attribute were lost in the intersection process. Before doing a vector to raster conversion I edited the attribute ‘BUFFNUM’ and added the values back in. The Grid-Gridding/Shapes to Grid module is used for the vector to raster conversion. Figure 6-19 displays the settings used. Figure 6-19. The settings used for the Grid-Gridding/Shapes to Grid vector to raster conversion. The input shapes data layer is the one displayed in Figure 6-18, ‘GVOpSpEDGE’. The attribute chosen to provide the data values for the grid cells on the output grid data layer is ‘BUFFNUM’. This attribute contains values that link the border zone polygons to the site polygon. The output grid data layer is re-named ‘GVOpSpEDGEgr’. 166 The next step is to use the ‘GVOpSpEDGEgr’ grid data layer with the grid layer containing land use classes for Mason County (‘MasonGenLU’). I will use the Grid Calculator module to add the two layers together. The data values on the ‘GVOpSpEDGEgr’ layer are from 100 to 1000. The data values on the ‘MasonGenLU’ layer are 1 through 12. Adding the two layers results with a set of new, unique, data values that can easily be related to the sites. For example, the data value of 100 from the ‘GVOpSpEDGEgr’ layer added to any of the land use classes would result with potential output values ranging from 101 through 112. The second and third digits in the value identify the land use class from the ‘MasonGenLU’ layer. Thus, the 100 series of data values relate to site #1 and border land use classes, the 400 series to site #4, etc. The coding of the ‘GVOpSpEDGEgr’ layer is crucial in this approach in order to distinguish the intersected classes. I did encounter one issue with this process. Site numbers 2 and 3, in one nearby area, are separated by a single grid cell. This was not enough separation for a buffer to be created for each site; instead, one buffer was created around both sites. This is not a big problem and was corrected by re-coding the border zone around site #2 on the grid data layer. Thus, the border zone around site #2 contained the value 200 and the zone around site #3 contained the value 300. This re-coding was accomplished using the Grid – Tools/Change Cell Values [interactive] module. The versatile Grid Calculator module is executed to create a new grid data layer using the ‘MasonGenLU’ and ‘GVOpSpEDGEgr’ grid data layers showing only the land use classes for the one-cell wide border zones. Figure 6-20 displays the settings for the Grid Calculator used for this objective. Figure 6-20. The Grid Calculator settings for creating a grid data layer for land use classes along the site borders. 167 The formula used is “ifelse(gt(a,0),(a+b),0)”. Variable “a” is the reference to the ‘GVOpSpEDGEgr’ layer and variable “b” refers to the ‘MasonGenLU’ layer. Layer “a” is checked to see if a found data value is greater than 0. If the condition is true, the data value would represent a border zone. If true, the data value for the cell in layer “a” (the border zone data value) is added to the corresponding cell data value in “b” (the land use class data value) and output to the corresponding cell on the output data layer. Again, the equation portion ‘(a+b)’ is adding the border zone number (from layer “a”) to the land use class data value (from layer “b”). The output in the grid data layer (‘GVOpSpEDGElanduse’) captures the land use adjacent to each of the open space sites being evaluated. The Geostatistics-Grids/Zonal Grid Statistics module can be used to produce a table that summarizes the number of cells for each data value that occurs on the output grid data layer (‘GVOpSpEDGElanduse’). Each of the data values includes the site number so the data can be summarized and displayed by site number. Figure 6-21 displays the settings for the execution of the Geostatistics-Grids/Zonal Grid Statistics module for the output ‘Result Table’. Figure 6-21. The ‘Zonal Grid Statistics’ module settings. I saved the output ‘Zonal Statistics’ table and opened it in Excel. I multiplied the cell counts by 90 to figure out the number of feet of adjacency for each land use by site. I had to do a visual interpretation for “water” adjacency as that is not a land use class included on the ‘MasonGenLU’ grid data layer. The final table is displayed in Figure 6-22. 168 Figure 6-22. The table for land use adjacency for each of the sites. The three criteria can now be evaluated. The above process could have been shortened by using the Geostatistics-Grids/Zonal Grid Statistics module at the beginning with the ‘GVOpSpEDGEgr’ layer for the ‘>>Zone Grid’ parameter and the ‘MasonGenLU’ layer for the ‘>Categorical Grids’ parameter. The lengthier approach described above illustrates several SAGA modules used together to achieve the objective. Criterion 5: The park parcel cannot be adjacent to these land uses: industry, commercial, mineral extraction, and utilities There are four land use classes that the committee feels a community park should not be adjacent. These are: industry, commercial, mineral extraction, and utilities. Viewing the table in Figure 6-22, site #5 has some type of utility facility adjacent to it. Criterion 6: The park parcel may be adjacent to these land uses: agriculture, commercial recreation, public, rural and residential There are five land use classes that the committee feels a community park is moderately compatible. These are: agriculture, commercial recreation, public, rural, and residential. At this time, this can be viewed on a presence or absence basis. All of the sites have residential land use adjacent to their boundaries. None of the sites are adjacent to a public land use. Site #4 borders some agriculture land use. Commercial recreation is adjacent to site #9. 169 Criterion 7: The park is most suitable adjacent to these land uses: open space, forest, rural, and water Open Space is the land use for each of the 10 sites. There are not any sites that will have adjacent open space. Any adjacent open space is part of the greater than 5-acre group. Sites 1, 2, 4, 6, 7, 8, and 9 all have adjacent lands in forestry land use. The Rural land use is found adjacent to site #’s 1, 2, 3, 6, 7, 9, and 10. The last one, water, is near site #’s 1, 4, and 5. A final step in evaluating criteria 5, 6, and 7 will be to calculate the adjacent land use lengths as a percentage of the total border length for the site. The sites can then be rated accordingly. Sites with higher percentages of border length in these land uses would be more acceptable. Criterion 8: The highest amount of ground surface with slopes between 0 and 3 degrees A table showing the acreage having 0 to 3 degree slopes will be generated. The percentage of acreage in 0 to 3 degree slope, relative to the total acreage in the site, will be calculated. The sites can be ranked accordingly. There is a grid data layer representing the 10 sites. The grid data values range from 1 through 10. It is named ‘GVOsites’. This layer will be used with a second layer that has two classes; less than 3 degree slopes (1’s) and slopes greater than 3 degree (‘0s’). The ‘MasonSlpDeg’ grid data layer contains slope data, in degrees, for the entire county. The Grid Calculator module will be used to create a grid data layer that contains less than 3 degree slopes only. The settings for the Grid Calculator are displayed in Figure 623 for creating the ‘MasonSlopeLT3’ grid data layer. Figure 6-23. Grid Calculator settings for creating a grid layer containing less than 3 degree slopes. 170 The ‘MasonSlpDeg’ grid data layer was chosen for the ‘>>Grids’ parameter. The equation checks the cell value on the input layer to see if it is less than 3 degrees. If it meets the condition (i.e., the condition is true), the value “1” is saved to the corresponding grid cell on the output grid data layer. Otherwise, if it does not meet the condition (it is false), a 0 value is output to the grid data layer. I name the output grid data layer ‘MasonSlopeLT3’. The Geostatistics-Grids/Zonal Grid Statistics module can be used with the ‘GVOsites’ layer chosen to represent zones and the ‘MasonSlopeLT3’ layer as a categorical grid. The settings for this module are displayed in Figure 6-24. Figure 6-24. The Zonal Grid Statistics module settings. The output from the Zonal Grid Statistics module is a table with the default name ‘Zonal Statistics’. I saved this table in a text file and then opened it with Excel. I used features in Excel to re-organize the ‘Zonal Statistics’ table output and to calculate percentages. Figure 6-25 displays my spreadsheet version of this table. 171 Figure 6-25. Acres and percent of total site acreage with slopes less than 3 degrees. Using the data in the last column of this table, it is relatively easy to rank the sites in order of preference related to less than 3 degree slope. The higher the percentage of less than 3 degree slope, the more suitable. The sites rank from best to worst as follows: 1, 8, 2, 6, 9, 3, 7, 5, 10, 4. Criterion 9: The viewing of these land use classes in the foreground and middleground will be minimized: industry, commercial, mineral extraction, utilities. This criterion will be used to compare parcel choices The viewshed, or area that can be seen, from each of the ten open space land use sites will be calculated. A centroid point in each site will serve as the observation point for the viewshed determination. The foreground, middle-ground, and background visual zones surrounding each site (observation point) will be delineated using distance definitions for the three zones. A summary table will be prepared listing the acres for industry, commercial, mineral extraction, and utilities land uses in the viewshed for each site and by the three distance zones. This analysis will be a worst-case analysis because vegetation will not be a factor. The delineated viewshed will be larger in size compared to the same viewshed delineated taking vegetation heights into consideration. Vegetation heights tend to shield or block views. There is not an easy way to select an observer point. An approach involving elevation, aspect, roads, etc., could be developed. Rather than develop a complex approach, I am going to use the Shapes-Polygons/Polygon Centroid module. This module will determine a centroid point or mean center for each site. This point is calculated using the site x and y boundary coordinates. One advantage of using a centroid is that the identification procedure for the centroid will be consistent for all sites. Figure 6-26 displays the settings used for this execution. A note of caution. Using the centroid approach with these polygons produced acceptable observation points. However, some polygon shapes could result with a centroid outside of the polygon. 172 Figure 6-26. The settings for the Shapes-Polygons/Polygon Centroid module. The single input is the polygon shapes data layer with the sites as polygon objects. This is layer ‘GVOpSpPoly’. The output point shapes data layer containing the centroid points will be renamed ‘GVOpSpObserve’. Figure 6-27 displays a zoomed in portion of the output for the centroid in site number 9. The ‘GVOpSpPoly’ shapes data layer is overlain on the point layer in the figure. Figure 6-27. A zoomed in view of the centroid (red dot) in site #9. The following discussion will focus on how a single site is processed. The procedures described were applied to all of the sites. 173 The Terrain Analysis-Lighting, Visibility/Visibility (single point) [interactive] module is used to calculate areas that can be seen from an observation point. A digital elevation model (DEM) provides the model for the topography. The name of the existing DEM for the county is ‘MasonDEM’. The DEM is displayed in Figure 6-28 with the ‘GVOpSpObserve’ point and ‘GVOpSpPoly’ polygon shapes data layers overlain. Figure 6-28. The observer point shapes data layer displayed with the DEM grid data layer. The settings for executing this interactive module are displayed in Figure 6-29. The interactive part of this module is collecting the x and y coordinate for the observation point. I will use the ‘GVOpSpObserve’ point shapes data layer as the source for the observation point coordinate. As noted above, Figure 6-28 displays the point shapes data layer overlain on the ‘MasonDEM’ grid data layer. This map view window can be displayed before executing the Terrain Analysis-Lighting, Visibility/Visibility (single point) [interactive] module or after. The on-screen point digitizing described below uses this map view window. 174 Figure 6-29. The settings for the Terrain Analysis-Lighting, Visibility/Visibility (single point) [interactive] module. The ‘MasonDEM’ grid data layer has been chosen for the ‘>>Elevation’ parameter. The ‘Height’ parameter is for specifying the height of the observation point above the terrain. I am placing this observer 5’ above the terrain. The ‘Unit’ parameter is for choosing the method. The “Visibility” method is one of four options supported. When the ‘Okay’ button is clicked on, the module begins execution. After execution is started, before the module can start delineating the viewshed or seen area for the observation point, you must digitize, on-screen, the location of the ) on the Tool Bar is used for digitizing this observation point. The ‘Action’ tool ( point. Once the ‘Action’ tool is chosen, the mouse pointer turns into a plus symbol with a lower-case ‘i’ at the bottom right. I am going to move the ‘Action’ tool pointer over the observation point that is the source of the viewshed for this execution. In order to make it easier to “digitize” the observer point location, I zoom in on the point in the map view window displayed in Figure 6-28. With the ‘Action’ tool positioned on the zoomed in point, I click the left button of the mouse. The x and y coordinate pair for the point is collected and passed to the module. Note that any grid data layer can be used as the coordinate source for the digitized point. It is best to use a layer that will provide detail if you zoom in on where the observer point is located. The module proceeds to calculate and delineate the area that can be “seen” from the observer point. The output ‘Visibility’ grid data layer will automatically pop up as a map view window when the calculation ends (you can view the progress of the calculation at the bottom right of the work space window). Even though the delineation of the seen area is complete, the module is still running. You can end the program execution by choosing 175 the module a second time or you can click on the ‘Modules’ title on the Menu Bar and click on the module name toward the bottom of the pop-up list of options. You will see a check mark to the left of the module name. When you click on the module name the check mark will disappear (it also disappears if you use the ‘Modules’ pop-up list approach). If you have the module selected for display in the ‘Object Properties’ window, you can also stop execution by clicking a second time on the ‘Execute’ button at the bottom of the ‘Settings’ tab display. The delineated viewshed is for the observer point in site number 1. I rename the ‘Visibility’ grid data layer to ‘Site01Vis’. The three distance zones, foreground, middle-ground, and background are delineated using the Grid-Tools/Grid Buffer. Before the three zones can be delineated for site number 1, the point shapes data layer with the point objects representing the observers must be converted from its vector format to raster, i.e., to a grid data layer. The GridsGridding/Shapes to Grid module is used for that conversion. The new grid data layer is named ‘GVOpSpviewers’. Each observer point also must be in a grid data layer as the only “feature” on the layer, i.e., all other cell data values should be “no data”. I use the Grid Calculator module to create ten new grid data layers, one for each observer point. The layer for observer number 1 is named ‘Viewer01’. Now I can develop the grid data layer for the observer point that identifies three distance zones around the point: foreground, middle-ground, and background. 176 The Grid-Tools/Grid Buffer module is executed once for each of the three zones. The distance for the foreground is from the observation point to ½ mile or 2640’. The middleground is from ½ mile to 3 miles or 2640’ to 15840’. Background is from 3 to 6 miles (15840’ to 31680’). The output grid data layer from each of these executions uses a ‘1’ data value for the “buffer”. The Grid Calculator is used to add the three layers creating a fourth one that has three data values: 3 for the foreground zone, 2 for the middle-ground zone, and 1 for the background zone. The settings for the Grid-Tools/Grid Buffer module for calculating the foreground zone are displayed in Figure 6-30. Figure 6-30. An example of the settings for the Grid-Tools/Grid Buffer module. The ‘Viewer01’ grid data layer containing a single feature (one grid cell with a data value of 1) is chosen for the ‘>>Features Grid’ parameter. The ‘Buffer Distance’ option has two choices, “Fixed” or “Cell Value”. The “Fixed” option is chosen. The width of the buffer is a specified distance. The distance is entered in the ‘Distance’ parameter value field: 2640. This is feet and represents the ½ mile width of the foreground. As noted earlier, this module is executed three times for each observer point, once for each of the distance zones. Two changes are made for the second execution of the module. The “[create]” option for the ‘Grids’ parameter are temporary, interim layers that were generated for the three distance zones. The output grid data layer from execution of the Grid Calculator was renamed ‘Vzones01’. The ‘Vzones01’ grid data layer is displayed in Figure 6-32 with the ‘GVOpSpPoly’ and ‘MasonGV’ polygon shapes data layers overlain. Figure 6-32. The ‘Vzones01’ grid data layer displaying the three distance zones for observer point 1. The next processing step is to develop a grid data layer for observation point number 1 that shows the foreground, middle-ground, and background areas within its seen area (i.e., the viewshed). The seen area for observation point number 1 is the grid data layer ‘Site01Vis’ (Figure 6-33). 178 Figure 6-33. The viewshed grid data layer for observation point number 1. The areas that can be seen are in blue. The ‘GVOpSpObserve’ and ‘MasonGV’ shapes data layers are overlain. This seen area is for observation point 1. The Grid Calculator will be used to produce a grid data layer for observation point number 1 showing foreground, middle-ground, and background areas of the seen area or viewshed. Figure 6-34 displays the settings used to produce this new grid data layer. Figure 6-34. The Grid Calculator settings used to produce a seen area grid data layer showing foreground, middle-ground, and background areas. 179 The output grid data layer produced by these settings is displayed in Figure 6-35. Figure 6-35. The output grid data layer showing foreground, middle-ground, and background areas for the viewshed for observation point number 1. The output grid data layer will be re-named ‘SeenZones01’. The above set of procedures was applied to each observation point. The last set of grid data layers, ‘SeenZones*’, are used with the Geostatistics-Grids/Zonal Grid Statistics module to develop the data for the table listing acres for land uses (industry, commercial, mineral extraction, and utilities) by the three distance zones (foreground, middle-ground, and background). There will be one table for each of the three distance zones. The tables are displayed in Figure 6-36. Figure 6-36. The visual impact tables. 180 Site number 5 is the only one of the ten sites from which an “incompatible” land use is visible in the foreground. Site number 8 can see more “incompatible” land use acres in the middle-ground than any of the other sites. There are more acres visible in the background from site number 10. Discussion of Results Criterion 1: Identifying the boundary for Grapeview. The boundaries for the Grapeview Fire District and Grapeview School District were merged together to define a physical boundary for the hamlet of Grapeview. Criterion 2: Identifying Open Space land use in Mason County and Grapeview. The land use map for Mason County includes the category of land use “Open Space”. This map is available as a grid data layer. The Grapeview hamlet polygon shapes data layer was converted to a grid data layer. The output from this process was a grid data layer for the Grapeview hamlet that identifies all the open space areas in the hamlet. Criterion 3: Open space parcels must be a minimum size of 5 contiguous acres. The Grid-Filter/Filter Clumps module was applied to the grid data layer for the Grapeview hamlet to remove from consideration groups of open space less than 5 acres in size. The output from this process was similar to Criterion 2 with the exception that only groups of open space 5 acres or larger are on the output grid data layer. There are 10 groupings of cells that meet this criterion. Criterion 4: Open Space land parcels with road access. Road access is access by paved or gravel road up to the boundary of the open space area or into the area. This could have been determined visually but tools in SAGA were used to make the determination. The area referred to as site #4 is the only one of the 10 sites without road access. Criterion 5: The park parcel cannot be adjacent to these land uses: industry, commercial, mineral extraction, and utilities. A table was developed identifying adjacent land use categories for each open space site. Site #5 is adjacent to a utility land use. This is the only site adjacent to an incompatible land use0. Criterion 6: The park parcel may be adjacent to these land uses: agriculture, commercial recreation, public, rural and residential. Using the total boundary length for each site, the percent of land use, by length, was determined. Using percent, Site #5 is adjacent to more of these land uses than any of the sites followed by #10, #8, #4, #2, #3, #7, #6, #9, and #1. Criterion 7: The park is most suitable adjacent to these land uses: open space, forest, rural, and water. 181 Using the total boundary length for each site, the percent of land use, by length, for the sites was determined. Using percent, site #4 is adjacent to more compatible land uses followed by #6, #9, #7, #1, #2, #3, #10, #5 and #8. Criterion 8: The highest amount of ground surface with slopes between 0 and 3 degrees. The Geostatistics-Grids/Zonal Grid Statistics module was used to develop a table displaying the acres of slope less than 3 degrees and the percent of area for each site with slopes less than 3 degrees. Parcels with higher percentage of low slopes are more suitable for a park. The sites rank from best to worst as follows: 1, 8, 2, 6, 9, 3, 7, 5, 10, and 4. Criterion 9: The viewing of these land use classes in the foreground and middle-ground will be minimized: industry, commercial, mineral extraction, utilities. This criterion will be used to compare parcel choices. The only site determined to have an “incompatible” land use visible in the foreground was site #5. More “incompatible” land use acres in the middle-ground are visible from site #8 and in the background from site number 10. The community park committee was able to clearly eliminate site #’s 4 and 5 from consideration. The results of applying the other criteria helped them become familiar with the available open space areas and the environmental situation of each site relative to the other sites. They realize that they need to develop criteria that can be combined together and produce an integrated set of values to base a decision. A first attempt at accomplishing this involved assigning weights based on the criteria. The attribute table linked to the ‘GVOpSpPoly’ polygon shapes data layer was opened with an older version of the Microsoft Excel program. Attributes were added for criteria 4 through 9. This could have been done using tools in SAGA supporting tables, but it was felt it might be easier using Excel. The committee decided to use a weighting scheme for each criteria based on 0 through 100, with 100 being the best. They reviewed the results of applying the various criteria and assigned weights to each site accordingly. Figure 6-37 displays this portion of the attribute table. 182 Figure 6-37. Weights assigned to the sites based on applying the criteria. The first column displayed on the left is the site number. Criteria 9 involved the three visual zones (foreground, middle-ground, and background). Rank weights were assigned for each zone and then averaged for the site. The Grid – Gridding/Shapes to Grid module was then used to create a grid data layer for each criteria using the “weights” as data values. Figure 6-38 displays the settings page used for creating the grid data layer for criteria #4 (the layer was renamed ‘GVrec4’). Figure 6-38. The Grid – Gridding/Shapes to Grid module settings page for creating the grid data layer for criteria #4. You can see that the ‘GVOpSpPoly’ polygon shapes data layer is chosen for the ‘>>Shapes’ parameter and the attribute for criteria #4 (“CRIT4”) is chosen for the 183 ‘Attribute’ parameter. Using the same module, a new grid data layer was developed for each criteria. The Grid – Calculus/Grid Calculator module “integrated” the six new grid data layers into a single data layer using a simple addition equation and dividing the sum for each grid cell by 6 (the number of grid data layers). Figure 6-39 displays the Grid Calculator settings for this step. Figure 6-39. Using the Grid Calculator to display a composite score for all criteria for each site. The six criteria grid data layers are chosen for the ‘>>Grids’ parameter. The equation refers to the six layers using letters “a” through “f”. After adding the six data values together, the sum is divided by six to determine the overall average for each site. Figure 6-40 displays the output grid data layer from the Grid Calculator execution. 184 Figure 6-40. The site grid data layer displaying the average weight for each site. Using this approach the preferred site for a park is #6 followed by 2, 1, 7, 8, 9, 3, 10, 4 and 5. The final weights for each site are displayed in the legend for the map in Figure 640. Their initial application of SAGA to help them explore criteria has been successful. They are starting to understand that they need to refine their preferences within a criterion as well as between criteria. Based on the descriptions included in SAGA, other modules they intend to explore, to refine their results and improve the analysis, include the Grid – Analysis/Analytical Hierarchy Process and Ordered Weighted Averaging (OWA). Summary SAGA modules have been applied to provide outputs in support of hypothetical criteria for evaluating proposed community park sites. The same modules as well as others not used in this example could probably be applied more efficiently. As you become familiar with how to use a module you often discover more effective analysis approaches; approaches that save time or ones that produce more specific results. 185 Chapter 7 – Using SAGA With Digital Satellite Images The topic of this chapter is using SAGA functions and modules to accomplish standard tasks related to digital satellite imagery. There are core functions related to the display of grid data layers that lend themselves also to the display of digital satellite imagery. Several of the module libraries contain modules developed for use with remotely sensed data. In addition, there are modules developed as utility tools that can effectively be used with image data. This chapter will expose you to six specific areas related to the display and processing of digital imagery, and more specifically, digital satellite imagery. This does not mean the discussion does not have relevance with non-image grid data layers, it just means the focus in this chapter will be on the use of SAGA modules with digital satellite imagery. The six sections in this chapter, and the modules addressed in each, are: Contrast Enhancements Contrast enhancements are techniques used to change the way the information contained in an image displays. The goal is to alter image contrast or use color to emphasize a particular characteristic or information category of interest. Filters Filters are used to selectively emphasize or de-emphasize information in an image. Grid – Filter/User Defined Filter (3 x 3) Grid – Filter/Simple Filter Grid – Filter/Gaussian Filter Grid – Filters/Laplacian Filter Grid Calculus/Grid Calculator Grid – Filters/Majority Filter Ratio Images A ratio image is an image created by dividing the image pixel values of one satellite band by the corresponding image pixel values of another band in the same scene. In SAGA, image bands become grid data layers and pixels become grid cells. Grid – Calculus/ Grid Calculator Shapes – Tools/Create New Shapes Layer Shapes – Grid/Get Grid Data for Shapes Composite Images You can create a composite image by assigning combinations of spectral ranges or bands to the primary display colors red, green, and blue (RGB). Basically, creating a color composite is an enhancement involving multiple bands rather than a single band. Grid – Visualisation/RGB Composite Vegetation Indices 186 Vegetation indices are used with multi-spectral imagery to estimate the likelihood or probability that vegetation was actively growing at a particular location when vegetation reflectance values were collected. Grid – Calculus/Grid Calculator Image Classification Image classification is the process of converting continuous reflectance values into data categories. This process is commonly associated with multi-spectral satellite digital imagery but can also be applied to non-image grid data layers. Imagery - Classification/Cluster Analysis for Grids Imagery - Classification/Supervised Classification Shapes-Tools/Create New Shapes Layer Much of the content in this chapter is derived from the module documentation you can view in Volume 3 of this updated User Guide for SAGA 2.0. If you seek more detailed information on how to execute a particular module used in this chapter, it is likely that you can locate the module documentation in Volume 3. Imagery Background Information The sensor on board the Landsat-7 satellite is the Enhanced Thematic Mapper Plus (ETM+) that was constructed by Santa Barbara Remote Sensing under contract to NASA. This sensor was a variation of the Thematic Mapper flown on Landsats 4 and 5. A Landsat-7 scene covers 185 km on a side. The spectral bands, spectral bandwidths, and resolution for the ETM+ sensor are listed here. Band # Band 1 Band 2 Band 3 Band 4 Band 5 Band 6 Band 7 Band 8 Bandwidth (μm) 0.45 – 0.52 0.53 – 0.61 0.63 – 0.69 0.78 – 0.90 1.55 – 1.75 10.4 – 12.5 2.09 – 2.35 0.52 – 0.90 Name Blue Green Red Near IR Mid-IR Thermal IR Mid-IR Panchromatic Resolution 30 m 30 m 30 m 30 m 30 m 60 m 30 m 15 m Band 6, the thermal IR band, includes both high and low gain settings. Band 8, the panchromatic band, covers most visible to near-IR in a single band. The ground resolution for data collection (except for bands 6 and 8) is 30 meters. This means that the ETM+ sensor records 6 spectral values for each 30-meter “area” on the earth’s surface it passes over. Three of these values are reflected visible light (i.e., blue, green and red). Three of the values are reflected invisible infrared light (i.e., near IR and two in the mid IR range). A seventh value is for emitted thermal infrared. This value can be associated with ground temperature levels providing ground truth is available. The 187 ground resolution for the thermal infrared band (6) is 60 meters. The eighth band is for reflected light from visible blue through visible red through near infrared. The ground resolution for band 8 is 15 meters. The Olympic Peninsula Landsat Scene The Olympic Peninsula Landsat scene is displayed on the left in Figure 7-1. The date of this scene is September 22, 2001. A sub-scene of this scene was created that covers Mason County, Washington. The sub-scene area is displayed on the right. The polygon shapes data layer for Mason County is displayed with the red boundary on both the full scene and sub-scene. Figure 7-1. The full Olympic Peninsula Landsat scene (on left) and Mason County sub-scene on right. Descriptive information comparing the full scene and sub-scene is displayed in Figure 72. Full Scene Sub-scene # of cells 8341 (x) * 7441 (y) = 62065381 1872 (x) * 2032 (y) = 3803904 Cell size 30 meters 30 meters West-East 345615 – 595815 = 250200 460552 – 516682 = 56130 South-North 5141985 – 5365185 = 223200 5213286 – 5274216 = 60930 Mem. Size 59.190MB 3.628MB Figure 7-2. Comparing information for the full scene and sub-scene. The sub-scene is one of two sub-scenes that will be used in many of the examples in this chapter. The full scene contains too much data for effective use with my desktop system. The sub-scene is much smaller in data volume at less than 1/10th the size of the full scene. The Central Arizona Landsat Scene The second scene is one that covers the central east portion of the state of Arizona. This scene was captured on September 26, 1999. The full central Arizona Landsat scene is displayed on the left in Figure 7-3. A polygon shapes data layer for the sub-scene boundary is overlain using red for the boundary on both the full scene and sub-scene. 188 Figure 7-3. The full central Arizona Landsat scene (on left) and sub-scene on right. Descriptive information for the full scene and sub-scene is displayed in Figure 7-4. Full Scene # of cells 8081 (x) * 7061 (y) = 57059941 Cell size 30 West-East 558615 – 801015 = 242400 South-North 3725985 – 3937785 = 211800 Mem. Size 54.417MB Sub-scene 2595 (x) * 2437 (y) = 6324015 30 628995 – 706815 = 77820 3757425 – 3830505 = 73080 6.031MB Figure 7-4. Descriptive information for the full scene and sub-scene. The sub-scene for Arizona is little larger than the Mason County sub-scene for the Olympic Peninsula but still significantly smaller in memory space than the full scene. These two satellite images provide coverage of two very different geographic areas. The Olympic Peninsula of Washington State includes the Olympic Mountains, the Olympic National Park, the Olympic National Forest, the Olympic rain-forest, a portion of the Pacific Ocean coastline, the Straits of Juan de Fuca and portions of Puget Sound. It has a moderate mid-latitude west coast climate, etc., having some precipitation extremes. Vegetation ranges from lowland to higher elevation conifer forests to tundra in the higher elevations of the Olympics. The slender fingers of the Puget Sound achieve their southern-most reach toward the capital city of the state of Washington (Olympia). East of the peninsula, on the eastern edge of the image, is the large metropolitan area consisting of the cities of Olympia-Tacoma-Seattle-Everett (from south to north). The majority of the population of Washington State resides in this urban environment. In contrast to the Olympic Peninsula is the central east Arizona sub-scene. The area it covers is semi-arid and generally high-elevation that includes portions of the Little Colorado Plateau and the White Mountains. Much of the landscape is above 7000’ of elevation and characterized as being “above” the rim, i.e., the Mogollon Rim that crosses Arizona in generally an east-west direction. One of the largest expanses of Ponderosa Pine forest is found in this area. The lower elevations are used for open range for cattle 189 grazing, although vegetation is quite sparse. There are some open range areas where it takes over 20 acres to support one animal. Higher elevations support mixed conifer forests; lower elevations support pinyon and juniper. Sparse population is the rule. The twin-cities of Eagar and Springerville have a total population of around 6000. The area is roughly equidistant between Albuquerque, New Mexico and Flagstaff and Phoenix, Arizona; literally three or four hours driving time to anywhere. The Apache-Sitgreaves National Forests dominate the vegetated areas with over 2 million acres, the majority covered by this sub-scene. Contrast Enhancements Contrast enhancement is often necessary and essential to visual analysis involving digital satellite imagery. The topics in this chapter all fit within the general definition of contrast enhancement to support interpretation of digital imagery. They all assist in visually interpreting digital imagery. The simplest of the techniques involve the stretching of the data value range or “slicing” of the data value range so that certain features in the image can be more easily distinguished from their surroundings. There is basic descriptive statistical information available for all grid data layers in SAGA. The first place to view this basic information is using the ‘Description’ tab at the bottom of the ‘Object Properties’ window for an active grid data layer. For example, Figure 7-5 displays the descriptive information available for the grid data layer containing spectral reflectance data for band 3 of a Landsat TM sub-scene. 190 Figure 7-5. Basic descriptive information for a grid data layer containing spectral data for band 3 of a Landsat TM sub-scene. Three pieces of information, in Figure 7-5, important to contrast enhancement, are the ‘Value Range’, ‘Arithmetic Mean’, and ‘Standard Deviation’. Referring to Figure 7-5, the ‘Value Range’ defines the minimum and maximum data values for data in the grid data layer. In this case, the grid data layer contains values for spectral reflectance in band 3. You can see that the range is 238 when the minimum data value is 16 and the maximum value is 254. The possible range for the data is 256 since the lowest possible reflectance value is 0 and the highest is 255. The ‘Arithmetic Mean’ and ‘Standard Deviation’ statistics explain the distribution of the data values around the arithmetic mean for the data values. These two statistics can be of more value than the ‘Value Range’. In the case of reflectance values, the distribution is assumed to be a normal distribution. Certain assumptions are valid with a normal distribution of data. Approximately 68% of the data values, in a normal distribution, will be within plus or minus 1 standard deviation of the arithmetic mean. Using this example means that 68% of the data values will be within the data range of 20.66 to 36.83. Eighty-seven percent of the data values will be within 1.5 standard deviations of the mean. Further, approximately 95.45% of the data values will be within plus or minus 2 standard deviations of the arithmetic mean, i.e., within the data range of 12.57 to 44.92. Viewing the descriptive statistics in Figure 7-5, we can see that the lower value for the plus or minus 2 standard 191 deviation range is below the minimum data value for the data. This means the data values for the grid data layer are skewed toward the lower end of the distribution. Figure 7-6. The histogram for spectral reflectance values for band 3 (grid data layer ‘MCL720020922_3’). The histogram for band 3 reflectance values displayed in Figure 7-6 supports the conclusion that the data values are skewed to the lower end of the value range. The form of the histogram also approaches the form for a normal distribution, i.e., a bell-shaped curve. Figure 7-7 displays a tabular version of the histogram. 192 Figure 7-7. A portion of the tabular version of the histogram in Figure 7-6. Using the full table of the tabular data of the histogram, we can do a quick check to see what the approximate percentage is for data values falling within 2 standard deviations of the mean. This will be the total count from row 1 in the table to row 31 divided by the total number of cells in the table. The percentage comes out to 94.9%. That is fairly close to the prediction of 95.45%. There are statistical tests that can be applied to test this conclusion. This helps us to understand the overall data distribution. It means that a very small portion of the data values range between 44 and 254 data values. 193 SAGA has some basic functions in its’ core software that can be used to change the visual appearance of a grid data layer. These core functions do not change data values; rather, they change the appearance of the data through the use of color methods using knowledge of the arithmetic mean and standard deviation of the data distribution. Several of these functions are accessible in the ‘Classification’ sub-menu when you rightclick with the mouse on the grid data layer name appearing in the ‘Data’ tab area of the Workspace. Figure 7-8 displays the ‘Classification’ sub-menu. Figure 7-8. Functions available in the ‘Classification’ sub-menu. Create Normalised Classification This option applies a normalization routine to the data range and produces a lookup table using a normalized distribution. The normalized classification applies a contrast stretch emphasizing the portion of the spectrum having the majority of the data. The “Lookup Table” option in the Display: Color Classification: Type setting must be chosen and the ‘Apply’ button clicked on (at the bottom of the ‘Object Properties’ window) in order for the effect of the normalized lookup table to be viewed. The “normalized” table can be reviewed by clicking on the Display: Color Classification: Lookup Table: Table: Table (columns…) setting. When the ellipsis appears, click on it and the table will display in the work area. 194 Figure 7-9. Comparing corresponding zoomed in areas on the input and output following applying the ‘Create Normalised Classification’ option to a grid data layer. The data values do not change. The colors assigned to various ranges of data values change. Create Lookup Table The second option in the ‘Classification’ sub-menu, ‘Create Lookup Table’, builds a color look-up table for displaying the data values for the grid data layer. When you choose this option, the ‘Create Lookup Table’ window in Figure 7-10 appears. Figure 7-10. The ‘Create Lookup Table’ window. The output using the setting for the ‘Create Lookup Table’ in Figure 7-10 is displayed in Figure 7-11. Again, the data values do not change; the colors assigned for displaying the data value change because of the new lookup table color assignments. 195 Figure 7-11. Comparing zoomed in portions of input and output grid data layer after using the ‘Create Lookup Table’ option. The entry to the right of the ‘Colors’ parameter (in Figure 7-10) displays the thumbnail representation of the selected color palette that will be used for the color assignments. The number displayed, in this example it is 100, is the number of classes the data values will be divided. A different number can be entered. First click in the value field to the right of the ‘Colors’ parameter. Next, click on the ellipsis on the right side of the value field. The ‘[CAP] Colors’ window in Figure 7-12 will display. Figure 7-12. The ‘[CAP] Colors’ window. 196 When you click on the ‘Count’ button (see Figure 7-12) an ‘Input’ window will appear allowing you to change the ‘Count’ numeric value to a new value. This is the parameter that controls the number of classes that the lookup table will contain. Using an entry of 100 means that the range of data values will be divided into 100 classes of equal intervals. Thus, if the data range is from a minimum value of 16 and a maximum value of 254 the range is 238. Dividing 238 by 100 means that there will be 100 classes and each class interval will be 2.38 data values wide. Figure 7-13 displays the first few records of the lookup table created using the setting in Figure 7-10. Figure 7-13. A portion of the lookup table created using an entry for the ‘Count’ parameter of 100. The lookup table is accessible in the ‘Settings’ tab section of the ‘Object Properties’ window for the grid data layer. The table displays when you click with the mouse pointer on the ellipsis that appears to the right of the Lookup Table: Table: Table (Columns …) parameter in the settings. Looking at the table in Figure 7-13, you can see that the class interval for each of the classes is 2.38 data units. This approach is a very versatile function for contrast or density slicing. Let’s look at another example. In this example, I will use a grayscale tone palette and change the entry for the ‘Count’ parameter to 20. The new ‘Count’ entry means that each class interval will be 23.8. Figure 7-14 displays corresponding zoomed in areas of the before and after views of the grid data layer. Figure 7-15 displays a portion of the lookup table created using the grayscale palette and a ‘Count’ entry of 20. 197 Figure 7-14. Comparing zoomed in portions of the before and after grid data layers using a grayscale palette option with 20 classes. Looking closely, you can see more detail in the before display. The before display uses 238 of the available 256 grayscale levels to display data values. The after display uses 20 grayscale levels (or classes) to display the data values. You can see the subtle differences in the two displays. The one on the right does not smoothly display gradations between values. This is because the display classes have been generalized from 238 classes to 20. Figure 7-15. The lookup table created using an entry for the ‘Count’ parameter of 20. 198 This is a flexible tool for automatically assigning a palette of colors (or as in the example above, grayscale tones) to a range of data in a grid data layer. Once the lookup table has been created, you can make further changes. The lookup table information can be modified when you display the table. Figure 7-16. Adding color to the grayscale lookup table. Here is an example where I created a grayscale lookup table with 40 display classes. Figure 7-16 shows the upper portion of the lookup table where I added color to the first four classes (these four classes account for approximately 95% of the grid data values). Figure 7-17. Comparing the original grid data layer with the one using the modified lookup table. 199 It is easy to see how this function can be used to emphasize specific data ranges in a grid data layer using the lookup table. The color assigned to a class can easily be changed. You can also change the class intervals by editing the values that were automatically selected for class “Minimum” and “Maximum”. Set Range to Minimum/Maximum Set Range to Standard Deviation (1.5) Set Range to Standard Deviation (2.0) These three options are used to set the minimum and maximum data values that display for a grid data layer. The first option uses the actual minimum and maximum of the data range to establish the display range for the layer. When you choose this option, all data values are displayed. The options using 1.5 and 2.0 standard deviations will set the minimum and maximum data values for display based on the distribution of the data around the arithmetic mean. Using ‘Set Range to Standard Deviation (1.5)’ will identify a minimum data value that is the arithmetic mean minus 1.5 times the standard deviation. The maximum data value to display will be the arithmetic mean plus 1.5 times the standard deviation. For example, using the earlier example (the grid data layer for the spectral reflectance values for band 3 of a Landsat sub-scene), the arithmetic mean for the data is 28.74 and the standard deviation is 8.09. One and a half times the standard deviation (8.09) is 12.135. This means that the display range, when this option is selected, will include data values between 16.61 and 40.88. Data values below this range will be displayed using the palette color on the far left of the ramp and values above this range will use the palette color on the right of the ramp. The same approach is used with the ‘Set Range to Standard Deviation (2.0)’ option. Any one of the three options can be set as the default ‘Display Range’ parameter for viewing grid data layer values. This parameter is set when you click on the ‘Data’ title in the ‘Data’ tab area of the Workspace. The settings for ‘Data’ will display in the ‘Object Properties’ window in the ‘Settings’ tab area. Toward the bottom of the settings is an area titled “Grid Display Defaults” with a parameter named ‘Display Range’. When you click with the mouse pointer positioned in the value field, the list of three options displays. You choose the option you prefer by moving your mouse pointer over the option and click the mouse button. This option will remain in effect until you change it again and is retained from one work session to the next. Additional functions for changing the appearance of data values using color methods are located in the grid data layers ‘Object Properties’ window. The Display: Color Classification: Type section options are accessed using the ‘Settings’ tab portion of the ‘Object Properties’ window. Most of the options for ‘Type’ can be used for changing the visual appearance of a grid data layer. Two of these options are particularly relevant for use with digital satellite grid data layers. These are “Shade” and “RGB Overlay”. 200 Shade The “Shade” ‘Type’ has eight color choices that can be used with grid data layers. While all of the options can be used for data display, the first two are grayscale options that can be useful for displaying spectral reflectance using 256 levels of gray. These choices are “bright-dark” and “dark-bright”. The “bright-dark” choice assigns white to the “0” data value, “255” to the black value and distributes values in-between in gray scale levels from white to black. The “darkbright” choice inverts the display. That is, the “0” value is assigned to black, the “255” value to white and values in-between are assigned to gray levels from black to white. These two choices are the equivalent of assigning a grayscale palette for the Graduated Colors: Colors setting. Figure 7-18 displays band 4 of a Landsat TM sub-scene using these two “Shade” ‘Type’ choices. Figure 7-18. The “bright-dark” shade choice on the left and “dark-bright” on the right using band 4 of a Landsat TM sub-scene. In Figure 7-18, the map view on the left uses the “bright-dark” shade ‘Type’ choice and the one on the right the “dark-bright” choice. RGB Overlay The other option, in the Display: Color Classification: Type setting, of interest is “RGB Overlay”. This is a relatively quick and easy way to display a color composite without creating a new grid data layer for a color composite. When you choose the “RGB Overlay” option, information identifying grid data layers to assign to the primary colors (red, green, and blue) must be entered in the ‘RGB Overlay’ section of the settings. Figure 7-19 displays the ‘RGB Overlay’ section of the grid data layer settings. The option chosen in the value field to the right of the ‘Coloring’ parameter determines assignment options for “this”, ‘>Overlay 1’ and ‘>Overlay 2’. The ‘Coloring’ parameter has six choices. The one displayed in Figure 7-19 means that the ‘R’ or red portion of ‘RGB’ is assigned in the “composite” display to the current grid 201 data layer. The ‘Overlay 1’ and ‘Overlay 2’ settings are where you choose the grid data layers that will be assigned to the ‘G’ or green color and the ‘B’ or blue color for the RGB overlay or composite. The other five choices for ‘Coloring’ are variations of the combinations and dependent on the “this” or active layer. In this example, I am using the default choice: “red=this, green=1, blue=2”. The color red is assigned to the active grid data layer. This is the data layer with its settings displayed in the ‘Object Properties’ window. Figure 7-19. Example settings for “RGB Overlay”. Figure 7-20 (on the left) displays the composite display using the settings in the ‘RGB Overlay’ section in Figure 7-19. The active grid data layer happens to be band 4 (near infra-red spectral reflectance). The composite on the right was produced with the Grid – Visualisation/RGB Composite module. Figure 7-20. Example composite display using the “RGB Overlay” choice. Both of the composites in Figure 7-20 use band 4 for the red color, band 3 for the green color, and band 2 for the blue color. This combination of band assignments is referred to as a false color infra-red (FCIR) composite. It is one of the more popular satellite image composites. The two composites were generated using the same spectral bands for input assigned to the same RGB colors. You can see that there are some visual similarities. For example, vegetation in both composites has a red tint. The composite on the right is an output grid data layer from the Grid – Visualisation/RGB Composite module. The module has more functionality for applying a contrast stretch to each input band during the composite 202 creation process. The ‘RGB Overlay’ setting is quick and easy but you have very little control for improving the visual display other than for assigning the grid data layers to the three colors. The ‘RGB Overlay’ image looks too bluish. None of the display options explained in this section change data values (except the RGB Composite module referred to above). The options are used to change the visual appearance of a grid data layer with the intent of making it easier to visually interpret or to isolate a feature or features of interest. Filters A filter is applied to a grid data layer (for example one that represents an image) to create a new image by applying a mathematical function on the original grid cell data values and their neighbors to calculate new values. The objective of the function is determined by the values stored in a template or kernel. Kernels will normally be a 3 x 3, 5 x 5, or 7 x 7 or larger (but generally always symmetrical) matrix of cells. The kernel moves from cell to cell. It is centered over each grid cell as it is processed and a new value for the central cell is computed from all cells within the kernel. Mean Filter One of the least complex filters is the mean filter. The mean filter tends to smooth data by removing noise such as speckles in digital image data by averaging the values for a neighborhood of cells. I will use the Grid – Filter/User Defined Filter (3 x 3) module to introduce you to the mean filter. Figure 7-21 displays the settings page for this module. Figure 7-21. The settings page for the User Defined Filter (3 x 3) module. The grid data layer containing the spectral values for band 1 for the Mason County Landsat sub-scene has been chosen for the input ‘>>Grid’ parameter. 203 The 3 column by 3 row table representing the ‘Filter Matrix’ parameter is where you enter values for the kernel. In this case, for a filter to calculate the mean of a central cell, the 3 x 3 kernel is loaded with 1’s. This means that each cell value in the 3 x 3 kernel receives equal weighting. The grid cell data values for the 9 cells identified by the kernel will be totaled and divided by 9, the result being the mean for the 9 cells. The data entered for the ‘Filter Matrix’ parameter is displayed in Figure 7-22. Figure 7-22. The ‘Filter Matrix’ parameter with values to calculate mean for the User Defined Filter (3 x 3) module. Figure 7-23 displays a zoomed in portion of the input grid data layer and the corresponding area on the output ‘filtered grid’. Figure 7-23. Comparing zoomed in areas of the input and output grid data layers. 204 Let’s look at the 9 grid cells on the top left of the display. The input grid data layer values are 60, 53, 51, 58, 53, 51, 53, 51, and 50. The center cell is the one in the second row with the value 53. The 9 grid cell values total to 480. When you divide the total by 9 you get a mean value of 53.33. This is the value for the output grid cell at the center of the kernel. Figure 7-24. Visual comparison of the input and output grid data layers. The input grid data layer for the TM band 1 of the sub-scene is displayed on the left in Figure 7-24 and the output on the right. It is only when you zoom in on areas that you are able to see the “smoothing” effect on the data. Detail is lost because of the averaging done by the filter. The Grid – Filter/Simple Filter module can reproduce the same result using the settings in Figure 7-25. Figure 7-25. The Grid – Filter/Simple Filter module settings for a 3 x 3 “mean” kernel. 205 Viewing the settings in Figure 7-25 you can see that the Simple Filter module provides a different set of options than the User Defined Filter (3 x 3) module. For example, there are two options for the ‘Search Mode’ parameter: “Square” and “Circle”. The ‘Filter’ parameter provides choices to use this module as a filter for smoothing or for sharpening or edge enhancement. The size of the kernel is flexible and controlled by the ‘Radius’ parameter. The default of “1” defines a 3 x 3 kernel. A radius of “2” defines a 5 x 5 kernel, and so on. A desirable feature the Simple Filter does not support is the ability to enter values (or weights) in the kernel. Gaussian Filter The Grid – Filter/Gaussian Filter module is used to generalize or smooth data contained in an input grid data layer. This filter uses an averaging approach with different weights assigned to each grid cell in the kernel. Figure 7-26 displays settings for the Grid – Filter/Gaussian Filter module. Figure 7-26. Settings used for the example of the Grid – Filter/Gaussian Filter module. The same input grid data layer used for the earlier examples involving the Simple Filter and User Defined Filter (3 x 3) modules is used in this example. The options for this execution are the defaults of “1” for the ‘Standard Deviation’ parameter, “Square” for the ‘Search Mode’ and “2” for the ‘Search Radius’. A zoomed in portion of the output grid data layer is displayed on the right in Figure 7-27 alongside the corresponding zoomed in area of the input grid data layer. 206 Figure 7-27. Visual comparison of the input and output grid data layers. The filters discussed so far are characterized as low pass filters. Their objective is to capture the general background of an image. As noted, their primary application is to generalize or smooth data. In this respect they are used with digital images as well as with digital elevation models. High pass filters are used to isolate and emphasize the detail portion of an image. In particular, these filters are used to sharpen and enhance contrast, particularly, edge enhancement. Laplacian Filter The Laplacian filter is one of the most popular filters for edge enhancement. The SAGA Grid – Filters/Laplacian Filter module supports four user selectable options. Three of the options are referred to as “Standard kernel 1”, “Standard kernel 2” and “Standard kernel 3”. All are 3 x 3 matrices, see Figure 7-28. 0 -1 0 -1 4 -1 0 -1 0 Standard kernel 1 -1 -1 -1 -1 8 -1 -1 -1 -1 Standard kernel 2 -1 -2 -1 -2 12 -2 -1 -2 -1 Standard kernel 3 Figure 7-28. The “Standard kernel 1”, “Standard kernel 2” and “Standard kernel 3” options in the Laplacian Filter module. The three standard kernel matrices use a positive weight for the central cell value and negative or zero weights for all other cells in the matrix or neighborhood of cells. The fourth option supported in the module is called “User defined”. When you choose this option, you must use the default values or enter new values for the three parameters: ‘Standard Deviation (Percent of Radius)’, ‘Radius’, and ‘Search Mode’. Note that ‘Search Mode’ has two choices, circle or square. Figure 7-29 displays settings for the Grid – Filters/Laplacian Filter module with “Standard kernel 1” chosen. A second execution was made using “Standard kernel 2”. 207 Figure 7-29. Settings for the Grid – Filters/Laplacian Filter module. Zoomed in areas of the input grid data layer (the spectral values for band 1) and corresponding areas of the outputs using “Standard kernel 1” and “Standard kernel 2” are displayed in Figure 7-30. 208 Figure 7-30. Comparing the input grid data layer and two output grid data layers. The zoomed in areas in Figure 7-30 are centered on an airfield outside of Shelton, Washington. The graphic toward the top of the figure displays the image for band 1. This was the input grid data layer for two separate executions of the Grid – Filters/Laplacian Filter module. The bottom left output used “Standard kernel 1” and the bottom right “Standard kernel 2”. The airport runways and taxiway do not show up well in the band 1 image. However, they do show up a little better on both of the module output grid data layers. Comparing the outputs it appears that the “Standard kernel 2” option displays more contrast in the edge enhancement areas. In general, linear features on both outputs can be more easily discerned compared to the input layer. Sobel Detection Filter Another approach to edge enhancement is referred to as the “Sobel Detection” method. The first step in this approach is to smooth the image using one of the low pass filters, e.g., the mean or Gaussian filter. Next, two additional passes over the smoothed image are made using two directional kernels. These kernels are displayed in Figure 7-31. 209 1 2 0 0 -1 -2 Pass 1 1 0 -1 1 0 2 0 1 0 Pass 2 -1 -2 -1 Figure 7-31. The “Sobel Detection” directional filters. A new grid data layer is created for each pass. Next, the equation Sqrt ((Pass 1 * Pass 1) + (Pass 2 * Pass 2)) is used to create a merging of the two passes into a single grid data layer. I will use the Grid Calculus/Grid Calculator module for this last step. I use the Grid – Filters/Gaussian Filter module for the first step. The settings I use are displayed in Figure 7-32. Figure 7-32. The settings used for the Gaussian Filter module. The Grid – Filters/User Defined Filter (3 x 3) module is used for the passes involving the two directional kernels. Figure 7-33 displays the ‘Filter Matrix’ entries for the first kernel execution. 210 Figure 7-33. The values used for the ‘Filter Matrix’ for the first directional kernel. I name the output from the first execution of the Grid – Filters/User Defined Filter (3 x 3) ‘Pass1’ and ‘Pass2’ for the execution that used the second directional kernel. Execution of the Grid – Calculus/Grid Calculator module is the last step in generating the output for the Sobel Edge Detection method. Figure 7-34 displays the settings I used to merge the two outputs from the two directional kernels. Figure 7-34. The Grid – Calculus/Grid Calculator module settings. The input grid data layers for the Grid Calculator module are ‘Pass1’ and ‘Pass2’. These are the renamed output data layers from the two separate executions of the Grid – Filters/User Defined Filter (3 x 3) module. Note that in the Grid Calculator module the 211 input data layers are referred to with alpha characters when used in the ‘Formula’ parameter. Thus, “a” refers to ‘Pass1’ and “b” refers to ‘Pass2’. The formula directs the module to square the values of each data layer, add them together and then take the square root for calculating the data value for the corresponding cell on the output grid data layer. This output was renamed ‘Sobelgrid’. Figure 7-35. Viewing the output from applying the Sobel Edge Technique. The upper left graphic in Figure 7-35 is the input grid data layer for spectral values for band 1. The upper right graphic is the smoothed grid data layer produced by the Guassian 212 Filter module. The two map views in the middle are ‘Pass1’ on the left and ‘Pass2’ on the right. The bottom graphic is the merging of ‘Pass1’ and ‘Pass2’ created using the Grid Calculator module. Looking closely you can discern the emphasis within the output of each of the two kernels by looking at ‘Pass1’ and ‘Pass2’. The first kernel isolated edges with horizontal orientation while the second kernel emphasized vertical oriented edges. The output from the Grid Calculator module, at the bottom of the graphic, appears to exaggerate linear features. Majority Filter The SAGA Grid – Filters/Majority Filter module uses another smoothing technique to help eliminate noise. The kernel is used to identify the predominant or most frequent occurring data value for the neighborhood of cells and assign that value to the central cell. Figure 7-36 displays the settings for this example using the Majority Filter module. Figure 7-36. Majority Filter module settings page. The grid data layer containing spectral values for the Mason County Landsat band 1 is chosen as the input for the ‘>>Grid’ parameter. The ‘Search Mode’ parameter has two choices, “Square” or “Circle”. I chose “Square”. The “1” entry for the ‘Radius’ parameter means a 3 x 3 filter matrix is used. The default “0” was left for the ‘Threshold [Percent]’ parameter. The value used for the ‘Threshold [Percent]’ parameter sets a threshold for the number of same values needed to change the central cell value. There are 9 grid cells involved. Each grid cell accounts for a little over 11% of the total number of values (1 cell out of 9, 1/9 = 213 .111 or about 11%). If you want the central cell to change when at least 3 of the 9 cells have the same data value, you would enter “33” for the ‘Threshold [Percent]’ parameter. Figure 7-37 displays corresponding zoomed in portions of input and output grid data layers used with the Majority Filter module. Figure 7-37. Comparing the input and output grid data layers related to the Majority Filter module execution. The output grid data layer on the right appears smoothed. When you compare the airfield in the two layers, the airfield on the right has rough edges compared to the airfield on the left in the input grid data layer. Figure 7-38. Comparing data values for the input and output grid data layers related to the Majority Filter module execution. The displays in Figure 7-38 illustrate how the Majority Filter operates. The predominant value for the matrix of 9 grid cells on the left is 54. This is the value chosen by the module for the corresponding central cell on the output grid data layer on the right side of the figure. 214 “Relief” Filters The versatile Grid – Filters/User Defined Filter (3 x 3) module can be used with a digital elevation model to output a relief-like grid data layer. The ‘Filter Matrix’ settings are used to specify the illumination direction. Figure 7-39 displays two kernels that will be used to create relief-like grid data layers for two different illumination directions. -3 -2 -2 1 -1 2 Kernel 1 -1 2 4 4 2 2 1 -1 -2 Kernel 2 -1 -2 -3 Figure 7-39. Two kernels for different illumination directions. Figure 7-40. Comparing the outputs using the two kernels. It is easy to discern the differing illumination perspectives by looking at the “shadows”. The output using “Kernel 1” (on the left) has illumination that appears to be from a southerly direction while the output using “Kernel 2” has illumination from a northerly direction. An easier way to generate relief like output is to use the SAGA Terrain Analysis – Lighting, Visibility/Analytical Hillshading module. The module takes into account the slope and aspect of each cell in the computation. Ratios Introduction Using digital imagery, a ratio image is created by dividing the pixel values of one spectral band by the corresponding pixel values of another spectral band from the same image. The advantage of a ratio image over a single band is that the influence of variation within a band due to topography or shadow can potentially be minimized, as the same variation is present in all spectral bands of a scene. A ratio of two bands captures the spectral relationship regardless of topography or shadow. 215 Figure 7-41 displays typical spectral range reflectance curves for water, vegetation, and dry bare soil. The figure also plots the spectral ranges for TM bands 1, 2, 3, 4, 5 and 7 using the gray areas. These curves will be referred to in the following discussion. Figure 7-41. Spectral Reflectance Curves and TM+ Band Spectral Ranges. The Grid – Calculus/ Grid Calculator module will be used to produce examples of ratio images. The settings page for the Grid Calculator is displayed in Figure 7-42. Figure 7-42. The settings page for the Grid Calculator module. 216 Using the Grid Calculator to create ratio images is very easy. First, the grid data layers containing spectral reflectance information (grid data values) for the various TM+ bands must be loaded in the work session. Ratio images normally involve two bands, one divided by the other. Assuming two grid data layers will be chosen, the equation entered into the value field for the ‘Formula’ parameter is “a/b” (as displayed in Figure 7-42 above). The “a” and “b” variables represent the two grid data layers for the TM+ band spectral data. The two layers used in the formula are chosen for the ‘>>Grids’ parameter in the settings. You must make sure they are in the correct order to match the formula. Ratio TM1/TM2 Band 1 is in the blue spectral range (0.45 to 0.52 µm) and Band 2 is in the green range (0.52 to 0.60 µm). Band 1 is designed for penetrating water. It is useful for coastal water mapping and for discriminating between soil and vegetation, forest types, and cultural feature interpretation. Band 2 is useful for measuring green reflectance of vegetation for vegetation discrimination and vegetation health. It also is useful for cultural feature interpretation. Viewing Figure 7-41, you can see the spectral ranges for bands 1 and 2 compared to the reflectance curves for water, dry bare soil, and vegetation. The curves for all three surface features are in the low end of reflectance values. The two bands appear to be highly correlated. Figure 7-43. The Grid Calculator settings for the ratio TM1/TM2 image. The entry for the value field for the ‘Grid system’ parameter is the grid system the grid data layers representing the spectral values for the TM+ bands are a part. The grid data layers for bands 1 (‘MCL720020922_1’) and 2 (‘MCL720020922_2’) have been chosen for the ‘>>Grids’ parameter. The equation used is “a/b”. 217 The ratio of TM1 and TM2 (blue and green) results with a low contrast image. Figure 744 displays examples of this ratio using the two Landsat sub-scenes for Mason County and an area in central east Arizona. Figure 7-44. Comparing two TM1/TM2 ratioed images. The grayscale images in the upper part of the image are the two TM1/TM2 ratio images. The one on the left is a zoomed in area of the Mason County sub-scene and the one on the right a zoomed in area for the semi-arid central east Arizona. The color composite images (the bottom two map view windows) are true color composites (to be discussed later in this chapter). The composites can assist with interpretation of the ratio images. The true color composite images have good contrast compared to the TM1/TM2 ratio image. Features, such as roads and vegetation, do not show up well on the ratio images. The two towns of Springerville and Eagar, in the Arizona sub-scene, are in the center of the zoomed in area. On the true color composite you can see the greenness of vegetation as well as the airfield that is on the northwest side of the towns. These areas show up poorly on the ratio image. In order to provide a simple comparison of the spectral values and the ratio image values, I identified fifteen sample points, three sample points each representing five general cover classes: roads, water, bare soil, conifer vegetation, and deciduous vegetation. The 218 development of this point shapes data layer is described at the end of this section. Figure 7-45 displays spectral data for the sample points with the spectral values for bands 1 and 2 and the ratio image TM1/TM2 created using the Grid Calculator module. MC - 1 62 65 57 51 50 50 59 59 57 50 50 49 50 53 51 MC - 2 47 53 44 28 31 30 46 45 42 34 34 34 36 39 40 TM12 1.32 1.23 1.30 1.82 1.61 1.67 1.28 1.31 1.36 1.47 1.47 1.44 1.39 1.36 1.27 AZ-1 HWY 180 HWY 60 HWY 60 Water 1 Water 2 Water 3 Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 74 64 71 49 50 49 92 111 117 46 48 50 49 48 55 AZ-2 ARTM12 65 1.14 57 1.12 59 1.20 34 1.44 38 1.32 37 1.32 104 0.88 123 0.90 121 0.97 35 1.31 36 1.33 36 1.39 46 1.07 40 1.20 50 1.10 HWY 101 SR 3 SR 3 Water-Cush Water-Hood Water-Mason Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 Figure 7-45. Spectral values for bands 1 and 2 and the ratio image TM1/TM2. The spectral reflectance values for the Mason sub-scene in band 1 do not vary a lot nor do the ones in band 2. You can see how this explains the low contrast values for the TM1/TM2 ratio image. There are some small numeric differences but not enough for a real definitive differentiation between features. The Arizona sub-scene is a little bit different. The soils in particular are highly reflective. This shows up in the ratio image TM1/TM2. Otherwise, the contrast does not appear to vary significantly. You can see significant differences between the two sub-scenes. These differences are not just due to temporal differences. The land cover is quite different between the two areas. There is a lot more bare soil land cover in the Arizona sub-scene compared to the Mason County sub-scene. There is much more biomass in the Mason County sub-scene compared to the Arizona sub-scene. Vegetation types are different between the two. Ratio TM3/TM4 Band 3 is in the red range (0.63 to 0.69 µm) and Band 4 is in the near infrared range (0.78 to 0.90 µm). The red spectral range is the chlorophyll absorption region making this band useful in vegetation discrimination. It is also useful for cultural feature interpretation. Band 4 is in the near infrared range. It is useful for interpreting vegetation types, vigor, and biomass content. It can be useful for delineating water bodies and for soil moisture discrimination. 219 Viewing Figure 7-41, you can compare the spectral ranges for bands 3 and 4 to the reflectance curves for water, dry bare soil, and vegetation. Reflectance for water is moderate in band 3 and very low to zero for band 4. Vegetation will appear in darker tones due to low reflectance values in red and high reflectance in near infrared. Figure 7-46 displays examples of this ratio using the two Landsat sub-scenes for Mason County and the area in central east Arizona. Figure 7-46. Comparing two TM3/TM4 ratio images. The contrast for the TM3/TM4 ratio images is much better compared to the TM1/TM2 image (see Figure 7-44). The topography in the sparsely vegetated areas in the Arizona sub-scene is displaying some detail in the ratio image. Roads in both images have good contrast and vegetation has more definition. Figure 7-47 displays spectral numeric data for bands 3 and 4 and the ratio image TM3/TM4 using the same fifteen sample points described earlier. 220 HWY 101 SR 3 SR 3 Water - Cush Water - Hood Water - Mason Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 MC - 3 45 54 37 19 20 20 44 41 39 24 25 24 25 27 27 MC - 4 65 76 89 13 13 14 79 80 61 74 65 68 111 134 126 TM34 0.69 0.71 0.42 1.46 1.54 1.43 0.56 0.51 0.64 0.32 0.38 0.35 0.23 0.20 0.21 AZ-3 HWY 180 HWY 60 HWY 60 Water 1 Water 2 Water 3 Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 74 59 65 26 27 30 134 154 164 29 31 29 38 32 43 AZ-4 ARTM34 74 1.00 77 0.77 72 0.90 17 1.53 16 1.69 18 1.67 126 1.06 134 1.15 143 1.15 60 0.48 60 0.52 59 0.49 113 0.34 130 0.25 114 0.38 Figure 7-47. Spectral values for bands 3 and 4 and the ratio image TM3/TM4. Unlike the previous example using bands 1 and 2, in this example you see some variation in the spectral reflectance values. Notice the difference between roads and water, for example. Also, you can see the difference between vegetation and other features. The near infrared values (the MC-4 and AZ-4 columns in the table) differentiate between conifer and deciduous vegetation. Thus, you can see that the ratio image TM3/TM4 does quite well differentiating between water, roads, bare soil, and general types of vegetation. Ratio TM5/TM2 Band 5 is in the mid-infrared range (1.55 to 1.75 µm) and band 2 is in the green range (0.52 to 0.60 µm). Band 5 reflectance values are useful for interpretation of vegetation moisture content and soil moisture. This band can also assist in differentiating between snow and clouds. Comparing the reflectance curves for water, vegetation, and dry bare soil (see Figure 741) in these two spectral ranges you can see that dry bare soil and vegetation have relatively higher reflectance values for vegetation and soil in band 5 than they do in band 2. Vegetation will display with lighter tones in this ratio. Figure 7-48 displays examples of this ratio using the two Landsat sub-scenes for Mason County and the area in central east Arizona. 221 Figure 7-48. Comparing two TM5/TM2 ratio images. Figure 7-49 displays spectral data for bands 5 and 2 and the ratio image TM5/TM2 using the same fifteen sample points described earlier. MC - 5 61 81 54 9 9 9 87 67 63 28 30 28 37 59 57 MC - 2 47 53 44 28 31 30 46 45 42 34 34 34 36 39 40 TM52 1.30 1.53 1.23 0.32 0.29 0.30 1.89 1.49 1.50 0.82 0.88 0.82 1.03 1.51 1.42 AZ-5 HWY 180 HWY 60 HWY 60 Water 1 Water 2 Water 3 Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 82 67 66 12 12 11 188 191 156 37 39 36 48 51 81 AZ-2 ARTM52 65 1.26 57 1.18 59 1.12 34 0.35 38 0.32 37 0.30 104 1.81 123 1.55 121 1.29 35 1.06 36 1.08 36 1.00 46 1.04 40 1.27 50 1.62 HWY 101 SR 3 SR 3 Water - Cush Water - Hood Water - Mason Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 Figure 7-49. Spectral values for bands 5 and 2 and the ratio image TM5/TM2. 222 This example also displays some significant variation in the spectral reflectance values. The Mason sub-scene roads have low reflectance values in both bands that produce a low ratio in the ratio TM5/TM2 image. This ratio, however, is a lot lower than ratios for other features. Thus, roads, due to the low ratio, will stand out against the higher ratio values of other features. The Arizona sub-scene road surface material is more reflective in both bands but the contrast is low. The Arizona lake water is much lower in reflectance than water in the Mason sub-scene. The low contrast stands out against the higher contrast backgrounds. Ratio TM3/TM7 Band 3 is in the red area of the spectrum (0.63 to 0.69 µm) and band 7 is another band in the mid-infrared range (2.08 to 2.35 µm). The spectral range for band 7 is useful for discrimination of mineral and rock types. This range is also sensitive to vegetation moisture content. Roads and other cultural features appear in lighter tones due to the relatively high reflectance in the red band and lower reflectance in the mid-infrared band (TM band 7). Turbid water can be interpreted in this ratio image. Figure 7-50 displays examples of this ratio using the two Landsat sub-scenes for Mason County and the area in north central Arizona. 223 Figure 7-50. Comparing two TM3/TM7 ratio images. Figure 7-51 displays spectral data for bands 3 and 7 and the ratio image TM3/TM7 using the same fifteen sample points described earlier. MC - 3 45 54 37 19 20 20 44 41 39 24 25 24 25 27 27 MC - 7 39 53 32 9 9 9 51 41 38 14 17 14 19 26 25 TM37 1.15 1.02 1.16 2.11 2.22 2.22 0.86 1.00 1.03 1.71 1.47 1.71 1.32 1.04 1.08 AZ-3 HWY 180 HWY 60 HWY 60 Water 1 Water 2 Water 3 Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 74 59 65 26 27 30 134 154 164 29 31 29 38 32 43 AZ-7 66 49 53 13 10 12 149 159 129 23 27 21 24 25 44 ARTM37 1.12 1.20 1.23 2.00 2.70 2.50 0.90 0.97 1.27 1.26 1.15 1.38 1.58 1.28 0.98 HWY 101 SR 3 SR 3 Water - Cush Water - Hood Water - Mason Soil - 1 Soil - 2 Soil - 3 Conif - 1 Conif - 2 Conif - 3 Decid - 1 Decid - 2 Decid - 3 Figure 7-51. Spectral values for bands 3 and 7 and the ratio image TM3/TM7. 224 This ratio image is similar to the first example, TM1/TM2. Due to low contrast data values, it is difficult to differentiate between features. The reflectance values for the sample points show some variation but not a lot. The Mason sub-scene road sample point and the Arizona water sample point ratio values are significantly different from the ratio values for other features. Again, there are some small numeric differences but not enough for a real definitive differentiation between most features. The Sample Point Shapes Data Layer As noted earlier, I identified fifteen sample points with three points in each of five general cover classes: roads, water, bare soil, conifer vegetation, and deciduous vegetation. I digitized these sample points on-screen using SAGA vector support tools. This process is described in several places including Volume 1, Chapter 4 in the ‘Edit’: Tools Available for Shapes Data Layers section. First I executed the Shapes – Tools/Create New Shapes Layer module to create a blank point shapes data layer. Figure 7-52 displays the settings I used for this module execution. Figure 7-52. The Shapes – Tools/Create New Shapes Layer module settings. I entered a name for the ‘Name’ parameter, ‘MasonPoints2’. There are four choices for the ‘Shape Type’ parameter: Point, Multipoint, Lines, and Polygon. I chose “Point” and clicked the ‘Okay’ button. I displayed a false-color infrared (FCIR) composite image with the new point shapes data layer as an overlay. Since the point shapes data layer is transparent, the FCIR image is visible. In the absence of ground truth, I used my interpretation of the content of the FCIR image to on-screen digitize three sample points each for cover classes roads, water, bare soil, conifer vegetation, and deciduous vegetation. 225 The attribute table linked to the point shapes data layer will have two attributes. One will be a unique number for each sample point. The second one will be a simple text name for each point. I added this attribute information as I collected the digitized sample points. Once the four ratio images were created, and with the grid data layers for the spectral bands available in the work session, I executed the Shapes – Grid/Get Grid Data for Shapes module. This module uses a point shapes data layer as input. Grid data layers are chosen as sources for attribute information to be added to the attribute table for the point objects of the input layer. Each grid data layer theme is added as an attribute column to the existing attribute table linked to the input data layer. The settings page for executing this module is displayed in Figure 7-53. Figure 7-53. The Shapes – Grid/Get Grid Data for Shapes module settings page. The input for the ‘>>Shapes’ parameter is my ‘MasonPoints2’ point shapes data layer containing the 15 point objects. The output parameter (‘>Red’, ‘>>Green’, and ‘>>Blue’. The output parameter is the one named ‘Grids’ label is important. The first input in the list will be referred to as variable “a” in the equation entered in the ‘Options’ section. The second input will be referred to as variable “b”. Although I won’t be including more than two variables in the formulae in these examples, additional inputs would be referred to as variables “c”, “d”, etc. The equation entered for the ‘Formula’ parameter is (a-b)/(a+b). Interpreting the variables, the equation reads: (Near infrared – red) divided by (Near infrared + red) The output produced by the Grid Calculator module for the NDVI is displayed in Figure 7-60. The graphic on the left is a zoomed in area of the Mason County sub-scene and the one on the right is a zoomed in area in the Arizona sub-scene. Figure 7-60. A zoomed in portion of the sample grid data layer outputs for the NDVI. The data values range from –0.404255 to 0.764706 for the Mason County sub-scene and –0.636364 to 0.74026 for the Arizona sub-scene. Figure 7-61 displays data related to the sample points for both sub-scenes. The data includes reflectance values for the red and infrared bands and the calculated NDVI index values. 234 Figure 7-61. Sample point NDVI values. Generally, the NDVI values for the sample points tend to follow the guides for interpreting the values. Numbers above zero indicate the presence of vegetation and numbers less than zero non-vegetation. The index values for vegetation are on the higher end, above zero, and most are greater than 0.4. Ratio The RATIO vegetation index separates green vegetation from soil background by dividing the reflectance values in the near infrared band by those contained in the red band (Rouse et al. 1974). The equation for the RATIO is: RATIO = NIR Red The output grid data layer shows the contrast between the red and infrared bands for vegetated grid cells. Combinations of low red band values (due to absorption by chlorophyll) and high infrared reflectance values (as a result of leaf structure) reflectance produce high index values. 235 The Grid Calculator settings page for generating the Mason County sub-scene RATIO grid data layer is displayed in Figure 7-62. Figure 7-62. The Grid Calculator setup parameters for creating a grid data layer for the Ratio vegetation index (RATIO). The equation entered in the value field to the right of the ‘Formula’ parameter is “a/b”. The layer for band 4 (‘MCL720020922_4’) is referred to in the equation as variable “a” and band 3 (‘MCL720020922_3’) is referred to as variable “‘b”. Thus, the equation divides the near infrared layer data values by the data values in the red layer. The outputs produced by the Grid Calculator module for the RATIO vegetation index are displayed in Figure 7-63. The graphic on the left is a zoomed in area of the Mason County sub-scene and the one on the right is for the Arizona sub-scene. Figure 7-63. The sample grid data layer output for the RATIO. The data values range from 0.424242 to 7.5 for the Mason County sub-scene and 0.222222 to 6.7 for the Arizona sub-scene. 236 Figure 7-64. Sample point RATIO values. RATIO index values less than 1.0 are likely to represent non-vegetation cells and a ratio value greater than 1.0 is likely to be vegetation. The road and soil sample points for the Mason County sub-scene appear to be in the vegetation range. However, if you look at the relative values between vegetation and non-vegetation, vegetation values are higher than the values in non-vegetation areas. If you use a differentiation value closer to 1.5 rather than 1.0, most of the vegetation would be greater than that value. Using the 1.0 value with the Arizona sub-scene seems to work all right except for the rock sample points. The Ratio Vegetation Index (RVI) This was another of the early vegetation indices (Richardson and Wiegand, 1977). The equation for the RVI is: RVI = Red NIR An RVI value is calculated by dividing the reflectance value in the red band (band 3) by the value in the near infrared band (band 4). The low reflectance of red light from green 237 vegetation combined with the strong reflectance of near infrared will result with high index values associated with vegetation biomass. Comparing the RVI equation to the RATIO index equation described earlier in this section, you can see that the RVI equation is the inverse of the RATIO index. The Grid Calculator settings page for creating a RVI grid data layer for the Mason County sub-scene is displayed in Figure 7-65. Figure 7-65. The Grid Calculator setup parameters for creating a grid data layer for the Ratio Vegetation Index. The equation entered in the value field to the right of the ‘Formula’ parameter is: a/b. The grid cell data values in the ‘MCL720020922_3’ layer are divided by the corresponding cell data values in the ‘MCL720020922_4’ layer. The output produced by the Grid Calculator module for the RVI is displayed in Figure 766. The graphic on the left is a zoomed in area of the Mason County sub-scene and the one on the right is for the Arizona sub-scene. 238 Figure 7-66. The sample grid data layer output for the RVI. The data value range for this vegetation index is.1333333 to 2.357143 for the Mason County sub-scene and 0.222222 to 6.7 for the Arizona sub-scene. Figure 7-67. Sample point RVI values. The suggested index value for differentiating between vegetation and non-vegetation with the RVI is 1. Values less than 1.0 are interpreted as vegetation and values greater than 1.0 as non-vegetation. Vegetation in the Mason County sub-scene does have index values 239 less than 1 (as do the vegetation sample points in the Arizona sub-scene). However, the road and soil sample points also are below 1. Based on the sample point values, using a value around .35 would seem to discriminate better between vegetation and nonvegetation. The 1.0 value works fine for the Arizona sub-scene; the one exception being two of the rock sample points. Transformed Vegetation Index (TVI) This index was developed in 1975 and is documented by Deering et al (1975). A constant of .5 is added to NDVI values and then the square root taken of the results. This was done to reduce the probability of getting negative values. The square root was intended to produce a normal distribution rather than the approximate Poisson distribution of most NDVI index output. The equation for the TVI is: ⎛ NIR − Red ⎞ TVI = sqrt ⎜ + .5 ⎟ ⎝ NIR + Red ⎠ The Grid Calculator settings page for creating a TVI grid data layer is displayed in Figure 7-68. Figure 7-68. The Grid Calculator setup parameters for creating a grid data layer for the Transformed Vegetation Index (TVI). The equation entered in the value field to the right of the ‘Formula’ parameter is: sqrt(a+.5). The value of .5 is added to each grid cell value in the input NDVI (‘MCLndvi’) grid data layer and then the square root of the result is calculated and output to the corresponding cell in the output grid data layer. 240 Output for the two sub-scenes produced by the Grid Calculator module for the TVI is displayed in Figure 7-69. The graphic on the left is a zoomed in area of the Mason County sub-scene and the one on the right is for the Arizona sub-scene. Figure 7-69. The sample grid data layer output for the TVI. The data values range from 0.309426 to 1.124591 for the Mason County sub-scene and 0 to 1.11367 for the Arizona sub-scene. 241 Figure 7-70. Sample point TVI values. The differentiation value between vegetation and non-vegetation with this index is suggested to be 0.71; less than 0.71 is non-vegetation and greater than 0.71 is vegetation. This guide does not seem to work with the Mason County sub-scene but works well with the Arizona scene except for two of the rock sample points. Corrected Transformed Vegetation Index (CTVI) The CTVI was developed to minimize the probability of negative values in the TVI. The equation for the CTVI is: CTVI = ( NDVI + 0.5) * sqrt (absolute( NDVI + 0.5) absolute( NDVI + 0.5) The Grid Calculator settings page used to generate an example CTVI grid data layer is displayed in Figure 7-71. Figure 7-71. The Grid Calculator setup parameters for creating a grid data layer for the Corrected Transformed Vegetation Index (CTVI). The equation entered in the value field to the right of the ‘Formula’ parameter is: ((a+0.5)/(abs(a+0.5)))*sqrt(abs(a+0.5)). The value of .5 is added to each grid cell value in the input NDVI (‘MCLndvi’) grid data layer. This is divided by the absolute value of the cell value plus 0.5 and then multiplied by the square root of the absolute value of the cell value plus 0.5. The result is output to the corresponding cell on the output grid data layer. Output produced by the Grid Calculator module for the CTVI is displayed in Figure 772. The graphic on the left is a zoomed in area of the Mason County sub-scene and the one on the right is for the Arizona sub-scene. 242 Figure 7-72. The sample grid data layer output for the CTVI. The data values range from 0.309426 to 1.124591 for the Mason County sub-scene and –-0.369274 to 1.11367 for the Arizona sub-scene. MC - 3 MC - 4 CTVI Water - Cush 20 13 0.5365 Water - Hood 22 13 0.4928 Water - Mason 20 13 0.5365 Road - 1 46 73 0.8526 Road - 2 40 83 0.9217 Road - 3 35 88 0.9648 Soil - 1 48 74 0.8445 Soil - 2 59 86 0.8284 Soil - 3 61 90 0.8319 Agri - 1 48 113 0.9506 Agri - 2 40 126 1.0090 Agri - 3 64 86 0.8042 Grass - 1 41 138 1.0207 Grass - 2 36 108 1.0000 Grass - 3 27 130 1.0752 Decid - 1 28 139 1.0792 Decid - 2 27 126 1.0710 Decid - 3 28 150 1.0888 Conif - 1 23 58 0.9655 Conif - 2 25 78 1.0073 Conif - 3 26 67 0.9700 Rock - 1 82 77 0.6845 Rock - 2 74 54 0.5863 Rock - 3 70 76 0.7356 AL - 3 Water - 1 Water - 2 Water - 3 Road - 1 Road - 2 Road - 3 Soil - 1 Soil - 2 Soil - 3 Agri - 1 Agri - 2 Agri - 3 Grass - 1 Grass - 2 Grass - 3 Decid - 1 Decid - 2 Decid - 3 Conif - 1 Conif - 2 Conif - 3 Rock - 1 Rock - 2 Rock - 3 26 27 84 73 79 87 119 135 182 45 39 36 46 44 36 40 38 36 30 29 29 52 37 44 AL - 4 CTVI 17 0.5392 16 0.4942 28 0.0000 67 0.6761 74 0.6836 73 0.6423 118 0.7041 120 0.6642 148 0.6301 188 1.0553 166 1.0581 142 1.0467 150 1.0152 150 1.0229 157 1.0616 130 1.0146 125 1.0167 142 1.0467 86 0.9913 80 0.9838 65 0.9397 42 0.6274 57 0.8443 68 0.8452 Figure 7-73. Sample point CTVI values. 243 The same differentiation value for the TVI, 0.71, is suggested for use with the CTVI to discriminate between vegetation and non-vegetation. This guide does not appear to work with the CTVI data for my sample points. Normalized Ratio Vegetation Index (NRVI) The NRVI is a modification of the RVI (Baret and Guyot, 1991). The result of RVI –1 is normalized over RVI + 1. The formula for the NRVI is: Re d − 1) NIR NRVI = Re d ( + 1) NIR ( The ‘Grid Calculator’ settings page used to generate the example NRVI grid data layer is displayed in Figure 7-74. Figure 7-74. The Grid Calculator setup parameters for creating a grid data layer for the Normalized Ratio Vegetation Index (NRVI). The equation entered in the value field to the right of the ‘Formula’ parameter is: (a-1)/(a+1). The value of 1 is subtracted from the RVI value for a grid cell and divided by the RVI value plus 1. The result is output to the corresponding cell on the output grid data layer. The output produced by the Grid Calculator module for the NRVI is displayed in Figure 7-75. The graphic on the left is a zoomed in area of the Mason County sub-scene and the one on the right is for the Arizona sub-scene. 244 Figure 7-75. The sample grid data layer output for the NRVI. The data values range from -0.764706 to 0.404255 for the Mason County sub-scene and – 0.74026 to 0.636364 for the Arizona sub-scene. Figure 7-76. Sample point NRVI values. The guide for differentiating between vegetation and non-vegetation is just the opposite as the one for the NDVI. Values less than zero indicate vegetation and values greater than zero non-vegetation. If you look at Figures 7-77 and 7-78 you can quickly compare the 245 index values for sample points for the NDVI and NRVI. It looks like the indices are the inverse of each other. Discussion SAGA has two modules specifically designed for supporting the production of grid data layers depicting a variety of both slope based and distance based vegetation indices. In addition, the Grid Calculator module can be used with many of the calculation formulas. Other tools can be used to collect statistical data to assist with interpretation of map data. SAGA has a lot of resources to apply to remotely sensed data. The sample points have been valuable in comparing vegetation index output for the eight broad cover classes within a specific index and between indices. The samples for roads are probably not too helpful. Roads have a lot of linear area but are quite narrow. Consequently, pixels characteristic of roads will generally reflect a combination of reflectance for surface and adjacent cover. Using sample points for roads is probably not too effective. MC - 3 MC - 4 NDVI RATIO RVI TVI CTVI NRVI 20 13 -0.2121 0.6500 1.5385 0.5365 0.5365 0.2121 22 13 -0.2571 0.5909 1.6923 0.4928 0.4928 0.2571 20 13 -0.2121 0.6500 1.5385 0.5365 0.5365 0.2121 46 73 0.2269 1.5870 0.6301 0.8526 0.8526 -0.2269 40 83 0.3496 2.0750 0.4819 0.9217 0.9217 -0.3496 35 88 0.4309 2.5143 0.3977 0.9648 0.9648 -0.4309 48 74 0.2131 1.5417 0.6486 0.8445 0.8445 -0.2131 59 86 0.1862 1.4576 0.6860 0.8284 0.8284 -0.1862 61 90 0.1921 1.4754 0.6778 0.8319 0.8319 -0.1921 48 113 0.4037 2.3542 0.4248 0.9506 0.9506 -0.4037 40 126 0.5181 3.1500 0.3175 1.0090 1.0090 -0.5181 64 86 0.1467 1.3438 0.7442 0.8042 0.8042 -0.1467 41 138 0.5419 3.3659 0.2971 1.0207 1.0207 -0.5419 36 108 0.5000 3.0000 0.3333 1.0000 1.0000 -0.5000 27 130 0.6561 4.8148 0.2077 1.0752 1.0752 -0.6561 28 139 0.6647 4.9643 0.2014 1.0792 1.0792 -0.6647 27 126 0.6471 4.6667 0.2143 1.0710 1.0710 -0.6471 28 150 0.6854 5.3571 0.1867 1.0888 1.0888 -0.6854 23 58 0.4321 2.5217 0.3966 0.9655 0.9655 -0.4321 25 78 0.5146 3.1200 0.3205 1.0073 1.0073 -0.5146 26 67 0.4409 2.5769 0.3881 0.9700 0.9700 -0.4409 82 77 -0.0314 0.9390 1.0649 0.6845 0.6845 0.0314 74 54 -0.1563 0.7297 1.3704 0.5863 0.5863 0.1563 70 76 0.0411 1.0857 0.9211 0.7356 0.7356 -0.0411 Water - Cush Water - Hood Water - Mason Road - 1 Road - 2 Road - 3 Soil - 1 Soil - 2 Soil - 3 Agri - 1 Agri - 2 Agri - 3 Grass - 1 Grass - 2 Grass - 3 Decid - 1 Decid - 2 Decid - 3 Conif - 1 Conif - 2 Conif - 3 Rock - 1 Rock - 2 Rock - 3 Figure 7-77. Band 3 and 4 spectral and vegetation index values for sample points in the Mason County sub-scene. 246 Water - 1 Water - 2 Water - 3 Road - 1 Road - 2 Road - 3 Soil - 1 Soil - 2 Soil - 3 Agri - 1 Agri - 2 Agri - 3 Grass - 1 Grass - 2 Grass - 3 Decid - 1 Decid - 2 Decid - 3 Conif - 1 Conif - 2 Conif - 3 Rock - 1 Rock - 2 Rock - 3 AL - 3 AL - 4 NDVI RATIO RVI TVI CTVI NRVI 26 17 -0.2093 0.6538 1.5294 0.5392 0.5392 0.2093 27 16 -0.2558 0.5926 1.6875 0.4942 0.4942 0.2558 84 28 -0.5000 0.3333 3.0000 0.0000 0.0000 0.5000 73 67 -0.0429 0.9178 1.0896 0.6761 0.6761 0.0429 79 74 -0.0327 0.9367 1.0676 0.6836 0.6836 0.0327 87 73 -0.0875 0.8391 1.1918 0.6423 0.6423 0.0875 119 118 -0.0042 0.9916 1.0085 0.7041 0.7041 0.0042 135 120 -0.0588 0.8889 1.1250 0.6642 0.6642 0.0588 182 148 -0.1030 0.8132 1.2297 0.6301 0.6301 0.1030 45 188 0.6137 4.1778 0.2394 1.0553 1.0553 -0.6137 39 166 0.6195 4.2564 0.2349 1.0581 1.0581 -0.6195 36 142 0.5955 3.9444 0.2535 1.0467 1.0467 -0.5955 46 150 0.5306 3.2609 0.3067 1.0152 1.0152 -0.5306 44 150 0.5464 3.4091 0.2933 1.0229 1.0229 -0.5464 36 157 0.6269 4.3611 0.2293 1.0616 1.0616 -0.6269 40 130 0.5294 3.2500 0.3077 1.0146 1.0146 -0.5294 38 125 0.5337 3.2895 0.3040 1.0167 1.0167 -0.5337 36 142 0.5955 3.9444 0.2535 1.0467 1.0467 -0.5955 30 86 0.4828 2.8667 0.3488 0.9913 0.9913 -0.4828 29 80 0.4679 2.7586 0.3625 0.9838 0.9838 -0.4679 29 65 0.3830 2.2414 0.4462 0.9397 0.9397 -0.3830 52 42 -0.1064 0.8077 1.2381 0.6274 0.6274 0.1064 37 57 0.2128 1.5405 0.6491 0.8443 0.8443 -0.2128 44 68 0.2143 1.5455 0.6471 0.8452 0.8452 -0.2143 Figure 7-78. Band 3 and 4 spectral and vegetation index values for sample points in the Arizona subscene. The soil and rock sample points can also experience a lot of variation within a sub-scene and between sub-scenes. Most vegetation will tend to be shades of green. Soil and rock reflectance will vary depending on parent material and geology. The reflectance can vary from bright to dark and include various colors. I think using sample points along with other SAGA modules greatly assists the use of remotely sensed data. My next set of sample points would focus on vegetation cover (after all, vegetation indices are for predicting the existence of vegetation) and include more examples of deciduous, coniferous, riparian vegetation types. In addition, there may be vegetation types unique to the geographic area covered by the satellite scene. Image Classification Image classification is the process of converting continuous reflectance values into categories representing land cover or surface condition. This process is commonly associated with satellite digital imagery. Please note that, although the focus of this section is on remotely sensed data, and in particular on Landsat satellite digital imagery, the techniques described here are also used with non-image data. 247 Landsat Thematic Mapper Satellite digital imagery measures reflectance from the earth’s surface using eight spectral bandwidths. These bandwidth ranges (referred to as band 1, band 2, band 3, etc.) include several ranges representing visible light and several ranges in the near and middle infrared areas. In addition, one spectral range of emitted or thermal radiation is measured. More information on Landsat imagery is provided in the “Imagery Background Information” section at the beginning of this chapter. Image classification is the process of analyzing multiple band reflectance values for a specific small area of the earth’s surface in order to identify the land cover or a condition that it represents. There are two general image classification approaches used to classify satellite digital imagery: supervised and unsupervised. Supervised classification starts with a set of training sites. The training sites represent “ground truth”; these are areas that have a known land cover or condition. The spectral return values for these cells, i.e., the digital values representing the reflectance values for each spectral band of the sensor, provide the basic data for establishing each land cover class mean and variance. The procedure evaluates the rest of the image cells to determine which of the known land cover classes each cell fits best. The unsupervised classification procedure uses a cluster analysis algorithm to identify clusters of values that seem distinct based on the spectral reflectance values for each spectral band of the sensor. Mathematically, the clustering approach is performed using an iterative series of passes to minimize the within cluster homogeneity and maximize the between cluster heterogeneity. Unsupervised Classification The Imagery - Classification/Cluster Analysis for Grids module in SAGA supports the unsupervised classification process. The unsupervised classification relies on the software to statistically analyze the input multiple bands (or multiple grid data layers in the case of using non-image files) without benefit of what the different spectral or numeric combinations represent. Using the SAGA Cluster Analysis for Grids module the user specifies the number of classes (the default is 10) to analyze. The classification proceeds using statistical rules that are applied to the various multi-band pixels or multi-grid data layer data values. When module execution finishes, the resulting classes must be interpreted and related to ground features or grid data layer themes to assign class names. 248 Figure 7-79. Settings page for the module Cluster Analysis for Grids. The choices for the ‘>>Grids’ parameter will be the set of grid data layers representing spectral reflectance in various bandwidth areas. These grid data layers must be a part of a grid system loaded in the current work session. The ‘Grid system’ parameter is used to choose the grid system. You can see in Figure 7-79 that the grid system loaded in the Workspace to the left of the ‘Cluster Analysis for Grids’ window has been chosen. Notice the list of grid data layers below the grid system name. These grid data layers represent spectral reflectance for bands 1 through 7 of a Landsat scene, in this case, actually a subscene. This is the Mason County sub-scene used in examples throughout this chapter. Once the grid system has been chosen, the multiple grid data layers representing spectral reflectance can be chosen for the ‘>>Grids’ parameter. When you click with the mouse pointer in the value field to the right of the ‘>>Grids’ label, a small ellipsis appears on the right side. Clicking with the mouse pointer on the ellipsis displays the ‘Grids’ window, i.e., a list of the grid data layers available for the chosen project in the current work session (See Figure 7-80). Figure 7-80. The ‘Grids’ window displaying a list of grid data layers available for the active project selected in the ‘Grid system’ parameter. 249 When you move the mouse pointer over a grid data layer name in the list and click the mouse button, the name will become selected and appear highlighted. You can continue choosing layers in the same manner. Clicking on another data layer name will select it also. Both layers will be highlighted. If you click with the mouse on an already highlighted grid data layer name, it will de-select. Once one or more grid data layers are highlighted, you move them to the right side of the window by clicking on the button that displays the ‘>>’ symbols in the center of the window. You can move layers from the right side back to the left side by selecting and using the ‘Grids’ parameter. The chosen method for the ‘Method’ parameter uses the training site spectral values along with techniques or decision rules related to the methods 257 approach to identify grid cells that fall within the land cover classes supported by the training sites. The ‘>>Training Areas’ parameter value field is for identifying the polygon shapes data layer containing the “training sites” objects. Prior to execution of the Supervised Classification module, you must make sure that the polygon shapes data layer has either been loaded as part of a project in the current work session or has been specifically loaded using the Shapes/Load Shapes command in the Menu Bar ‘File’ drop down menu. When you click the mouse in the value field to the right of the ‘Training Areas’ parameter input is a polygon shapes data layer. The ‘Class Identifier’ parameter is where you choose the field in the attribute table that provides the numeric cover type identifiers. When you click in the value field to the right of the ‘Class Identifier’ label, a list of the attributes available for the chosen polygon shapes data layer in the ‘>>Training Areas’ parameter displays. The attribute you choose should be a numeric attribute containing data values for the land cover training sites. If you select a text or string attribute, the results will be unpredictable. There are three outputs created by this module: a mandatory ‘