Locating the maximum values of a raster intersection with a vector of polygons using Python

The rasterstats python module is based in rasterio, however, it has crop = True option activated by default (when I use raster_out = True for getting all unmasked array values) and I couldn't find an easy way for accessing to raster indices for converting them as points in map units. So, I used directly rasterio tools in following script (combined with your code lines for comparison purpose) and it was possible to identify the coordinates of the found maximum values.

import fiona
import json
import numpy as np
import rasterio
import rasterio.mask
from shapely.geometry import box, Polygon
from rasterstats import zonal_stats

def max_round(x):
    m = x.max()
    if m is np.ma.masked:
        return None
    else:
        return m

vectorbac = "/home/zeito/pyqgis_data/polygon9.shp"
raster = "/home/zeito/pyqgis_data/raster3.tif"

stats = zonal_stats(vectorbac, 
                    raster,
                    copy_properties = True, 
                    geojson_out = True, 
                    prefix = "preci_",
                    stats = ['max'],
                    raster_out = True, 
                    add_stats = {"max": max_round})

loopy = range(0, len(stats))

geom = []

for i in loopy:
    print(stats[i]['properties']['id'], 
          stats[i]['properties']['preci_max'])
    geom.append(stats[i]['geometry']['coordinates'])

with fiona.open("/home/zeito/pyqgis_data/polygon9.shp", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

with rasterio.open("/home/zeito/pyqgis_data/raster3.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, 
                                                  shapes, 
                                                  crop = False, 
                                                  nodata = -999 )
    out_meta = src.meta

xminR = out_transform[2]
xmaxR = xminR + out_meta['width']*out_transform[0]
ymaxR = out_transform[5]
yminR = ymaxR + out_meta['height']*out_transform[4]

for i in range(len(geom)):
    poly = Polygon(geom[i][0])
    bounds = poly.bounds
    xmin, ymin, xmax, ymax = bounds
    rowR1 = int((ymaxR - ymax)/(-out_transform[4]))
    colR1 = int((xmin - xminR)/out_transform[0])
    rowR2 = int((ymaxR - ymin)/(-out_transform[4]))
    colR2 = int((xmax - xminR)/out_transform[0])
    width = colR2 - colR1 + 1
    height = rowR2 - rowR1 + 1

    print(rowR1, colR1, rowR2, colR2, width, height)

    values = []
    indices = []

    for i in range(rowR1, rowR2):
        for j in range(colR1, colR2):
            values.append(out_image[0][i][j])
            indices.append([i,j])

    ma = np.ma.masked_equal(values, -999, copy=False)

    print(ma.max())
    idx = ma.argmax()
    print(indices[idx])
    print(xminR + indices[idx][1]*out_transform[0] + out_transform[0]/2, ymaxR + indices[idx][0]*out_transform[4] + out_transform[4]/2)

After running above script with layers of following image:

it was printed in Python Console following values:

0 2097  #first maximum with rasterstats
1 2000  #second maximum with rasterstats
2 2004  #third maximum with rasterstats
2 2 9 13 12 8
2097  #first maximum with rasterio
[3, 3] #indices for first max
457704.6911268293 4455264.170870967  #first point with rasterio
0 22 8 35 14 9
2000 #second maximum with rasterio
[4, 22] #indices for second max
459110.47684390243 4455190.182148387  #second point with rasterio
12 15 24 32 18 13
2004 #third maximum with rasterio
[19, 16]  #indices for third max
458666.54451219516 4454080.351309677  #third point with rasterio

With the help of QuickWKT plugin in QGIS 3, I generated point layers observed in above image and with Value Tool plugin it was corroborated that all points are effectively maximum points.

Editing Note:

I spent a lot of time checking each masked arrays produced with raster_out = True option of rasterstats python module and, despite it produces correct maximum values and indices, they are in a wrong position in masked arrays. So, it is impossible to identify the exact coordinates of the found maximum values with rasterstats module. I think that there is probably a bug in it.

rasterstats.zonal_stats() also returns affine transformation (mini_raster_affine) of each masked array (mini_raster_array), thus enabling us to locate XY position of the maximum values.

import geopandas as gpd
import pandas as pd
import numpy as np
from rasterstats import zonal_stats

vectorbac = "/home/zeito/pyqgis_data/polygon9.shp"
raster = "/home/zeito/pyqgis_data/raster3.tif"

def pixel2coord(col, row):
    """Returns global coordinates to pixel center using base-0 raster index"""
    xp = a * col + b * row + a * 0.5 + b * 0.5 + c
    yp = d * col + e * row + d * 0.5 + e * 0.5 + f
    return(xp, yp)

stats = zonal_stats(vectorbac, raster , stats = ['max'], raster_out = True, geojson_out = True, nodata=-999)

df = pd.DataFrame(columns=list('XY'))
for idx, i in enumerate(range(0, len(stats))):
    arr = stats[i]['properties']['mini_raster_array']
    col = np.argmax(np.max(arr, axis=0))
    row = np.argmax(np.max(arr, axis=1))
    c, a, b, f, d, e = stats[i]['properties']['mini_raster_affine'].to_gdal()
    df.loc[idx] = pixel2coord(col, row)

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.X, df.Y))