Calculating raster variance within an area of interest in Python

Question

The zonal statistics tool within QGIS is able to calculate the average raster variance within an area of interest. I am now trying to automate this in Python.

I typically use zonal_stats from the rasterstats library for calculating zonal statistics in Python, but zonal_stats is only able to calculate min, max, mean, count, sum, std, median, majority, minority, unique, range, nodata, and percentile.

How would I calculate raster variance over an area of interest in Python?

Answer 1

I figured it out- instead of using rasterstats, you can use rasterio and numpy together to calculate raster variance. Below is the code I ended up using.

import numpy as np
import rasterio


def calculate_raster_variance(raster_path: str):
    raster = rasterio.open(raster_path)
    arr = raster.read()
    out = np.var(arr)
    return out