How could I perform a Tabulate Area using GDAL and Numpy? I am trying to summarize how many pixels of each class in a Landcover raster are within each polygon of a Census Blocks shapefile.

I've tried playing with numpy.histogram2d, but didn't have any luck.

  • 2
    Have a look at the code of my QGIS plugin LecoS (or use it directly). It can extract raster values per polygon using scipy and PIL.
    – Curlew
    Oct 7, 2013 at 22:17

1 Answer 1


try something like this:

import gdal
import numpy as np

area_per_pixel = 100 #example...you change to suit

r = gdal.Open(rasterfile)
band = 1
raster_arr = np.array(r.GetRasterBand(band).ReadAsArray())

for cover in np.unique(raster_arr):
    tot_num_pixels = np.sum(raster_arr == cover)
    area = tot_num_pixels * area_per_pixel
    print cover, area

as for you shapefile and exactly what's going on there, I'm not sure since we cannot see it.

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.