TransWikia.com

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

Geographic Information Systems Asked by a.fornasiero on December 6, 2020

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'])

2 Answers

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.

Answered by xunilk on December 6, 2020

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))

Answered by YuriTheFury on December 6, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP