Tabular el área usando GDAL y Numpy

Question

¿Cómo podría realizar una Tabulate Area usando GDAL y Numpy? Estoy tratando de resumir cuántos píxeles de cada clase en un raster de Landcover están dentro de cada polígono de un shapefile de Census Blocks.

He intentado jugar con numpy.histogram2d pero no ha tenido suerte.

Answer 1

prueba algo como esto:

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

En cuanto a tu shapefile y a lo que ocurre exactamente allí, no estoy seguro ya que no podemos verlo.