I have a very large raster data and a shape file with 79,867 polygons. I want to efficiently extract statistics from the raster within each of the polygons. To do so, I have combined rasterio and geopandas. I first convert the shape into a geopandas dataframe. Then I iterate over each row of the dataframe and cookie cut the raster with each of the rows (as one-row geopandas dataframe), generate a numpy array and then I can get any stats I want with it (see code below). I save the stats into a dictionary and convert the final dictionary into a pandas DF. This is a slow process, of course, because iterating over dataframe rows is naturally slow. One feature of my code below is that I have a provision to save partial results into a csv file and give me a progress report. My goal now is to parallelize this algorithm. I have been searching alternatives for that, and most of the suggestions include dividing the dataframe into chunks and process each chunk in a vectorized fashion in a different core. This can be done with multiprocessing pool and there is plenty of examples available. However, if I do that, I will lose (at least as far as my coding skills go) the ability to create progress reports and save partial results. In addition, because the cookie cut process is rather slow, I am not convinced that vectorizing it will show gains in speed, and I am not sure the geopandas dataframe will handle such a complex process well. So, my inclination is to retain the row by row processing and distribute the individual rows into different cores instead of partitioning the entire dataframe as suggested. I would like some advice on what is the best alternative to parallelize this process. Unfortunately, most examples I came across to parallelize pandas dataframe processes apply very simple functions, like x*x, which is far from being this case.
import pandas as pd
import geopandas as gpd
import rasterio as rio
from rasterio.mask import mask
from time import time
StartTime = time()
geodf = gpd.read_file("data/polygon.shp")
statsDict = {}
with rio.open("data/raster.tif") as src:
profile = src.profile
n=0
for gpd in geodf.itertuples():
if n > 0 and n%5000 == 0:
perc = 100*n/geodf.shape[0]
print (int(perc), ' %', (time() - StartTime)/ 60, ' minutes')
forecastFinish = round(((100 * (time() - StartTime)/perc)/3600), 1)
print(' estimated time to finish: %s hrs' %forecastFinish)
df = pd.DataFrame(statsDict.items(), columns=['feature_id', 'meanScore'])
fileName = 'results/partial_' + str(n) + '.csv'
df.to_csv(fileName, index=False)
statsDict = {}
feature_id = (gpd.__getattribute__('feature_id'))
polygon = gpd.__getattribute__('geometry')
IterPolygon = [gpd.__getattribute__('geometry')]
mask_feature_id, mask_feature_id_Transform = mask(xscores, \
IterPolygon, invert=False)
flat_mask_feature_id = np.ravel(mask_feature_id)
meanScore = np.mean(convertNpArrayIntoPreppedList(flat_mask_feature_id))
statsDict[feature_id] = meanScore
n+=1
df = pd.DataFrame(statsDict.items(), columns=['feature_id', 'meanScore'])
fileName = 'results/final' + '.csv'
df.to_csv(fileName, index=False)
