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).