2

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()
    if m is np.ma.masked:
        return None
    else:
        return m

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

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


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

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:

enter image description here

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.

| improve this answer | |
  • 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

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.