1

I am very new to python coding and have been struggling with this for a while I cannot see the reason for this error!

Here is the code that is causing this error:

import arcpy
import numpy as np
import math
from arcpy.sa import *

# Get the input arguements

# Input DEM
input_raster = arcpy.GetParameterAsText(0)

# Output 
output_raster = arcpy.GetParameterAsText(1)

# Input Discharge
mean_discharge = float(arcpy.GetParameterAsText(2))

# Input number of timesteps
input_timestep = int(arcpy.GetParameterAsText(3))

#Other stuff to put outside while loop

# The below text takes the input raster and calculates the bottom left corner
extent_xmin_result = arcpy.GetRasterProperties_management(input_raster, "LEFT")
extent_xmin = float(extent_xmin_result.getOutput(0))
extent_ymin_result = arcpy.GetRasterProperties_management(input_raster, "BOTTOM")
extent_ymin = float(extent_ymin_result.getOutput(0))

# Turns the corner into a point
bottom_left_corner = arcpy.Point(extent_xmin, extent_ymin)

arcpy.overwriteOutput = True

arcpy.env.workspace = r"C:\SedMod\SedMod\Workspace.gdb"
arcpy.env.scratchWorkspace = r"C:\SedMod\SedMod\Workspace.gdb"

actual_timestep = 0

while actual_timestep <= input_timestep: # Slope calculation

# Set local variables
out_measurement = "PERCENT_RISE" 

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

#Execute Slope
out_slope = Slope(input_raster, out_measurement)

# Convert to gradient
slope_percentage = out_slope / 100

#Save slope
save_location = r"C:/SedMod/SedMod/Workspace.gdb/test" + str(actual_timestep)
slope_percentage.save(save_location)

#Execute Flow Direction

output_flowdirection = FlowDirection(input_raster)

#Execute Drainage Area

# Method for describing the raster and then using the cellsize to determine the Drainage factor for calculating Drainage area
desc_raster = arcpy.Describe(input_raster)
raster_cell_width = desc_raster.MeanCellWidth # Cell width
raster_cell_height = desc_raster.MeanCellHeight # Cell height
spatial_reference = desc_raster.spatialReference # Spatial reference
multi_factor = (raster_cell_height / 1000) * (raster_cell_width / 1000)

# Run the flow direction
output_flowaccumulation = FlowAccumulation(output_flowdirection)
output_drainage_area = output_flowaccumulation * multi_factor

#Execute Discharge Calculation

mean_discharge_float = float(mean_discharge)

#Calculating the discharge based on the average taken from the Avon catchment to bath (1605km2) and a mean flow of 20.762 m3/s giving a discharge per km2 of 0.01294 m3/s
discharge_raster = output_drainage_area * mean_discharge_float

#Execute StreamPower

# Set variables for calculating stream power    
density_of_water = 1000
acceleration_gravity = 9.81

# Stream power calculation = stream power, g is the density of water (1000 kg/m3), g is acceleration due to gravity (9.8 m/s2), discharge (m3/s) and channel slope.
unit_width_stream_power = (density_of_water * acceleration_gravity * discharge_raster * slope_percentage) / raster_cell_width

#Execute calculating Sediment Transport        

# Set variables for calculating dimensionless stream power    
pw = 1000 # Density of water
ps = 2650 # Density of sediment
g = 9.81 # Acceleration due to gravity 
Di = 0.01 # Sediment size

# Dimensionless unit width stream power
dimensionless_stream_power = unit_width_stream_power / ((16186.5) * math.sqrt((1.65 * 9.81 * (Di * Di *Di ))))

# Additional variables to calculating sediment transport
timestep_hour = 3600

# Calculates dimensionless sediment transport       
dimensionless_sediment_transport = Con(dimensionless_stream_power < 0.25, ((dimensionless_stream_power ** 6) * 100), (0.2 * (dimensionless_stream_power ** 1.5)))

# Calculates bed transport rate
bed_material_transport = dimensionless_sediment_transport * ((16186.5) * math.sqrt((1.65 * 9.81 * (Di * Di *Di ))))            

# Calculates it for the cell size in kg per second        
bed_material_transport_kg = (bed_material_transport * raster_cell_width)  / 9.81

# Converts it into volume m3
bed_material_transport_volume = (bed_material_transport_kg / 2650) * timestep_hour



# Convert the sediment transport and the flow direction rasters into Numpy arrays
sediment_transport_np = arcpy.RasterToNumPyArray(bed_material_transport_volume, '#', '#', '#', -9999)
flow_direction_np = arcpy.RasterToNumPyArray(output_flowdirection, '#', '#', '#', -9999)

[rows,cols] = flow_direction_np.shape
elevation_gain = np.zeros((rows,cols), np.float)
elevation_gain_32 = np.zeros((rows,cols), np.float)
elevation_gain_64 = np.zeros((rows,cols), np.float)
elevation_gain_128 = np.zeros((rows,cols), np.float)
elevation_gain_16 = np.zeros((rows,cols), np.float)
elevation_gain_1 = np.zeros((rows,cols), np.float)
elevation_gain_2 = np.zeros((rows,cols), np.float)
elevation_gain_4 = np.zeros((rows,cols), np.float)
elevation_gain_8 = np.zeros((rows,cols), np.float)

# Main body for calculating elevation change        

# Attempt 4

# Move the sediment transport in the direction for 32 direction        
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 32:
            elevation_gain_32[i - 1, j - 1]  = sediment_transport_np[i, j]
    except IndexError:
            elevation_gain_32[i - 1, j - 1] = 0

# Move the sediment transport in the direction for 64 direction
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 64:
            elevation_gain_64[i - 1, j]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_64[i, j] = 0

# Move the sediment transport in the direction for 128 direction                  
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 128:
            elevation_gain_128[i - 1, j + 1]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_128[i, j] = 0                

# Move the sediment transport in the direction for 16 direction
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 16:
            elevation_gain_16[i, j - 1]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_16[i, j] = 0

# Move the sediment transport in the direction for 1 direction
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 1:
            elevation_gain_1[i, j + 1]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_1[i, j] = 0

# Move the sediment transport in the direction for 2 direction
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 2:
            elevation_gain_2[i + 1, j + 1]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_2[i, j] = 0

# Move the sediment transport in the direction for 4 direction
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 4:
            elevation_gain_4[i + 1, j]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_4[i, j] = 0

# Move the sediment transport in the direction for 8 direction
for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 8:
            elevation_gain_8[i + 1, j - 1]  = sediment_transport_np[i, j]
    except IndexError:
        elevation_gain_8[i, j - 1] = 0                                   

# Covert numpy arrays back into ArcGIS raster format as I do not know how to successfully add together multidimensional numpy arrays                  
elevation_gain_32_ras = arcpy.NumPyArrayToRaster(elevation_gain_32, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_64_ras = arcpy.NumPyArrayToRaster(elevation_gain_64, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_128_ras = arcpy.NumPyArrayToRaster(elevation_gain_128, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_16_ras = arcpy.NumPyArrayToRaster(elevation_gain_16, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_1_ras = arcpy.NumPyArrayToRaster(elevation_gain_1, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_2_ras = arcpy.NumPyArrayToRaster(elevation_gain_2, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_4_ras = arcpy.NumPyArrayToRaster(elevation_gain_4, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
elevation_gain_8_ras = arcpy.NumPyArrayToRaster(elevation_gain_8, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)

# Calculate  volume in raster for whole catchment        
volume_in = elevation_gain_32_ras + elevation_gain_64_ras + elevation_gain_128_ras + elevation_gain_16_ras + elevation_gain_1_ras + elevation_gain_2_ras + elevation_gain_4_ras + elevation_gain_8_ras

#State what volume out is and the elevation           

net_volume = volume_in - bed_material_transport_volume

elevation_change = net_volume / (raster_cell_width * raster_cell_height)

new_elevation = input_raster + elevation_change

actual_timestep = actual_timestep + 1

new_elevation.save(output_raster)   

The error I am getting is:

Traceback (most recent call last):
  File "C:\SedMod\SedMod\Slope1.py", line 47, in <module>
    out_slope = Slope(input_raster, out_measurement, )
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 5523, in Slope
    z_factor)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 5518, in wrapper
    z_factor)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: ERROR 010240: Could not save raster dataset to C:\SedMod\SedMod\Workspace.gdb\numpy_ras6 with output format FGDBR.

Failed to execute (SedModv1).
1
  • 1
    Do input_timestep and mean_discharg need to be integers? arcpy.GetParameterAsText returns the value as a string. I would guess the while loop is looking like this : while 0 < '1', which doesn't make sense. Maybe i'm wrong, but that's my guess.
    – dchaboya
    Commented Jun 17, 2013 at 17:59

1 Answer 1

1

It looks like you are telling Arc that you are trying to input a z-factor in the slope tool, but not actually defining the variable. Removing the extra comma from that line of code should help, unless of course you do want to specify the variable in which case you need to put something in there after the comma:

out_slope = Slope(input_raster, out_measurement )

or

out_slope = Slope(input_raster, out_measurement, slope )

As @dchaboya has pointed out, you'll need to explicitly cast the parameter-as-text variables as integers or your script may not work as expected:

mean_discharge = int(arcpy.GetParameterAsText(2))
input_timestep = int(arcpy.GetParameterAsText(3))
9
  • I thought this too. I am afraid this is not a fix still getting the error on the 2nd loop round.
    – Nick Jones
    Commented Jun 17, 2013 at 17:42
  • Ok, the last bit of the error where is says Could not save raster dataset to C:\SedMod\SedMod\Workspace.gdb\numpy_ras6 with output format FGDBR looks to be signalling the fact that you are trying to save multiple rasters to the same GDB with the same name ("PERCENT_RISE"). See help.arcgis.com/EN/arcgisdesktop/10.0/help/index.html#//…. You must either devise a way to come up with unique names (e.g. add your counter-value to the name using string formatting) or set acrpy.overwriteOutput = True. Commented Jun 17, 2013 at 17:53
  • Thanks for the suggestions I have implemented those changes too, but still not working :( I will edit the updated code above.
    – Nick Jones
    Commented Jun 17, 2013 at 20:08
  • 1
    I didn't read the documentation on the slope tool well enough...looks like the output_measurement variable you had was good, needs to be one of "DEGREE" or "PERCENT RISE". Sorry to confuse you there. Not sure where the problem is - can you update your question with your latest code? I will edit my answer to remove the part where it changes the output_measurement variable to a non-acceptable input value for the slope tool. Commented Jun 17, 2013 at 20:46
  • 1
    Are still getting exactly the same error message? Commented Jun 17, 2013 at 21:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.