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)