import gdal,ogr, osr

in_path = 'C:\Test\Multiband Image.tif'
out_path = 'C:\Test\'

ds = gdal.Open(in_path)
ds_proj=ds.GetProjection()

raster_proj = osr.SpatialReference()
raster_proj.ImportFromWkt(ds_proj)

for i in range (1,ds.RasterCount+1): ##Begin at band 1 through to to the last one
    #Save each band individually
    srcband = ds.GetRasterBand(i)
    out_ds = gdal.Translate(out_path + 'band' + str(i) + '.tiff', ds, format='GTiff', bandList=[i])
    out_ds=None

    #Prepare shapefile
    outShapefile = "C:\Test\polygonized"
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outDatasource = driver.CreateDataSource(outShapefile+ str(i) + ".shp")
    outLayer = outDatasource.CreateLayer("polygonized", srs=raster_proj)
    #Add the DN field
    newField = ogr.FieldDefn('DN', ogr.OFTInteger)
    outLayer.CreateField(newField)
    gdal.Polygonize(srcband, None, outLayer, 0, [], callback=None )
    outDatasource.Destroy()
    sourceRaster = None

print "Done"