**2. Materials and Methods**

area for each cell exceeds a threshold (i.e., maximum cross grading area), the FD8 procedure switches to D8 flow, allowing the modeling of both divergent and convergent flow [12]. Lea [14] developed an aspect driven kinematic routing algorithm that was no longer restricted to the eight nearest neighbors (i.e., cardinal and diagonal directions). This algorithm inspired both the D∞ and DEMON methods. The DEMON stream tube method developed by Costa-Cabral and Burges [15] is computationally intensive and involves routing flow downstream along tubes that expand and contract in a way where the tubes do not necessarily coincide with the cell boundaries [12]. The D∞ algorithm first calculates flow direction from the infin‐ ite set of possible flow directions around each cell, not just in the 8 cardinal directions.

The terrain attributes in references [5-7] were calculated with TAPES-G utilizing the FD8 method. TAPES-G still exists but it is no longer being supported and updated by the devel‐ opers. The Windows installation only works on legacy versions of ArcGIS. The preferred ap‐ proach for terrain analysis today has become TauDEM with the D∞ flow direction

Agriculture contributes 80% of the excessive phosphorus (P) flowing to the Gulf of Mexico from the Mississippi River Basin. Among the states contributing to this environmental im‐ pact, Kentucky is the sixth largest contributor of nitrogen and phosphorus to the Mississippi River Basin [16]. Sediment is the leading cause of impairment of Kentucky's rivers and streams, impacting over 4,329 linear km with agriculture being the primary contributing source [17]. The problem of sedimentation, often underestimated, is enormous. It not only adversely impacts aquatic life [18] and impaires ecological and economic functioning of drainage systems, wetlands, streams, and rivers but also increases the risk of flooding and the need for water treatment. Properly constructed grassed waterways (NRCS Code 412), with appropriately sized vegetative filters on either side reduce considerable runoff, nu‐ trient, and sediment delivery to surface waters. Specifically, they reduce ephemeral gully erosion, which is a substantial but often overlooked problem that accounts for about 40% of all erosion in agricultural fields [19]. In one study, the installation of grassed waterways led to reductions of 39 and 82% in runoff and sediment, respectively [20]. In another study, grassed waterways reduced dissolved reactive phosphorus losses by factors ranging be‐ tween 4 to 7 and particulate phosphorus by factors ranging between 4 and 10 [21]. Because grassed waterways are so effective, the USDA Conservation Reserve Program (CRP) and Environmental Quality Incentives Program (EQIP) provide funding for this and other con‐

The objective of this study was to compare how predictions of concentrated flow erosion performed with LiDAR, survey grade RTK GPS, and 10-m USGS DEMs, and using the D8,

**1.3. Concentrated Flow Erosion, Grassed Waterways and The Environment**

algorithm.

48 Research on Soil Erosion Soil Erosion

servation practices.

**1.4. Research Objective**

FD8, DEMON, and D∞ flow direction models.

This study was conducted in five fields in Shelby County located in the Outer Bluegrass physiographic region of Kentucky. Field A (38°17´ N, 85°9´ W), B (38°18´ N, 85°11´ W), C (38°20´ N, 85°12´ W), D (38°20´ N, 85°11´ W), and Field E (38°20´ N, 85°14´ W), were 23, 36, 11, 57, and 33 ha in size, respectively. The fields had been in a no-till, corn (Zea mays L.), wheat (Triticum aestivum L.), and double-crop soybean [Glycine max (L.) Merr.] or cornwheat rotation for more than 20 yr. Soils in this region developed primarily from limestone residuum overlain with pedisediment from limestone weathered materials and loess [22].

### **2.1. Field observations**

One of the co-owners of the Worth and Dee Ellis Farms was trained by the NRCS to identify eroded zones from concentrated water flow. The farmer delineated and mapped the eroded ephemeral gullies in the five study fields using a DGPS system. Subsequently, the farmer in‐ stalled grassed waterways in these areas but did not reshape them as is typically the case when installing these conservation structures. This is an important distinction because re‐ shaping these locations would have changed the terrain attributes, making the analyses pre‐ sented in this chapter difficult or impossible to interpret.

Because these field observations were conducted by the farmer, we asked the NRCS district conservationist to validate the erosion delineations in the field. The conservationist visited the fields in 2008 and examined 42% of the eroded channels delineated by the farmer in the five study fields. He determined that all the channels would have been eligible for cost sup‐ port for grassed waterways through the continuous USDA Conservation Reserve Program CRP program if they did not already have waterways installed.

### **2.2. Acquisition of Elevation Data, DEM Creation, and Terrain Analysis**

Survey grade RTK GPS equipment was used to collect elevation data from Field D in 2000 [23], Fields E and B in 2004 [24], and in Fields A and C in 2007 [5]. Dual frequency Trimble AgGPS 214 (base station) and Trimble 5800 (rover) RTK receivers were used to create sur‐ veys in Fields A, B, C, and E. The survey for Field D was created with two single-frequency Trimble 4600 receivers and elevation was determined during post processing. The surveys were created at varying intensities with elevation measurements spaced approximately ev‐ ery 3 (Field D) and 4 m (Fields A, B, C, and E) along the direction of travel during the sur‐ vey. There were 7.5 (Field D) and 12 m (Fields A, B, C, and E) between survey passes. According to the Trimble data specifications, vertical errors for all receivers were expected to be <2.2 cm because the baselines were <1200m [25,26]. The LiDAR survey was created (courtesy of Photo Science) in November of 2009 after harvest with soybean residue and stubble on the ground. The airplane height above ground was 1.219 km and speed was 185 km hr-1. The estimated horizontal and vertical accuracies were 16 and 11 cm, respectively. The average elevation point density was 2.82 points m-2. PhotoScience converted the raw Li‐ DAR data into a triangulated irregular network (TIN). Then they converted the TIN into a 1 m raster file using natural neighbor interpolation. After the 1-m LiDAR DEM was delivered

to the University of Kentucky from PhotoScience, the grid was sampled to a 4 by 4 meters. Elevation data on 9.1-m grids for these fields were obtained from Kentucky Division of Geo‐ graphic Information (KDGI) (http://technology.ky.gov/gis/). The USGS DEMs had been pre‐ viously re-interpolated by KDGI from 10-m to 9.1 meters so that they could distribute in the Kentucky State Plane coordinate system. The 9.1-m USGS raster was converted to a point file format. The raw RTK GPS and USGS points data were converted to 4 by 4-m grids with the ArcGIS (ESRI, Redlands, CA) TOPOTORASTER command (drainage enforcement op‐ tion was not used for the reasons described in reference [5]). To smooth the data, 1-m con‐ tour maps were created from the 4 by 4-m RTK, LiDAR, and USGS DEMs with the ArcGIS spatial analyst extension. Then the TOPOTORASTER command was used to convert the contours back into new 4-m DEMs.

erosion [27]. Unfortunately, this index does not estimate the true length-slope well on longer

Terrain Analysis for Locating Erosion Channels: Assessing LiDAR Data and Flow Direction Algorithm

http://dx.doi.org/10.5772/51526

51

A new variable, "ErWater," was created that was assigned a value of 0 if the grid point obser‐ vations were not from areas with evidence of concentrated flow and a value of 1 if there was evidence. The terrain attributes for each 4 x 4 grid from all five fields (n= 99,505) were export‐ ed from ArcGIS in a database format so they could be read by SAS. The sample module in SAS Enterprise Miner 4.3 were used to random subsample the dataset with an equal num‐ ber of observations from each field and areas within the fields with and without concentrat‐ ed flow erosion. So after subsampling, the final dataset included 1450 observations (145 from each of the eroded areas and non eroded areas in each of the five fields studied). SAS PROC LOGISTIC was used for logistic regression with the 1450 point subset used as the input data‐ set. The SCORE option was used to obtain full dataset predictions with all 99,505 records. The topographic wetness index, the estimated length-slope for the Universal Soil Loss Equation, and plan curvature were included as predictor variables. The dependent variable was the new‐ ly created ErWater variable. Proc FREQUENCY was used to create confusion tables for the score datasets and the results were reported in percentages. The confusion tables presented

in this chapter can be interpreted according to the guide presented in Table 1.

non-eroded areas

Areas (Type II errors)

The scored logit data (i.e., 0's and 1's) were then exported from SAS into ArcGIS where the

Logistic regression models and statistical tests are given in Table 2. The topographic wetness index, length-slope, and plan curvature values were all highly significant in the models. The discretized output of the LiDAR D8 (Figure 3) and D∞ (Figure 4) models indicated a high correspondence between the eroded waterway boundaries and the black shaded areas (i.e.,

NE (non-eroded) Correctly classified

point data were converted to raster format for display in GIS.

E (eroded) Incorrectly classified eroded

**3.1. Performance of LiDAR data with D∞ and D8 Flow Direction Models**

**Predicted NE (non-eroded) E (eroded)**

> Incorrectly classified un-eroded areas (Type 1 errors)

> > Correctly classified eroded areas

and steeper slopes [22] but still has adequate utility.

**Actual Field Status**

**Table 1.** Confusion table interpretation guide.

**3. Results and Discussion**

**2.3. Statistical analyses**

Next, TAPESG was used to remove sinks (depressions) and calculate terrain attributes for the RTK and USGS datasets using the FD8 and DEMON flow direction algorithms. TauDEM was used to remove pits and calculate terrain attributes for the RTK, LiDAR, and USGS da‐ tasets using the D8 and D∞ flow direction algorithms. Primary terrain attributes included slope (β), plan curvature, upslope contributing area, and specific catchment area. Slope was calculated using the finite distance algorithm in TAPES and the D∞ algorithm (i.e., slope could be in any direction, not just in the cardinal and diagonal directions) with TauDEM. For the FD8 method, the maximum cross-grading area [12], was hardcoded to a value of 50,000 m2 by the authors of the TAPESG for Windows software. Each cell in the DEM had a flow width dependent on the direction of flow entering that cell from the eight neighbors. For the TauDEM D8 and D∞ procedures, flow width was assumed to be 1 grid increment in all directions. The flow width for FD8 and DEMON algorithms is described in detail in ref‐ erence [12]. The specific catchment area adjusts the upslope contributing area to account for flow width and was calculated as follows:

$$\text{Specific Catchment Area} = \frac{\text{Upslope Continuing Area}}{\text{Flow Width}} \tag{1}$$

Secondary terrain attributes (i.e., those computed from two or more primary attributes) were also determined with TAPESG. This included the topographic wetness index calculat‐ ed as

$$\text{Topographic wetness index} \newline \text{In} \left( \frac{\text{Specific attachment area}}{\tan \beta} \right) \tag{2}$$

and the length-slope factor calculated as

$$\text{Length-Slope=1.4} \left( \frac{\text{Specific Catchment}}{22.13} \right) \left( \frac{\sin \beta}{0.0896} \right) 1.3 \tag{3}$$

The length-slope terrain attribute (Eq. [3]) was developed to be an estimate of the lengthslope factor from the Revised Universal Soil Loss Equation (RUSLE), and index of potential erosion [27]. Unfortunately, this index does not estimate the true length-slope well on longer and steeper slopes [22] but still has adequate utility.

#### **2.3. Statistical analyses**

to the University of Kentucky from PhotoScience, the grid was sampled to a 4 by 4 meters. Elevation data on 9.1-m grids for these fields were obtained from Kentucky Division of Geo‐ graphic Information (KDGI) (http://technology.ky.gov/gis/). The USGS DEMs had been pre‐ viously re-interpolated by KDGI from 10-m to 9.1 meters so that they could distribute in the Kentucky State Plane coordinate system. The 9.1-m USGS raster was converted to a point file format. The raw RTK GPS and USGS points data were converted to 4 by 4-m grids with the ArcGIS (ESRI, Redlands, CA) TOPOTORASTER command (drainage enforcement op‐ tion was not used for the reasons described in reference [5]). To smooth the data, 1-m con‐ tour maps were created from the 4 by 4-m RTK, LiDAR, and USGS DEMs with the ArcGIS spatial analyst extension. Then the TOPOTORASTER command was used to convert the

Next, TAPESG was used to remove sinks (depressions) and calculate terrain attributes for the RTK and USGS datasets using the FD8 and DEMON flow direction algorithms. TauDEM was used to remove pits and calculate terrain attributes for the RTK, LiDAR, and USGS da‐ tasets using the D8 and D∞ flow direction algorithms. Primary terrain attributes included slope (β), plan curvature, upslope contributing area, and specific catchment area. Slope was calculated using the finite distance algorithm in TAPES and the D∞ algorithm (i.e., slope could be in any direction, not just in the cardinal and diagonal directions) with TauDEM. For the FD8 method, the maximum cross-grading area [12], was hardcoded to a value of

by the authors of the TAPESG for Windows software. Each cell in the DEM had a

Flow Width (1)

tan <sup>β</sup> ) (2)

0.0896 )1.3. (3)

flow width dependent on the direction of flow entering that cell from the eight neighbors. For the TauDEM D8 and D∞ procedures, flow width was assumed to be 1 grid increment in all directions. The flow width for FD8 and DEMON algorithms is described in detail in ref‐ erence [12]. The specific catchment area adjusts the upslope contributing area to account for

Specific Catchment Area= Upslope Contributing Area

Topographic wetness index=ln( Specific catchment area

Length-Slope=1.4( Specific Catchment

Secondary terrain attributes (i.e., those computed from two or more primary attributes) were also determined with TAPESG. This included the topographic wetness index calculat‐

The length-slope terrain attribute (Eq. [3]) was developed to be an estimate of the lengthslope factor from the Revised Universal Soil Loss Equation (RUSLE), and index of potential

22.13 )0.4( sinβ

contours back into new 4-m DEMs.

flow width and was calculated as follows:

and the length-slope factor calculated as

50,000 m2

50 Research on Soil Erosion Soil Erosion

ed as

A new variable, "ErWater," was created that was assigned a value of 0 if the grid point obser‐ vations were not from areas with evidence of concentrated flow and a value of 1 if there was evidence. The terrain attributes for each 4 x 4 grid from all five fields (n= 99,505) were export‐ ed from ArcGIS in a database format so they could be read by SAS. The sample module in SAS Enterprise Miner 4.3 were used to random subsample the dataset with an equal num‐ ber of observations from each field and areas within the fields with and without concentrat‐ ed flow erosion. So after subsampling, the final dataset included 1450 observations (145 from each of the eroded areas and non eroded areas in each of the five fields studied). SAS PROC LOGISTIC was used for logistic regression with the 1450 point subset used as the input data‐ set. The SCORE option was used to obtain full dataset predictions with all 99,505 records. The topographic wetness index, the estimated length-slope for the Universal Soil Loss Equation, and plan curvature were included as predictor variables. The dependent variable was the new‐ ly created ErWater variable. Proc FREQUENCY was used to create confusion tables for the score datasets and the results were reported in percentages. The confusion tables presented in this chapter can be interpreted according to the guide presented in Table 1.


**Table 1.** Confusion table interpretation guide.

The scored logit data (i.e., 0's and 1's) were then exported from SAS into ArcGIS where the point data were converted to raster format for display in GIS.
