Understanding percentage drop-raster (flow direction-tool, TWI)

[HenriRiihimäki]

I'm creating a Topographic Wetness Index using a python script (http://arcscripts.esri.com/details.asp?dbid=16750)

As I don't want this just to be a blackbox, I'm trying to understand the different rasters this script produces. One of them is outdrop-raster, which is a by-product of flow direction raster. According to ESRI:


"The output drop raster is calculated as the difference in z-value divided by the path length between the cell centers, expressed in percentages. For adjacent cells, this is analogous to the percent slope between cells. Across a flat area, the distance becomes the distance to the nearest cell of lower elevation. The result is a map of percent rise in the path of steepest descent from each cell."

..which seems reasonable. The problem is, that there are several numbers which pop out. For example 16,666.. % and 25 % are by far the most common numbers, why? I'm sure there is an explanation for this, but I just can't figure it out. It doesn't seem logical to have continuous elevation (thus slope) data and results that are not (that) continuous. Original data is low-density LiDAR DEM (2 m pix), values are 32 bit float (197.132, 198.013 m etc.).

[HenriRiihimäki]

I'll answer to myself:

It seems that all the values are rounded to the nearest 0,5. So 0,401 m / 2.8284 m = 0,14177, would instead be 0,5/3 = 16,6667. Or if distance is 2 metres, then 0,5/2 = 25 %. In the manual it says the values are rounded in FLAT areas, which seems to be a little off considering it does that where slope is 25 %. Another remark is that in flat areas (lakes) the drop percentages are indeed continuous!

[ATTACH=CONFIG]28643[/ATTACH]

Here's a sample showing what I mean. Darker gray equals flow direcetion to lower left and percentage 16,666 % drop while lighter gray here has 25 % drop and flow direction to the left.

[XanderBakker]

I'll answer to myself:

It seems that all the values are rounded to the nearest 0,5. So 0,401 m / 2.8284 m = 0,14177, would instead be 0,5/3 = 16,6667. Or if distance is 2 metres, then 0,5/2 = 25 %. In the manual it says the values are rounded in FLAT areas, which seems to be a little off considering it does that where slope is 25 %. Another remark is that in flat areas (lakes) the drop percentages are indeed continuous!

[ATTACH=CONFIG]28643[/ATTACH]

Here's a sample showing what I mean. Darker gray equals flow direcetion to lower left and percentage 16,666 % drop while lighter gray here has 25 % drop and flow direction to the left.

Hi Henri,

I had a look at the original script and your data noticed a few things:

• the script assumes a cell size of 5m


• you are working with a DEM with 2m cell size


• it uses the old arcgisscripting library

The assumption of 5m cell size is based on the following lines in the script:

• gp.times_sa(Output_Plus_raster, "25", Output_Times_raster) # calculates area of cell 5x5=25m


• Input_true_raster_or_constant_value = "5" # distance between cell centers horizontally and vertically


• false_raster_or_constant_value = "7.07" # distance between cell centers diagonally

I adopted the code to arcpy (10.2) and added some code to retrieve the cell size from the input DEM. No additional checks are added though. I used it on a DSM (including buildings and flat terrain giving high and low values) and couldn't reproduce the effects you describe. Give it a try and let me know how it went.

Example usage:
[ATTACH=CONFIG]28648[/ATTACH]

#   Created By: Prasad A Pathak
# Purpose: Calculate Topographic Wetness Index: natural logarithm of area divided by slope
# It indicates the probable water saturation level of the ground
# By default, ArcGIS calculates slope by considering just 3 X 3 neighborhood
# For this index a slope is calculated from the topmost pixel to any particular pixel is needed
# This is acheived by calculating the minimum elevation raster

# small adaptions to ArcGIS 10.2 by Xander Bakker
# - adapted to arcpy
# - cellsize is determined from input DEM
# - output TWI is name (string) not raster dataset


# Import system modules
import arcpy
import os
import math
from arcpy import env
env.overwriteOutput = True

# Get the name of folder
env.workspace = arcpy.GetParameterAsText(0)

#Input DEM
InDEM = arcpy.GetParameterAsText(1)

#Output TWI raster NAME
OutTWI = arcpy.GetParameterAsText(2) # note this was a raster, is now changed to a name

# Check out license for Spatial Analyst
arcpy.CheckOutExtension("Spatial")

# determine cellsize of input raster
desc = arcpy.Describe(InDEM)
cellsize = desc.meanCellHeight # new


# Variables used
Output_surface_raster = env.workspace + os.sep + "eldodem_fill"
Output_flow_direction_raster = env.workspace + os.sep  + "eldodem_flw"
Output_drop_raster = env.workspace + os.sep + "eldodem_drp"
Output_accumulation_raster = env.workspace + os.sep  + "eldodem_acc"
Output_Plus_raster = env.workspace + os.sep  + "eldodempls"
Output_Times_raster = env.workspace + os.sep  + "eldodemtim"
Output_Mean_elev_raster = env.workspace + os.sep  + "mean_elev"
Output_edrop = env.workspace + os.sep  + "edrop"
Out_change = env.workspace + os.sep  + "change_elev"
Out_distance = env.workspace + os.sep  + "distance"
Out_slope = env.workspace + os.sep  + "slope"
Out_preatan = env.workspace + os.sep  + "preatan"
Out_TWI_raster = env.workspace + os.sep  + OutTWI # new


# Fill DEM to make it depressionless
filledDEM = arcpy.sa.Fill(InDEM)
filledDEM.save(Output_surface_raster)
print "DEM filled"
arcpy.AddMessage("DEM filled")

# Flow Direction
flowDir = arcpy.sa.FlowDirection(filledDEM,"NORMAL",Output_drop_raster)
flowDir.save(Output_flow_direction_raster)
print "Flow direction raster created successfully"
arcpy.AddMessage("Flow direction raster created successfully")

# Flow Accumulation
flowAcc = arcpy.sa.FlowAccumulation(flowDir,"","FLOAT")
flowAcc.save(Output_accumulation_raster)
print "Flow accumulation raster created successfully"
arcpy.AddMessage("Flow accumulation raster created successfully")

# One is added to each pixel to get an count of how many pixel including the current are contributing the flow
flowAccPlus1 = flowAcc + 1
flowAccPlus1.save(Output_Plus_raster)
print "addition raster created successfully"
arcpy.AddMessage("addition raster created successfully")

# Contributing Area using the number of pixels and
flowAccPlus1Times = flowAccPlus1 * cellsize**2
flowAccPlus1Times.save(Output_Times_raster)
print "multiplication raster created successfully"
arcpy.AddMessage("multiplication raster created successfully")

# Process: Block Statistics...
meanElev = arcpy.sa.BlockStatistics(filledDEM, arcpy.sa.NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA")
meanElev.save(Output_Mean_elev_raster)
print "MIN elevation raster created successfully"
arcpy.AddMessage("MIN elevation raster created successfully")


eDrop = filledDEM - meanElev
eDrop.save(Output_edrop)
print "EDROP raster created successfully"
arcpy.AddMessage("EDROP raster created successfully")


Input_false_raster_or_constant_value = 0.005
outChange = arcpy.sa.Con(eDrop >= 0.005, eDrop, Input_false_raster_or_constant_value)
outChange.save(Out_change)
print "Change elevation raster created successfully"
arcpy.AddMessage("Change elevation raster created successfully")


inListras = arcpy.sa.InList(flowDir,[1,4,16,64])
Input_true_raster_or_constant_value = cellsize
false_raster_or_constant_value = cellsize * math.sqrt(2)
outDist = arcpy.sa.Con(arcpy.sa.IsNull(inListras),false_raster_or_constant_value,Input_true_raster_or_constant_value)
outDist.save(Out_distance)
print "Distance raster created successfully"
arcpy.AddMessage("Distance raster created successfully")


outSlope = outChange / outDist
outSlope.save(Out_slope)
print "Slope raster created successfully"
arcpy.AddMessage("Slope raster created successfully")


outPreatan = flowAccPlus1Times / outSlope
outPreatan.save(Out_preatan)
print "Pre-atan raster created successfully"
arcpy.AddMessage("Pre-atan raster created successfully")


outTWI = arcpy.sa.Ln(outPreatan)
outTWI.save(Out_TWI_raster)
print "TWI raster created successfully"
arcpy.AddMessage("TWI raster created successfully")

Kind regards,

Xander

[HenriRiihimäki]

Thanks Xander,

Good point. I also noted those and had them corrected in my script. I Should've mentioned that.
(I'm running ArcGIS 10.1.).

Here's the Pathaks code modified :

# Import system modules
import arcpy, sys, string, os, arcgisscripting
from arcpy import env  
from arcpy.sa import *

# Geoprocessor object
gp = arcgisscripting.create(10.1)

gp.overwriteoutput = 1

#############
# Input DEM #
#############
#set workspace (file directory)
env.workspace = "G:\\xxxx\\LiDAR\\TWI\\"
#set dem name (inside workspace)
InDEM = 'dem_ldm.tif'

arcpy.CheckOutExtension("Spatial")

#Output TWI
rasterOutTWI = "G:\\xxxx\\LiDAR\\TWI\\"

# Load toolbox
gp.AddToolbox("C:/ESRI/Desktop/Desktop10.1/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")

########################################################################

# Fill DEM to make it depressionless
outfill = Fill(InDEM)
print "DEM filled"


# Flow Direction
outdrop = "G:\\xxxx\\LiDAR\\TWI\\outdrop"
flow_dir = FlowDirection(outfill, "NORMAL", outdrop)
print "Flow direction raster created successfully"

# Flow Accumulation
flow_accu = FlowAccumulation(flow_dir, "",  "FLOAT")
print "Flow accumulation raster created successfully"

# One is added to each pixel to get an count of how many pixel including the current are contributing the flow
plus_raster = Plus(flow_accu, 1)
print "addition raster created successfully"

# Contributing Area using the number of pixels and 
Times_raster = Times(plus_raster, 4)
print "multiplication raster created successfully"

# Process: Block Statistics...
Output_Mean_elev_raster = BlockStatistics(outfill, "Rectangle 3 3 CELL", "MINIMUM", "DATA")
print "MIN elevation raster created successfully"

Output_Mean_elev_raster = Minus(Output_Mean_elev_raster, outdrop)
print "EDROP raster created successfully"

Input_false_raster_or_constant_value = "0,005"
Out_change = Con(outdrop, outdrop, Input_false_raster_or_constant_value, "\"VALUE\" >= 0,005")
print "Change elevation raster created successfully"

Input_true_raster_or_constant_value = "2"
false_raster_or_constant_value = "2,828427125"

Out_distance = Con(flow_dir, Input_true_raster_or_constant_value, false_raster_or_constant_value, "\"VALUE\" = 1 OR \"VALUE\" = 4 OR \"VALUE\" = 16 OR \"VALUE\" = 64")
print "Distance raster created successfully"

Out_slope = Divide(Out_change, Out_distance)
print "Slope raster created successfully"

Out_preatan = Divide(Times_raster, Out_slope) 
print "Pre-atan raster created successfully"

TWI = Ln(Out_preatan)
print "TWI raster created successfully"

TWI.save("G:\\xxxx\\LiDAR\\TWI\\twi_ykj")

So maybe I should just drop the arcscripting?

-Henri

[XanderBakker]

So maybe I should just drop the arcscripting?

Hi Henri,

I think it is a good idea to drop the arcgisscripting, since it requires loading the toolboxes and the location of those depend on the installation location, which makes it less interchangeable:

# Load toolbox
gp.AddToolbox("C:/ESRI/Desktop/Desktop10.1/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")

I think the code I attached in the previous post will work with 10.1 too (haven't tried though). Just change the lines:

# Get the name of folder
env.workspace = arcpy.GetParameterAsText(0)

#Input DEM
InDEM = arcpy.GetParameterAsText(1)

#Output TWI raster NAME
OutTWI = arcpy.GetParameterAsText(2)

into:

env.workspace = "G:\\xxxx\\LiDAR\\TWI\\"
InDEM = 'dem_ldm.tif'
OutTWI = "TWI01" # or a different name

and see if it runs...

Kind regards,

Xander

[HenriRiihimäki]

It ran trough, but drop-value problem remains.

(noticed a tiling effect, so there might be a projection problem here too -> looking into that currently)

best regards,
Henri

[XanderBakker]

It ran trough, but drop-value problem remains.

(noticed a tiling effect, so there might be a projection problem here too -> looking into that currently)

best regards,
Henri

Hi Henri,

You're right the Blockstatistics causes this effect. Change it to FocalStatistics and that problem is most likely solved (parameters remain the same):

meanElev = arcpy.sa.BlockStatistics(filledDEM, arcpy.sa.NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA")

becomes:

meanElev = arcpy.sa.FocalStatistics(filledDEM, arcpy.sa.NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA")

Kind regards,

Xander

[HenriRiihimäki]

Yep, that's right.

Many thanks.

It runs through and everything seems logical here and I think this modified script is safer (however all the honors to Pathak for making the original script availalbe in the first place).  The original problems (drop-values) remains to be solved. The drop-values are still rounded to 16.666 and to 25 % as I mentioned before. With this new python-script I get less unique values, but still values are rounded up. I sent a piece of DEM to ESRI support to check it out. So let's see.

best regards,

Henri