ArcPy – How to Access Raster Values to Create Histogram

arcpyhistogramraster-dataset

Using native Python and arcpy function what method is best to create a histogram of raster dataset? I know it's possible to use raterio and GDAL functions, but I would like to keep it a arcpy and native Python.

Specifically I am working on a script to delineate watersheds and I would like to create a raster based on the values in a flow accumulation raster dataset. From the calculated raster I would then convert to vector stream layers. The ultimate goal is to be able to create separate stream order layers.

Best Answer

I would recommend converting your raster to a numpy array, which you can then use to calculate statistics for the histogram. matplotlib is then used for the plotting. For example:

import arcpy
import numpy as np
import matplotlib.pyplot as plt

# This example uses an 8bit unsigned format TIFF image
tiff = r'C:\path\to\your\image.tif'

# Convert raster to numpy array and calculate histogram
array = arcpy.RasterToNumPyArray(tiff)
hist, bins = np.histogram(array, bins = [25,50,75,100,125,150,175,200,225,250,275])

# Plot the histogram
plt.bar(bins[:-1], hist, width = 25)
plt.xlim(min(bins), max(bins))
plt.ylabel('Pixel Count')
plt.xlabel('Pixel Values')
plt.show() 

enter image description here


However, I believe Rasterio does a better job plotting raster histograms. For example, here is a Rasterio generated histogram of 4 band NAIP imagery:

import rasterio
from rasterio.plot import show_hist

src = rasterio.open("path/to/your/image.tif")
show_hist(
    src, bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram")

enter image description here