Iterating through an numpy array and index to a value in another numpy array

Question

I am struggling to get this code to work I want to iterate through an numpy array and based on the result, index to a value in another numpy array and then save that in a new position based on that value.

 # Convert the sediment transport and the flow direction rasters into Numpy arrays
 sediment_transport_np = arcpy.RasterToNumPyArray(sediment_transport_convert, '#', '#', '#', -9999)
 flow_direction_np = arcpy.RasterToNumPyArray(flow_direction_convert, '#', '#', '#', -9999)
 [rows,cols]= sediment_transport_np.shape
 elevation_change = np.zeros((rows,cols), np.float)
 # Main body for calculating elevation change 
 # Attempt 1
 for [i, j], flow in np.ndenumerate(flow_direction_np):
 if flow == 32:
 elevation_change[i, j] = sediment_transport_np[i - 1, j - 1]
 elif flow == 16:
 elevation_change[i, j] = sediment_transport_np[i, j - 1]
 elif flow == 8:
 elevation_change[i, j] = sediment_transport_np[i + 1, j - 1]
 elif flow == 4:
 elevation_change[i, j] = sediment_transport_np[i + 1, j]
 elif flow == 64:
 elevation_change[i, j] = sediment_transport_np[i - 1, j]
 elif flow == 128:
 elevation_change[i, j] = sediment_transport_np[i - 1, j + 1]
 elif flow == 1:
 elevation_change[i, j] = sediment_transport_np[i, j + 1]
 elif flow == 2:
 elevation_change[i, j] = sediment_transport_np[i + 1, j + 1]
 #sediment_flow_np = np.concatenate((flow_direction_np, sediment_transport_np), axis = 0)
 #case1 = numpy.where((flow_direction_np == 32)) & (sediment_transport_np > 0), 
 # Convert the Numpy arrays back to rasters
 #sediment_transport_raster = arcpy.NumPyArrayToRaster(sediment_transport_np, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
 #flow_direction_raster = arcpy.NumPyArrayToRaster(flow_direction_np, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
 elevation_change_raster = arcpy.NumPyArrayToRaster(elevation_change, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
 elevation_change_raster.save(output_raster)

The error i get is:

Running script elevation_change...

Traceback (most recent call last): File "", line 606, in execute IndexError: index (655) out of range (0<=index<655) in dimension 0

Failed to execute (elevation_change)

• Jason Bellino: It looks like your flow direction array has different dimensions than one or more of the other two arrays which you are indexing with i and/or j.
• Jenny Jones: They should all be identical is there anyway to check?
• Jason Bellino: I see you've used the shape property of the Use the sediment_transport_np array to create another array of the same dimensions, flow_direction_np. Use the shape property of the flow_direction_np array to make sure it has the same dimensions.

Answer 1

This error means that you are trying to get data from the 656th row of the sediment_transport array which is only 655 rows long. This is probably happening at the outlet of the flow direction raster. At this spot, the flow direction raster is saying something like "water will flow south from this cell." When you try to index the cell in sediment_transport that is one cell south, it doesn't exist and is throwing the error.

To get around this, you could just add a try and except and fill elevation_change with some placeholder value like -1:

for [i, j], flow in np.ndenumerate(flow_direction_np):
 try:
 if flow == 32:
 elevation_change[i, j] = sediment_transport_np[i - 1, j - 1]
 elif flow == 16:
 elevation_change[i, j] = sediment_transport_np[i, j - 1]
 elif flow == 8:
 elevation_change[i, j] = sediment_transport_np[i + 1, j - 1]
 elif flow == 4:
 elevation_change[i, j] = sediment_transport_np[i + 1, j]
 elif flow == 64:
 elevation_change[i, j] = sediment_transport_np[i - 1, j]
 elif flow == 128:
 elevation_change[i, j] = sediment_transport_np[i - 1, j + 1]
 elif flow == 1:
 elevation_change[i, j] = sediment_transport_np[i, j + 1]
 elif flow == 2:
 elevation_change[i, j] = sediment_transport_np[i + 1, j + 1]
 except IndexError:
 elevation_change[i, j] = -1 #Placeholder for flow direction outlet

I also want to point out that if you are trying to do a flow accumulation of the sediment_transport data, the Flow Accumulation geoprocessing tool in ArcGIS has a parameter for an optional weight raster. The result would have the accumulated values from the in_weight_raster instead of just a count of upstream cells.