Add support to specify the bands to calculate statistics on to the python API

[theroggy]

Now all bands are calculated, or you need to choose one band when providing a RasterSource.

[dbaston]

Something like exact_extract(rast, poly, 'mean', bands=[1,2,3]) ?

It's clear what this would mean if rast is a filename. Less clear what it means if rast is a list of filenames, or a list of RasterSources.

[theroggy]

Ideally it would be possible to combine it with multiple images as input... So par example something like this:

rast=[("image1.tif", [1, 2, 3]), ("image2.tif", [4, 5, 6])]

That's also the signature I'm now using for a helper function to calculate zonal statistics in bulk: zonal_stats_bulk