Python - Looping over raster cells, extremely slow

Question

Being not the most savvy python user, I've been trying to find a solution to speed up my code. I have 2 raster files, both have the same dimension and extent. One raster is a river, extracted from a LIDAR image, with a narrow range of elevation values, all other values in the river file are 0. So basically, with the exception of the river width (the pixels with elevation data), 99% of the river raster has the value 0. The other raster file is extracted from a LIDAR image and has elevation data throughout its entire extent.

My plan is to go over each pixel in the river file that does not have a 0 value, add 1 to it (if the river water level rises by 1 meter) and compare it to the neighboring pixel elevation data (basically the pixel location of the LIDAR elevation). If it is higher, give it a value of 1 (which means that this pixel is flooded). The code that I came up with works, but it looks horrible, and extremely slow (paint dries faster than the code runs..) So I'm looking for a way to speed things up. I looked at vectorizing, numpy, etc, but can't make sense of it yet (could be because I've been at work for over 9 hours..)

So here it is, any suggestions would be very much appreciated:

river = arcpy.Raster(r'C:\Flood.gdb\Clip_River1')
lidar = arcpy.Raster(r'C:\Flood.gdb\Clip_Lidar_1')

arrayRiver = arcpy.RasterToNumPyArray(river,nodata_to_value=0)
(rHeight, rWidth)=arrayRiver.shape

arrayLidar = arcpy.RasterToNumPyArray(lidar,nodata_to_value=0)
(lHeight, lWidth)=arrayLidar.shape

# extent of the 2 rasters: columns 3822, rows 10129

flood = 1
while flood == 1:   
    flood = 0
    for row in range(0,rHeight-1):
        for col in range(0,rWidth-1):
            if arrayRiver.item(row,col) <> 0 and arrayRiver[row,col] <> 1:
                if arrayLidar[row,col-1] <> 1 and arrayLidar.item(row,col-1) < arrayRiver.item(row,col)+1:
                    arrayLidar[row,col-1] = 1
                    arrayRiver[row,col-1] = arrayRiver.item(row,col)
                    flood = 1

                [....doing the same concept for each possible neighboring pixel]

                if arrayLidar[row+1,col+1] <> 1 and arrayLidar.item(row+1,col+1) < arrayRiver.item(row,col)+1:
                    arrayLidar[row+1,col+1] = 1 
                    arrayRiver[row+1,col+1] = arrayRiver.item(row,col)
                    flood = 1

            arrayRiver[row,col] = 1

newRaster = arcpy.NumPyArrayToRaster(arrayLidar,lowerLeft,cellSize,value_to_nodata=0)
newRaster.save(r"C:\Flood.gdb\newRastaRiver_small")

newEmptyRaster = arcpy.NumPyArrayToRaster(arrayEmpty,lowerLeft,cellSize,value_to_nodata=0)
newEmptyRaster.save(r"C:\Flood.gdb\newEmptyRastaRiver_small")

Before doing anything you should do a profiling, throw out the whole if statements to see how long your code takes natively, and then you can see where the main bottleneck is. And you didn't add the dimensions of your cells, so its not clear where the bottleneck could be from. — user1767754, Apr 30 '15 at 22:06
I'd ask this on gis.stackexchange.com. In general, don't loop through arrays. As an example, here's a recent question I asked about getting all of the neighbors of each element without a loop: http://stackoverflow.com/questions/29925684/generalize-stacking-of-array-elements-neighbors-into-3-d-array — Paul H, Apr 30 '15 at 22:18
Ok, I removed the entire if statement, and the entire code, including raster to array statement, runs in 8.2 seconds. The dimensions of the raster are 3822 columns, 10129 rows — derBrain, May 01 '15 at 13:21
before I removed it, I had it running for 48 (!) hours, and finally terminated it after the 1530th iteration over the raster.. — derBrain, May 01 '15 at 14:59
@don_q if the pixel at that location is 1, then it is already 'flooded'. I put it in the code so I don't have to re-check the elevation value at arrayLidar[row+1,col+1] — derBrain, May 01 '15 at 15:11
@don_q I guess I should be using != instead of <> from now on to express 'not equal'... Good to know. — derBrain, May 01 '15 at 15:23

don_q · Answer 1 · 2015-04-30T22:21:42.093

1

This is only a partial answer but should get you started.

Use numpy to exclude (mask) all pixels of low values.

arrayRiver = arcpy.RasterToNumPyArray(river,nodata_to_value=0)
# Create a any array makring pixels less than zero.
mask = arrayRiver > 0
arrayRiver = arrayRiver[mask]
arrayLidar = arrayLidar[mask]

Your arrayRiver & arrayLidar should now be much smaller, and the rest of your code should be faster. Btw, your sample code doesn't appear to be making an arrayLidar array.

edited Apr 30 '15 at 22:21

answered Apr 30 '15 at 22:05

don_q

415
3
11

Thanks for the input Don. I tried your suggestion, but it appears that the `arrayRiver = arrayRiver[mask]` line creates a list (please correct me if I'm wrong), and I can't access the original location (row, col) of it anymore to compare it to the lidarArray[row,col] elevation. – derBrain May 01 '15 at 14:54

Python - Looping over raster cells, extremely slow

1 Answers1