I have two rasters, one is a classified image and the other the corresponding probabilities for a correct classification. Now I need to calculate zonal statistics for polygon segments in the following way:
- get majority of predicted class within each segment
- get all probability values within this segment ONLY for those pixels that have the majority class and calculate the median probability
Here is what I do (basic code only):
- rasterize polygons based on FID to create an FID-array
- read classification array
- read probability array
- create new fields in the polygon file that shall hold the final values
# id_array: holds all FIDs as pixel values # classi_array: holds classification (per pixel) # prob_array: holds probabilities (per pixel) # lyr: Layer object from ogr.DataSource().GetLayer() # all_ids: array of unique FIDs created with np.unique(id_array[id_array != NODATA]) # this is the slow part: for fid in all_ids: segment = np.where(id_array==fid) majority_class = np.bincount(classi_array[segment].flatten()).argmax() probs = prob_array[segment] majority_probs = probs[np.where(classi_array[segment] == majority_class)] med_majority_probs = np.nanmedian(majority_probs) feat = lyr.GetFeature(fid) feat.SetField(class_column, int(majority_class)) feat.SetField(prob_column, float(med_majority_probs)) lyr.SetFeature(feat) feat.Destroy()
It works, but takes forever. The question is: how can I speed up this? Is there any way to use something like multiprocessing here (with which I've never worked before) or even something completely different?
Versions (I'm bound to them due to company-wide installations, sorry):
- Python 2.7
- numpy 1.15.4
- gdal 1.11.3.
EDIT: Now that it has run, here is the performance for the loop:
Progress: 100%|██████████| 130833/130833 [5:43:43<00:00, 6.34it/s]
EDIT 2: I guess I could somehow utilize the answer from https://stackoverflow.com/questions/20190668/multiprocessing-a-for-loop, but just can't get it to work.