rasterstats zonal_stats - return cell position

Question

I've been using the rasterstats zonal_stats method to return the min pixel value (elevation) within a buffered polygon. To return the pixel's position I thought I could used the raster_out=True option. However, I've noticed the coordinate that is reflected in the affine is actually the extent of the pixels that the polygon touches (since I am using all_touch=True).

I am looking for a way to return any sort of pixel positional data (bounding box, centroid, xy, lat-lon, etc.) for a given stat method (minimum in my case). I'm currently thinking that rasterio might be a good canidate to help solve this problem, but I'm open to other options.

Suggestions?

from rasterstats import zonal_stats
import geopandas as gpd
import rasterio as rio
from shapely.geometry import Point

# Load point and raster
gpd_point = gpd.read_file("./gpd_points.shp")
with rio.open("./dem.tif") as src:
    affine = src.transform
    arr = src.read(1)

# Buffer point to 3 feet
buff = gpd_point.buffer(3)
# buff.to_file('./test_buff.shp', crs=src.crs, driver='ESRI Shapefile')

# Get min stats from all raster cells that touch buffer
stats = zonal_stats(buff, arr, 
                    stats="max", 
                    affine=affine, 
                    all_touched=True,
                    raster_out=True)

# Identify min and pixel coords
min_elevation_pixel = stats[0]['mini_raster_affine']
x_upper_left = min_elevation_pixel[2]
y_upper_left = min_elevation_pixel[5]

# Output
# gdf = gpd.GeoDataFrame(geometry=[Point(x_upper_left, y_upper_left)])
# gdf.to_file(f'./test_min_point.shp', crs=src.crs, driver='ESRI Shapefile')

# print(stats)
# print(src.crs)
# print(gpd_point.crs)
# print(buff.crs)

Answer 1

You can grab the indices of the mini raster where the pixel value == the statistic value (i.e min or max) using numpy.argwhere, then use the transform/affine to work out the map coordinates.

e.g.

from rasterstats import zonal_stats
import numpy as np
from rasterstats import zonal_stats
import geopandas as gpd
import rasterio as rio
from shapely.geometry import Point

stat = "max"

# Load point and raster
gpd_point = gpd.read_file("./point.shp")
with rio.open("./raster.tif") as src:
    affine = src.transform
    arr = src.read(1)

# Buffer point
buff = gpd_point.buffer(3)
buff.to_file('./buff.shp', crs=src.crs, driver='ESRI Shapefile')

# Get min stats from all raster cells that touch buffer
stats = zonal_stats(buff, arr,
                    stats=stat,
                    affine=affine,
                    all_touched=True,
                    raster_out=True)

# get indices of all pixels that match the stat value
indices = np.argwhere(stats[0]['mini_raster_array']==stats[0][stat])

# just the first match for demo
py, px = indices[0]

# calculate coordinates  
# for top left of pixel, use:
# mx, my = stats[0]['mini_raster_affine'] * (px, py)

# for centroid of pixel, use:
mx, my = stats[0]['mini_raster_affine'] * (px+0.5, py+0.5)

gdf = gpd.GeoDataFrame(geometry=[Point(mx, my)])
gdf.to_file(f'./{stat}_point.shp', crs=src.crs, driver='ESRI Shapefile')