GEOG 487
Environmental Challenges in Spatial Data Science

Part II: Customize the Data and Produce the DRASTIC Groundwater Vulnerability Layer

PrintPrint

Part II: Customize the Data and Produce the DRASTIC Groundwater Vulnerability Layer

In Part II, we will create a series of grids representing the DRASTIC Ratings for each parameter (D -Depth to Water Table, R- Net Recharge, A - Aquifer Media, S - Soil Media, T - Topography, I - Impact of Vadose Zone, and C- Hydraulic Conductivity). The dataset we will use to create each grid is shown in the graphic below. In this section, we will introduce two new spatial analyst concepts: creating slope grids from elevation and reclassifying ranges of values as opposed to unique values.

Acronym Letter and Input Dataset. D & R = provided, A = geology, S = Soil, T = elevation, I & C = geology
  1. Add DRASTIC ratings to the Soil Layer

    1. Open the Soil attribute table and examine the data. Pay particular attention to the "TEXTURE" field. We will use this field to convert the vector file to a raster grid.
    2. Go to the Analysis tab, Geoprocessing group > Tools > Toolboxes > Conversion Tools > To Raster > Feature to Raster.
      Feature to raster settings: Input features: soil, Field: TEXTURE, Output raster: soilgrid.tif, output cell size: 30
      stop icon

      Make sure you have the correct answer before moving on to the next step.

      The "soilgrid.tif" attribute table should have all of the attributes shown below. If your data does not match this, go back and redo the previous step. Be sure to go to Feature to Raster tool > Environments and double-check and the output coordinates and processing extent to the same as "LakeRaystown.” Also, be sure to expand the table columns to view all COUNT totals.
       

      Data shown respectively: (OID, value, count, texture).  (0,1,2488761, Silty clay loam), (1,2,220145, Loam), (2,3,60137, Loamy sand)
    3. Compare the "soilgrid.tif" map and attribute table to that of the "soil" shapefile.
    4. Now that the soil data is in grid format, we can reclassify the grid to assign DRASTIC Ratings to each soil type.
    5. Use the values in the table below to reclassify "soilgrid.tif" based on the "TEXTURE" field. Name the output file "s.tif" (This is the letter used in the DRASTIC acronym to represent the Soil Media). Be sure to confirm the Reclassify > Environments > output coordinates, processing extent and mask are the same as "LakeRaystown.”
      Table 1: DRASTIC Ratings for Soil Textures (S)
      Texture DRASTIC Rating
      Silty Clay Loam 3
      Loam 5
      Loamy Sand 6
      Input raster: soilgrid.tif, Reclass field: TEXTURE, Reclass: same data in table above. Output raster: S.tif. change missing values to nodata
      stop icon

      Make sure you have the correct answer before moving on to the next step.

      The "s.tif" attribute table should match the example below. If your data does not match this, go back and redo the previous step.

      Data shown respectively: (OID, value, count, texture).  (0,3,2488568 , Silty clay loam), (1,5, 220140, Loam), (2,6,60131, Loamy sand)
  2. Add DRASTIC ratings to the Geology Layer

    Three of the seven DRASTIC factors (A - Aquifer media, I - Impact of the vadose zone, and C - Hydraulic Conductivity) can be defined on the basis of geology. We will use the Reclassify Tool again to assign DRASTIC ratings corresponding to these three factors for the appropriate surface geology units contained in the geology layer.

    1. Open the Geology attribute table and examine the data. Pay particular attention to the "ROCK_TYPE" field. We will use this field to convert the vector file to a raster grid.
    2. Convert the Geology shapefile to a grid. Use "ROCK_TYPE" as the "Field." Name the grid "geologygrid.tif"
      stop icon

      Make sure you have the correct answer before moving on to the next step.

      The "geologygrid.tif" attribute table should have all of the attributes shown below. If your data does not match this, go back and redo the previous step.

      screenshot geologygrid.tif data. See accessible table below
      Click here for an accessible text version of the image above
      Accessible geologygrid.tif dataset
      OID Value Count Rock_type
      0 1 1480784 Interbedded Sedimentary
      1 2 643096 Sandstone
      2 3 388372 Shale
      3 4 256791 Carbonate

    3. Tables 2, 3, and 4 show the DRASTIC ratings for Aquifer Media, Vadose Zone, and Hydraulic Conductivity, respectively. Create three new grids from the "geology grid" raster .tif using the reclassify tool and "ROCK_TYPE" field. Name the new grids "a.tif," "i.tif," and "c.tif". Remember to confirm the Reclassify > Environments >output coordinates, processing extent and mask are the same as "LakeRaystown.”
      Table 2: DRASTIC Ratings for Aquifer Media (a)
      Rock Type DRASTIC Rating
      Interbedded Sedimentary 6
      Sandstone 6
      Shale 2
      Carbonate 10
      Table 3: DRASTIC Ratings for Vadose Zone (i)
      Rock Type DRASTIC Rating
      Interbedded Sedimentary 6
      Sandstone 6
      Shale 3
      Carbonate 10
      Table 4: DRASTIC Ratings for Hydraulic Conductivity (c)
      Rock Type DRASTIC Rating
      Interbedded Sedimentary 2
      Sandstone 1
      Shale 1
      Carbonate 10
      stop icon

      Make sure you have the correct answer before moving on to the next step.

      The "a," "i," and "c" attribute tables should have all of the attributes shown below. If your data does not match this, go back and redo the previous step. Again, be sure to expand the COUNT field to see all the complete values.

      a Data shown respectively: (OID, value, count, rock_type). (0,2,388363, shale), (1,6,2123687 sedimentary/sandstone), (2,10,256789, carbonate)
      i Data shown respectively: (OID, value, count, rock_type). (0,3,388363 shale), (1,6,2123687 sedimentary/sandstone), (2,10,256789, carbonate)
      c Data shown respectively: (OID, value, count, rock_type). (0,1,1031329 sandstone/shale) (1,2,1480721 sedimentary), (2,10,256789, carbonate)
  3. Create a Slope Map from the Digital Elevation Model and add DRASTIC ratings

    When you have data that represents elevation, you can create several different types of raster layers, one is a slope grid. Slope represents steepness, incline, or grade of a line or area. A higher slope value indicates a steeper incline. With Spatial Analyst, it is easy to create a slope layer from elevation data.

    1. Go to the Analysis tab, Geoprocessing group > Tools > Toolboxes > Spatial Analyst Tools > Surface > Slope.
    2. Choose "elev" as the "Input raster", select "PERCENT RISE" for the "Output measurement." Leave the "Z factor" at 1 and name the output raster "pctslope.tif". The resulting layer depicts steep slopes with high values and gentle slopes with low values. Remember to confirm that the Slope > Environments > output coordinates, processing extent, and mask are the same as "LakeRaystown.”
       
      Geoprocessing slope settings: Input raster: elev, output raster: pctslope.tif, output measurement: Percent rise, Method: Planar, Z factor: 1
      information icon

      Degree vs. Percentage

      Be careful when choosing the slope output measurement. There are two ways to express slope values, either as a percent or as a degree. "45 degrees" slope and "45 %" slope are NOT equivalent values.

      Degree slope (θ): angle created by a right triangle with sides of length "rise" and "run"

      Percent slope: length of "rise"/length of "run" * 100

      Rise Run right Triangle. Rise going up and run going to the right. Angle marked as far right acute angle
    3. Examine the "pctslope.tif" grid. Notice how the attribute table is greyed out. Remember from Lesson 5 that raster attribute tables are not created if the values contain decimals.
    4. We want to reclassify the "pctslope.tif" grid using the DRASTIC Ratings in Table 6.
      Table 6: Ranges and Ratings for Topography
      Topography Range DRASTIC Rating
      0-2 10
      2-6 9
      6-12 5
      12-18 3
      >18 1
    5. Open the Reclassify tool and select the "pctslope.tif" grid. Notice the default number of classes and break values listed in the "Start" and "End" columns. These are not particularly useful to us, since we want to use 5 break values (2, 6, 12, 18, and the largest number in the dataset).
    6. The quickest way to change these settings is to click on the "Classify" button. Change the number of classes to "5." Manually type in the break values. "
       Screenshot reclassify window highlighting the classify option
    7. Modify the "New Values" in the reclassify window based on the values in Table 6. Name the resulting grid "t.tif". Make sure you check the box "Change missing values to NoData" and confirm the Reclassify > Environments > output coordinates, processing extent and mask are the same as "LakeRaystown.”
      stop icon

      Make sure you have the correct answer before moving on to the next step.

      The "t" attribute table should have all of the attributes shown below. If your data does not match this, go back and redo the previous step.


            "t" Data shown respectively: (OID, value, count). Data: (0,1,847058), (1,3,544004), (2,5,780674), (3,9,468931), (4,10,116560)
      info icon

      Reclassifying Ranges of Numbers vs. Unique Values

      When you need to reclassify data based on ranges of values instead of unique values. For example, notice above that the old value of "2" is specified as the upper bound in the range "0-2" and the lower bound in the range "2-6." What new value, either "10" or "9," will be assigned to old values of "2" in the output grid?

      In this case, ArcGIS will assign the old value "2" to a new value of "10," and the old value of "2.0001" to a new value "9" in the output grid. The general rule is that ArcGIS will include the break values themselves in the group that it forms the upper range boundary. Notice that you will encounter this same issue for all break values (e.g., "6", "12", and "18" in the example above).

      This is particularly important when the break values themselves are meaningful in your analysis. The most common example of this situation is when you encounter specifications of "less than x" vs. "less than or equal to x" in your requirements. If you want to reclassify values "less than 5" to a new value, you would need to specify a break value of "4.99999999," so the value of "5" is not included in your new category. The particular number of decimals you need to specify will depend on the number of decimals in your input data. For example, if your data layer has five decimal places, then you would set the reclassification thresholds as follows: a.aaaaa - b.bbbbb, b.bbbbb - c.ccccc, and so forth.

      See the ArcGIS Help for further information regarding reclassification by range.

  4. Explore the DRASTIC Rating Output Grids

    1. Add the "r" and "d" grids from your L8 folder and open the attribute tables. The numbers in the "VALUE" fields correspond with DRASTIC ratings based on each cell’s recharge rate and depth to groundwater, respectively. These are the final input datasets required to calculate DRASTIC groundwater vulnerability ratings for the watershed.
    2. Change the symbology of the d, r, a, s, t, i, and c grid layers so that low vulnerability ratings (1-3) are green, medium ratings (4-7) are yellow, and high ratings (8-10) are red. You can classify similar values together using "Unique Values" on the Symbology pane > Press the Ctrl key and click to select rows under Value heading to group > Right-click select Group values.
    3. Update the label column as shown below so the results are easier to interpret.
      Table update: Green, 2; 3, low vulnerability, yellow, 5; 7, medium vulnerability, red, 9; 10, high vulnerability
      lightbulb

      Compare the "d" grid to the "streams_buffer" shapefile. Do areas near streams have high or low vulnerability?

      Which input datasets (d, r, a, s, t, i, c) have the highest DRASTIC rating values?

      Do you see any spatial patterns in the individual drastic grids?

    4. Remove all of the datasets other than the drastic grids, the steam buffers, and the watershed boundary so your map is easier to work with. Save your map.
  5. Calculate the DRASTIC Groundwater Vulnerability Index

    Now that you have the required data layers, you can create a DRASTIC groundwater vulnerability grid based upon the DRASTIC index equation. This will involve use of the Raster Calculator to combine several grids in a weighted overlay. The graphic below shows an example of how cell values are updated during the calculation.

    Shows how values are updated.  “S” x Factor of 2 = Weighted “S” & “A” x Factor of 3 = Weighted “A” and also Weighted “S” + Weighted “A”
    information icon

    Combining raster layers is a simple, yet very important process with Spatial Analyst. You will often find that it is necessary to create a single layer that is comprised of several data sets. The idea is similar to that of performing an overlay with vector layers, in that you are making one out of many, with the major exception that the cell values change based on the expression used.

    The addition (+) and multiplication (*) signs are the most common arithmetic operators used to combine raster layers. The plus (+) sign performs an addition with each cell, so the value in a given cell of one grid will be added to the value of the same cell in the next grid, and so on. The multiplication (*) sign, as expected, performs a multiplication based on the values in each cell.

    Either of these can be used when the purpose is to simply combine grids, although you should use the same operator for all grids. However, when forming an expression that includes additional operations on individual grids, as in the case above, it is important to understand the precedence that the operators will be performed. In mathematical order of operation rules, multiplication always takes precedence over addition. Hence, in the expression above, the values in the "D" grid will be multiplied by 5 before they are added to the values in the "R" grid. If an expression should occur that is out of precedence, enclose that expression with parentheses, as you would when using a calculator.

    1. Open the Raster Calculator Analysis tab, Geoprocessing group > Tools > Toolset > Spatial Analyst Tools > Map Algebra > Raster Calculator) and enter the expression below and name the grid "drastic_index". Be careful when choosing your input files. Also, the syntax for the raster calculator must be absolutely correct, or you will get a "syntax error."
      "d" * 5 + "r" * 4 + "a.tif" * 3 + "s.tif" * 2 + "t.tif" * 1 + "i.tif" * 5 + "c.tif" * 3
      Screenshot output raster, settings described above
    2. Change the color ramp so that high values are shades of red and low values are shades of green like the example below.
      Color Ramp, green to red
       
      stop icon

      Make sure you have the correct answer before moving on to the next step.

      The" drastic_index" grid should have the following information. The statistics from the "COUNT" field are also provided. If your data does not match this, go back and redo the previous step.

      screenshot drastic index value. high: 215, low: 53

           Mean: 17,410.1329113924, Median: 2891, Std. Dev. 33,926.6924616866, Count: 158, Min: 1, Max: 177359, Sum: 2750801, Nulls: 0
       
      lightbulb

      What do the numbers in the "VALUE" field of the "drastic_index" mean in the real world? For example, do high values represent areas with high or low vulnerability to groundwater pollution?

      Which parts of the watershed are most vulnerable to groundwater pollution?

      Do any of the parameters have a greater influence on the final results?