docs: def stats(self, shp, stats='mean', add_stats=None, raster_out=True, *args, **kwargs):
          '''
          Compute raster statistics for a given geometry in shape, where shape is either
          a GeoPandas DataFrame, shapefile, or some other geometry format used by
          python-raster-stats. Runs python-raster-stats in background
          (additional help and info can be found there)
  
          Returns dataframe with statistics and clipped raster
  
          Usage:
  
          df = geo.stats(shape, stats=stats, add_stats=add_stats)
  
          where:
          raster_out: If True (Default), returns clipped Georasters
          '''
          df = pd.DataFrame(zonal_stats(shp, self.raster, nodata=self.nodata_value,
                                        all_touched=True, raster_out=raster_out,
                                        affine=Affine.from_gdal(*self.geot),
                                        geojson_out=True, stats=stats, add_stats=add_stats))
          df['GeoRaster'] = df.properties.apply(lambda x: GeoRaster(x['mini_raster_array'],
                                                                    Affine.to_gdal(x['mini_raster_affine']),
                                                                    nodata_value=x['mini_raster_nodata'],
                                                                    projection=self.projection,
                                                                    datatype=self.datatype))
          statcols = list(set([i for i in df.properties[0].keys()]).difference(set(shp.columns)))
          cols = shp.columns.tolist()+statcols
          cols = [i for i in cols if i != 'geometry' and i.find('mini_raster') == -1]
          df2 = pd.DataFrame([df.properties.apply(lambda x: x[i]) for i in cols]).T
          df2.columns = cols
          df2 = df2.merge(df[['id', 'GeoRaster']], left_index=True, right_index=True)
          df2.set_index('id', inplace=True)
          return df2