1

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.

1
  • Need to stick to those "basic" modules, what takes so long is the loop over 130.000 ids. Thought about something like parallelization or so. – s6hebern Dec 2 '20 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.