# Efficient way to analyze multiple (large) rasters against a different raster using numpy

The situation I have run into involves four rasters. I will try to explain each one and what it is for to help you get a better idea of what I'm trying to achieve.

One of the rasters has been created using gdal rasterizeLayers on vector polygon data. The polygon data represented swath polygons from agriculture machinery. The way I rasterized the vector data was to give each individual swath polygon (which there are just over 80,000 of) a unique value (just 1-80k starting with the first polygon and ending with the last, incrementing by 1). I will refer to this raster as the "swath" raster.

One of the other three rasters is the result of running gdals IDW interpolation on point data that also came from an agriculture machine (Harvester). This raster has values between 1-350 assigned to each cell. Some cells are of course marked with nodata where no crop was harvested.

The next of the three rasters is very similar to the previous one, except it was from a planter instead of a harvester. This one recorded population values (the number of seeds planted at every location the machine logged a point). The values in this raster are 25-40k.

The final of the three rasters is the result of rasterizing a vector polygon map that represented soils data. There were three different soil types in the polygons, so each cell in this raster is either a 1, 2, or 3 depending on which soil type it represents.

All of these rasters are covering the same geospatial extent (a field where crops were planted and harvested using precision ag machinery). I have also cropped each of the rasters so that they use the same number of columns, rows, and pixel sizes. This way they will line up properly.

The general idea of what I'm trying to do is extract data from the three rasters (yield, planting, soils) into each individual swath polygon from the swath raster. I have a working implementation, but it takes upwards of 30 minutes per layer. I am hoping to find a way to achieve the same result in much less time.

The way I have gone about doing this is to load the raster data into numpy arrays for the swath raster and one other raster at a time. Once I've done that, I simply run a loop (starting with the lowest value in the swath raster and ending with the highest) and create a mask using the swath raster that will mask all pixels/cells except the ones that contain the value of the current iteration of the loop. I then apply that mask to the other raster to find the pixels/cells that are "within" or "overlapping" the swath polygon from the swath raster. In the case of the planting and yield rasters, I am getting the mean of all unmasked values. In the case of the soils raster, I am getting the mode of all unmasked values.

Aside from the method that I currently have working, I have also tried using the "raster" package in R to run zonal statistics (with the swath polygon layer being vector and all others still being raster). This worked, but it took about 4-5 times as long as the method I'm currently using.

Here is the code example for what I am currently doing. I am using python, python gdal bindings, numpy, and scipy. This example is only the function that I'm calling with the data which performs the actual masking/analysis. I call this function once with the swath polygon raster and each of the three other rasters (yield, soils, and planting). This function is where nearly all of the time is taken. I think it is taking so long because there are so many different unique values in the swath raster to iterate over. I'm hoping there is a faster way to achieve the same thing.

``````def masking(results, polygons, swath, layer, result_index):
# results is just a dictionary with key / values to store the results per unique swath polygon
# polygons is just an integer representing the number of unique values in the swath polygons raster
# swath is a class I've created that has the swath raster datasource and raster band loaded
# layer will be a class (identical to swath class, just helps load rasters) that has the soils, planting, or yield datasource and raster band loaded
# result index is a string that will be used as the key when we store the mean/mode in the results dictionary

# the swath polygon raster band
# the yield, planting, or soils raster band
### time - 0.0400419235229

# iterate each unique value in the swath raster
### start timer to time entire loop (1000 iterations, then break for timing run)
for i in range(1, polygons + 1):
# get some information from the swath polygon raster attribute table (RAT) using the current pixel value (current iteration value of i)
### start timer to time getting swath rat info
[value, product, row_spacing] = rats.findSwathPolygonInformation(swath.rat, i)
### end timer to time getting swath rat info
### time - 4.81605529785e-05

# create an entry in the dictionary for this value (i) if one doesn't exist
### start timer to time adding element to dictionary
if i not in results:
results[i] = {'product': product, 'id': value, 'row_spacing': row_spacing}
### end timer to time adding element to dictionary
### time - 9.53674316406e-07

# mask all pixels/cells that do not contain the value of the current iteration (i)
### start timer to time first masking operation
### end timer to time first masking operation
### time - 0.0231530666351

# use the previous mask to mask the other layer (soils, planting, or yield)
### start timer to time second masking operation
### end timer to time second masking operation
### time - 0.000151872634888

# remove any zeroes
### start timer to time third masking operation
### end timer to time third masking operation
### time - 0.0105729103088

# get the mode if we have soils
if result_index == 'soils':
# get mode
### start mode timer
### end mode timer
### time - 0.00011682510376

# if we got a value
if len(mode):
### start if block timer
# get the value from the array
mode = int(mode[0])
# look the value up in the soils rat and get a string (name of soil type)
soil = rats.findSoilType(layer.rat, mode)
# set the value on the current element in results
results[i][result_index] = soil
### end if block timer
### time - 1.59740447998e-05
else:
# there was no soils data for this swath polygon
results[i][result_index] = 0
else:
# get the mean and set it on the current element in results
### start mean timer
### end mean timer
### time - 6.5803527832e-05

# make sure we don't get any bad values (None, etc)
### start isinstance timer
if not isinstance(results[i][result_index], float):
results[i][result_index] = 0
### end isinstance timer
### time - 1.19209289551e-06

### end loop timer
### time - 35.7992780209

return results
``````

This is my first time posting, so please let me know if you need any more information. I'm mostly just looking for suggestions on what I might be able to do in order to make this faster. I am open to using other tools if needed. I can also provide a full working code example and test data if needed.

UPDATE : I have done some timing to determine where most of the time is being spent. I tested this while running 1000 iterations of the loop and then breaking. The way I did the timing was simply importing `time` and then storing `time.time()` in a variable at the start point, then storing `time.time()` again at the stopping point, and finally subtracting the stopping point time from the starting point time. If it would be better to do this another way, let me know what to do. I have inserted timer specific comments in the code example with the start/stop points and the time it took from start to stop. The timer related comments start with ###. It looks like the three masking options are taking the most time compared to any of the other operations within the loop. They don't take a long time in any single iteration (I added the times for each iteration and divided by 1000 (total iterations before breaking) to get the times for those), but we are iterating just over 80000 times in this case.

• Have you benchmarked your code to see what sections are taking the longest? – Paul Jan 13 '17 at 21:04
• I will do that now to confirm exactly where the problem is. I'll update the post with the answer. Thanks for the suggestion! – Dylan Hamilton Jan 13 '17 at 21:11
• What are the dimensions of the rasters? – Logan Byers Jan 13 '17 at 22:10
• I am not sure what you hope to accomplish. You are going to come up with some summary statistics for each "swath" and try to determine what factors are affecting yields? I think you should have a look at the tools in Spatial Analyst>Multivariate and see if you could use that as a starting point – Mike Jan 13 '17 at 22:13
• Did you have a look at the rasterstats package? I think it does exactly what you want to achieve – Loïc Dutrieux Jan 14 '17 at 11:01