import geopandas as gpd
geodf = gpd.GeoDataFrame.from_file('../tests/data/polygons.shp')
geodf.plot()

geodf

import rasterio
with rasterio.open("../tests/data/slope.tif") as src:
    transform = src.meta['transform']
    array = src.read_band(1)

print (transform)
array

array_chop = array[5:-5,5:-5]  # chop off the outer 5 pixels

import matplotlib.pyplot as plt
plt.imshow(array_chop, interpolation='nearest')
plt.colorbar()
plt.title('Slope Map')
plt.show()

from rasterstats import zonal_stats
zonal_stats(geodf, array)

transform

stats = zonal_stats(geodf, array, transform=transform)
stats

import pandas as pd

new_geodf = geodf.join(pd.DataFrame(stats))
new_geodf

## Hack around it for now
new_geodf.__class__ = gpd.GeoDataFrame
new_geodf.crs = {}
new_geodf.set_geometry('geometry')
## /hack

! rm /tmp/polygons_with_slope.geojson
new_geodf.to_file('/tmp/polygons_with_slope.geojson', driver="GeoJSON")

# Confirm that everything works according to OGR
! ogrinfo -ro -al /tmp/polygons_with_slope.geojson