1

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)

enter image description here

1 Answer 1

5

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

enter image description here

2
  • Perfect, I was just thinking about using the indices this way too. np.argwhere works nicely. Just so that I understand the documentation correctly np.argwhere returns [[row column]] position of the index?
    – Binx
    Commented Mar 22, 2023 at 14:35
  • 1
    Yes that's right, numpy is always row, col. Which is why I've used py, px = indices[0]
    – user2856
    Commented Mar 22, 2023 at 20:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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