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

arcpynumpypythonraster

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)

Best Answer

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.

Related Question