current community

your communities

more stack exchange communities

Stack Exchange
tour help
Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 3 down vote favorite

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)

share|improve this question
1  
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. –  Jason Jun 14 '13 at 17:47
    
They should all be identical is there anyway to check? –  Jenny Jones Jun 14 '13 at 18:26
    
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. –  Jason Jun 14 '13 at 19:35

1 Answer 1

up vote 3 down vote

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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