I have a slope raster spatial dataset and want to calculate the RUSLE equation's S value. The challenge is that for values with slope gradient <= 9% I have to apply the following equation:
S = 10.8 * Sin(θ) + 0.03 (Equation 1)
And for slope gradient > 9% I have to apply the following equation:
S = 16.8 * Sin(θ) - 0.50 (Equation 2)
Where θ is the angle in PERCENT-RISE of the slope.
So I need to identify within the raster spatial dataset where the slope gradient is 'equal/lower' than 9% and then apply Equation 1, and apply Equation 2 to the locations where the slope is greater than 9% and subsequently apply Equation 2.
So I've written an ArcPy code to try to carry out this task. It works but it is wrong.
Firstly, I "identify" slope gradients equal or lower than 9%, and RECLASSIFY the raster dataset using the 'Con' statement in the 'Raster calculator tool' what is out of this query to "zero", then save it to "outraster3".
Secondly, using the same routine as the above step, I "identify" the slope gradients higher than 9%, RECLASSIFY what is out of this range to "zero" and then save it to "outraster4".
So my problem is that when trying to apply the tool "SET NULL" to transform the values "zero" in both raster datasets to "null" it does not work, but it does not throw an error message either, letting the code run until the end.
Can someone throw some light on this problem?
Please find bellow the code:
# -*- coding: utf-8 -*-
# -*- Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 -*-
"""
Creation start date
Fri Apr 10 16:31:12 2020
Finish date
Wed Apr 21 14:45:12 2020
@author:
Filipe S C Pinhati
PhD candidate
University of Queesnland
School of Civil Engineering
Aquatic studies group
with the help of "BERA" and "Vince"
Geographic Information Systems website
"https://gis.stackexchange.com/questions/357950/arcpy-usle-s-value-set-null-in-raster-is-not-setting-values-to-null"
Description:
Using ArcPy to calculate values for the S parameter
from the RUSLE soil loss equation, according to
slope gradients whether being higher or lower/equal to 9 percent-rise.
The equation applied to calculate the "S" parameter for different
groups of slope is referenced on:
Lu H, Gallant J, Prosser I, Moran C, Priestley G (2001). Prediction
of sheet and rill erosion over Australian continent, incorporating
monthly soil loss distribution. National land & water resource audit.
CSIRO Land and Water, technical report 13/01.
Hrabalikova M, Janecek M (2017). Comparison of different approaches to LS
factor calculations based on a measured soil loss under simulated rainfall.
Soil & Water Res. 12(2):69–77. https://doi.org/10.17221/222/2015-SWR
Requirements: ArcMap
Requirements: Spatial Analyst Extension
"""
"""
Import system modules
"""
import arcpy
from arcpy import env
from arcpy.sa import *
"""
Check out the ArcGIS Spatial Analyst extension license
"""
arcpy.CheckOutExtension("Spatial")
"""
Set environment settings
"""
env.workspace = "C:/User_Folder"
arcpy.env.overwriteOutput = True
"""
Set local variables
"""
print("-"*50)
print(" Reading local variables")
print("-"*50)
inRaster = Raster("SLP_PCT.tif") # slope raster dataset with values in percent-rise
whereClause = "VALUE = 0" # clause to be used in the "SetNull" tool
print("Local variables read\n")
print(" ")
"""
Execute raster calculator to modify the slope raster dataset with values
in percent-rise, keeping the values <= 9% and the remaining pixels set to zero
"""
print("-"*50)
print(" Creating outRas1")
print("-"*50\n)
print("Creating raster dataset with values <= 9% and the remaining pixels set to zero\n")
outRas1 = Con(inRaster <= 9, 1, 0)
outRas1.save("C:/User_Folder/outRas1")
print arcpy.GetMessages()
print("outRas1 done - regions > 9 in raster slope_percent set to zero\n")
"""
Execute raster calculator to modify the slope raster dataset with values
in percent-rise, keeping the values > 9% and the remaining pixels set to zero
"""
print("-"*50)
print(" Creating outRas2")
print("-"*50)
print("Creating raster dataset with values > 9% and the remaining pixels set to zero\n")
outRas2 = Con(inRaster > 9, 1, 0)
outRas2.save("C:/User_Folder/outRas2")
print arcpy.GetMessages()
print("outRas2 done - regions <= 9 in raster slope_percent set to zero\n")
"""
Execute raster calculator to multiply the slope raster dataset with values in percent-rise by "outRas1"
so pixels with angles equivalent to slope <= 9% will be multiplied by 1, while the remaining pixels will be
multiplied by 0
"""
print("-"*50)
print(" Creating outRas3")
print("-"*50)
print("Creating slope raster dataset with pixels in percent-rise in regions where slope is")
print("<= 9% and the remaining pixels set to zero\n")
outRas3 = (inRaster * outRas1)
outRas3.save("C:/User_Folder/outRas3")
print arcpy.GetMessages()
print(" ")
print("outRas3 done - created slope raster dataset in degrees showing areas where the slope is <= 9%")
print("and the remaining areas area set to 0\n")
"""
Execute SetNull to finalize the process of creating the 'slope in percent-rise' raster dataset only in regions
where the slope is equivalent to <= 9% and the remaining regions are set to NO DATA
WARNING! THIS OPERATION WILL SET FLAT AREAS WITH 0° SLOPE TO ZERO
THESE AREAS NEED TO BE ADDED TO FINAL RUSLE`S S AS 0.03 PIXELS VALUE
"""
print("-"*50)
print(" Creating outRas4")
print("-"*50)
print("Applying SetNull to generate slope raster dataset with pixels in degrees in regions where slope is")
print(" <= 9% and the remaining pixels set to NO DATA\n")
outRas4 = SetNull(outRas3,outRas3,whereClause)
outRas4.save("C:/User_Folder/outRas4.tif")
print arcpy.GetMessages()
print(" ")
print("setnull completed as outRas4 - created slope raster dataset with values in percent in areas")
print("equivalent to slope <= 9% and the other areas set to no data\n")
"""
Execute raster calculator to multiply the slope raster dataset with values in percent-rise by "outRas2", so pixels
with angles equivalent to slope > 9% will be multiplied by 1, while the remaining pixels will be multiplied by 0
"""
print("-"*50)
print(" Creating outRas5")
print("-"*50)
print("Creating slope raster dataset with pixels in percent-rise in regions where slope is")
print(" > 9% and the remaining pixels set to zero\n")
outRas5 = (inRaster * outRas2)
outRas5.save("C:/User_Folder/outRas5")
print arcpy.GetMessages()
print(" ")
print("outRas5 done - created slope raster dataset in percent showing areas where the slope is > 9%")
print ("and the remaining areas are set to 0\n")
"""
Execute SetNull to finalize process of creating the 'slope in percent-rise' raster dataset only in the regions
where the slope is equivalent to > 9% while the remaining regions are set to NO DATA
"""
print("-"*50)
print(" Creating outRas6")
print("-"*50)
print("Applying SetNull to generrate slope raster dataset with pixels in percent-rise in regions where slope is")
print(" > 9% and the remaining pixels set to NO DATA\n")
outRas6 = SetNull(outRas5,outRas5,whereClause)
outRas6.save("C:/User_Folder/outRas6.tif")
print arcpy.GetMessages()
print(" ")
print("setnull completed as outRas6 - created slope raster dataset with values in percent in areas")
print("equivalent to slope > 9% and the other areas set to no data\n")
"""
--- note for generating outRas7 and outRas8 ---
Execute 'Math Toolset' to apply equation to calculate S values for slope areas of
<= 9% (outRas7) and >9% (outRas8).
The requested input to calculate the RUSLE`s S value must be slope in perecent-rise.
Hence, firstly the pixel values must be converted from percent-rise to radians using the formula:
angle_radians = atan*(angle_percent/100)
In ArcPy the input of the 'ATan' tool is interpreted to be unitless, and its output is in radians.
As the input values for the 'Sin" tool are interpreted to be in radians, no further conversion is needed.
In ArcPy the output of the 'Sin' tool is unitless
"""
"""
Calculating 'S' for outRas7 = (10.8*(sin(outRas4))+0.03)
"""
print("-"*50)
print(" Creating outRas7")
print("-"*50)
print("Calculating S values for areas with slope <=9%\n")
outDivide1 = Divide(outRas4,100) # first step to convert angle in percent-rise to radians
outATan1 = ATan(outDivide1) # finalizing the process of conversion from percent-rise to radians
outSin1 = Sin(outATan1) # initializing the S value calculation
outTimes1 = Times(10.8, outSin1) # calculating the second term of the S equation
outRas7 = Plus(0.03,outTimes1) # calculating the third term of the S equation
outRas7.save("C:/User_Folder/outRas7")
print arcpy.GetMessages()
print(" ")
print("outRas7 done - RUSLE S factor for angle <= 9% calculated\n")
"""
Calculating 'S' for outRas8 = (16.8*(sin(outRas6))-0.5)
"""
print("-"*50)
print(" Creating outRas8")
print("-"*50)
print("Calculating S values for areas with slope >9%\n")
outDivide2 = Divide(outRas6,100) # first step to convert angle in percent-rise to radians
outATan2 = ATan(outDivide2) # finalizing the process of conversion from percent-rise to radians
outSin2 = Sin(outATan2) # it results in the same outcome as calculating the sine of an angle in degrees
outTimes2 = Times(16.8,outSin2) #calculating the third term of the S equation
outRas8 = Minus(outTimes2,0.5)
outRas8.save("User_Folder/outRas8")
print arcpy.GetMessages()
print(" ")
print("outRas8 done - RUSLE S factor for angle > 9% calculated\n")
"""
Mosaic outRas7 and outRas8 creating a final raster showing "S" values for the USLE equation
"""
print("-"*50)
print(" Creating RUSLE S parameter")
print("-"*50)
print("Mosaicking S values raster datasets\n")
try:
arcpy.MosaicToNewRaster_management("outRas7;outRas8","C:/User_Folder","RUSLE_S", "",\
"", "", "1", "","")
print arcpy.GetMessages()
print("-"*50)
print(" S PARAMETER GENERATED"*1)
print("-"*50)
except:
print "Mosaicking outRas7 & outRas8 failed\n"
print arcpy.GetMessages()
