Finding maximum values (of one raster) which are inside each segment of another raster

Question

There are two rasters as below. One consisting of only four values [1,2,3,4]. The other, consisting of values between 800 to 2500. The problem is to go through all of the raster-1 regions and find the maximum values of raster-2 which are located inside each region or segment.

In theory, it seems simple but I can't find a way to implement it. I'm reading scikit image documentation and I'm getting more confused. In theory, it would be:

for i in raster1rows:
    for j in i:
        # where j is a part of closed patch, iterate through the identical
        # elements of raster-2 and find the maximum value.

There is another problem inherent to this question which I can't post as a different topic. As you can see, there are a lot of isolated pixels on raster-1, which could be interpreted as a region and produce a lot of additional maximums. to prevent this I used :

raster1 = raster1.astype(int)
raster1 = skimage.morphology.remove_small_objects(raster1 , min_size=20, connectivity=2, in_place=True)

But raster-1 seems to take no effect.

Answer 1

To remove the small object I've done

array_aspect = sp.median_filter(array_aspect, size=10)

And it gave me good results.

To find maximum elevation inside each closed part I've done:

# %%% to flood-fill closed boundaries on the classified raster
p = 5
    ind = 1
    for i in rangerow:
        for j in rangecol:
            if array_aspect[i][j] in [0, 1, 2, 3, 4]:
                print("{}. row: {} col: {} is {} is floodfilled with {}, {} meters".format(ind, i, j, array_aspect[i][j], p, array_dem[i][j]))
                array_aspect = sk.flood_fill(array_aspect, (i,j), p, in_place=True, connectivity=2)
                p = p + 1
            else:
                pass
            ind = ind + 1
# %%% Finds the max elev inside each fill and returns an array-based [Y,X, (ELEV #in meters)]
p = 5
maxdems = {}
for i in rangerow:
    for j in rangecol:
        try:
            if bool(maxdems[array_aspect[i][j]]) == False or maxdems[array_aspect[i][j]][-1] < array_dem[i][j]:
                maxdems[array_aspect[i][j]] = [i, j, array_dem[i][j]]
            else: 
                pass
        except: #This is very diabolical, but yeah :))
                maxdems[array_aspect[i][j]] = [i, j, array_dem[i][j]]
print(maxdems)`

I've got my desired results.