# Writing raster statistics in Python as fast as ArcGIS

I am trying to run a Zonal Statistics in Python. My zone raster has about 300,000 different zones, and in each zone I want to calculate several statistics (mean, median, stdev, count, etc.).

I've used the rasterstats package, xrspatial package, and written the code myself in numpy by looping through each zone value. However, for loops in Python are very slow and create a bottleneck, so all of these approaches take about an hour to run. I used the "Zonal Statistics as Table" tool in ArcGIS Pro and it took < 10 seconds. How did they write this code so fast?

Is there a way to do this in Python without looping through each zone value, and utilizing numpy's vectorized functions?

Here is an example of my implementation in numpy:

``````import numpy as np
import rasterio
vals = np.unique(zoneRas)
mean = []
for val in vals:
results = ndviRas[zoneRas == val]
mean.append(results.mean())
return (mean)
``````
• Load zones and values in 2 columns of pandas dataframe and use group by method. Commented Apr 11, 2023 at 4:23
• Thanks, this is what I was looking for. Commented Apr 11, 2023 at 16:14

Below is my solution using only pandas and numpy based on the suggestion from @FelixIP.

The shape of the arrays are (5525, 6563) with 150,000 zones. Runtime is < 2 seconds to calculate 6 stats at once (min, max, mean, median, std, count):

``````def zonal_stats(data, zones, nodata_val):
import pandas as pd
import numpy as np

# First, group all missing data into one zone for easy handling
zones = np.where(data != nodata_val, zones, nodata_val)

# Create pandas dataframe
f = pd.DataFrame({'ID': zones, 'Value': data})

# Remove no data zones
df = df[df.ID != nodata_val]

# Calculate stats and add results to table as a column
df = df.groupby('ID').agg(
{'Value': ['min', 'max', 'mean', 'median', 'std', 'count']}
)

# Drop multilevel column and rename columns if necessary
df = df.droplevel(0, axis=1).rename(
columns = {'min': 'MinValue',
'max': 'MaxValue',
'mean': 'MeanValue',
'median': 'MedianValue',
'std': 'StdevValue',
'count': 'CountValue'})

# Removes ID column from index and makes a column again
df.reset_index(inplace = True)

return(df)

# Function takes a 1D array, so flatten the rasters first using ravel
zonal_stats(ndviRas.ravel(), zoneRas.ravel(), -100)
``````
• Glad it works but using transform populates all 5525*6553 rows of dataframe. This statology.org/pandas-groupby-multiple-aggregations might work much faster. Commented Apr 12, 2023 at 1:15
• You're right. Edited the post to use the groupby().agg() method which is more efficient than groupby().transform() for this purpose. Commented Apr 13, 2023 at 0:17

I'm not sure if it is mandatory for your case that the input zones are in raster format, but I recently did some tests with calculating zonal stats on a polygon layer vs. a raster file.

Up to now, the fastest solution I found was the QgsZonalStatistics function in pyqgis. The code for the benchmarks used can be found in this github repo, in the benchmarks_zonalstats dir.

The results for 15.000 polygons on a 10 meter pixel size raster: