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

I have extracted the statistics of a raster in the polygons of a vector layer using the python `zonal_stats` function. I need to identify the coordinates of the found maximum values. How can I do it using Python?

This is the Python script which I used to extract zonal statistics. I need to retrieve the location of `'preci_max'` for each polygon.

``````import json
import numpy as np
from rasterstats import zonal_stats

def max_round(x):
m = x.max()
return None
else:
return m

vectorbac = "bacini.json"
raster = "preci.nc"

stats = zonal_stats(
vectorbac, raster,
copy_properties=True, geojson_out=True, prefix="preci_",
"max": max_round,
}
)

loopy=range(0,len(stats))
for i in loopy:
print(stats[i]['properties']['NAME'],stats[i]['properties']['preci_max'])
``````

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
from shapely.geometry import box, Polygon
from rasterstats import zonal_stats

def max_round(x):
m = x.max()
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,

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:
shapes,
crop = False,
nodata = -999 )
out_meta = src.meta

xminR = out_transform
xmaxR = xminR + out_meta['width']*out_transform
ymaxR = out_transform
yminR = ymaxR + out_meta['height']*out_transform

for i in range(len(geom)):
poly = Polygon(geom[i])
bounds = poly.bounds
xmin, ymin, xmax, ymax = bounds
rowR1 = int((ymaxR - ymax)/(-out_transform))
colR1 = int((xmin - xminR)/out_transform)
rowR2 = int((ymaxR - ymin)/(-out_transform))
colR2 = int((xmax - xminR)/out_transform)
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[i][j])
indices.append([i,j])

print(ma.max())
idx = ma.argmax()
print(indices[idx])
print(xminR + indices[idx]*out_transform + out_transform/2, ymaxR + indices[idx]*out_transform + out_transform/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.

• Thank you very much. I have inserted the masking operation into a loop on the polygons so I was able to find the maximum in each polygon. Otherwise not all the maximum values were correctly positioned. Thanks a lot again. – a.fornasiero Mar 20 at 12:08