# code
from descartes import PolygonPatch
import geopandas as gp
import pysal as ps

from matplotlib import pyplot as plt
from matplotlib.patches import Polygon as mpl_Polygon
from matplotlib.collections import PatchCollection

shapefile = "filepolygon.shp"

df_map_elements = gp.GeoDataFrame.from_file(shapefile)


df_map_elements["mpl_polygon"] = np.nan
print type(df_map_elements['mpl_polygon'])

df_map_elements['mpl_polygon'] = df_map_elements['mpl_polygon'].astype(object)
for self_index, self_row_df in df_map_elements.iterrows():
    m_polygon = self_row_df['geometry']
    poly=[]
    if m_polygon.geom_type == 'MultiPolygon':
        for pol in m_polygon:
             poly.append(PolygonPatch(pol))
    else:
        poly.append(PolygonPatch(m_polygon))
    df_map_elements.set_value(self_index, 'mpl_polygon', poly)

sourceDataGdalObj = gdal.Open("rgb_landsat.tif")
ulx, xres, xskew, uly, yskew, yres = sourceDataGdalObj.GetGeoTransform()
lrx = ulx + (sourceDataGdalObj.RasterXSize * xres)
lry = uly + (sourceDataGdalObj.RasterYSize * yres)
RasterBBox = [ulx, lrx, uly, lry]

outputArray = []
for band in range(sourceDataGdalObj.RasterCount):

    dataGdalArray = sourceDataGdalObj.GetRasterBand(band+1).ReadAsArray()
    outputArray.append(dataGdalArray)

dict_mapindex_mpl_polygon = df_map_elements['mpl_polygon'].to_dict()
fig, ax = plt.subplots()

for c_l ,patches in dict_mapindex_mpl_polygon.items():
    p = PatchCollection(patches,color='none',lw=2.5,edgecolor='black')
    ax.imshow(np.dstack(outputArray), interpolation='nearest', extent=RasterBBox)
    ax.add_collection(p)

ax.autoscale_view()