diff --git a/.travis.yml b/.travis.yml
index 6260062..93543b4 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -28,3 +28,4 @@ install:
 script: py.test --cov rasterstats --cov-report term-missing
 after_success:
   - coveralls
+
diff --git a/src/rasterstats/cli.py b/src/rasterstats/cli.py
index d365439..81cb9b1 100644
--- a/src/rasterstats/cli.py
+++ b/src/rasterstats/cli.py
@@ -24,10 +24,12 @@
 @click.option('--nodata', type=int, default=None)
 @click.option('--prefix', type=str, default='_')
 @click.option('--stats', type=str, default=None)
+@click.option('--keep_ids', default=True)
+@click.option('--keep_properties', default=False)
 @cligj.sequence_opt
 @cligj.use_rs_opt
 def zonalstats(features, raster, all_touched, band, categorical,
-               indent, info, nodata, prefix, stats, sequence, use_rs):
+               indent, info, nodata, prefix, stats,keep_ids, keep_properties, sequence, use_rs):
     '''zonalstats generates summary statistics of geospatial raster datasets
     based on vector features.
 
@@ -61,6 +63,8 @@ def zonalstats(features, raster, all_touched, band, categorical,
         nodata=nodata,
         stats=stats,
         prefix=prefix,
+        preserve_properties=keep_properties,
+        preserve_ids=keep_ids,
         geojson_out=True)
 
     if sequence:
diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py
index 540ed20..e656fb8 100644
--- a/src/rasterstats/main.py
+++ b/src/rasterstats/main.py
@@ -42,7 +42,9 @@ def gen_zonal_stats(
         zone_func=None,
         raster_out=False,
         prefix=None,
-        geojson_out=False, **kwargs):
+        geojson_out=False,
+        preserve_properties=False,
+        preserve_ids=False, **kwargs):
     """Zonal statistics of raster values aggregated to vector geometries.
 
     Parameters
@@ -108,6 +110,12 @@ def gen_zonal_stats(
         Original feature geometry and properties will be retained
         with zonal stats appended as additional properties.
         Use with `prefix` to ensure unique and meaningful property names.
+    
+    preserve_properties: boolean (default: False)
+        preserve the properties of each feature in the returned stats data
+    
+    preserve_ids: boolean (default: False)
+        Preserve the IDs of each feature in the returned stats data
 
     Returns
     -------
@@ -195,6 +203,13 @@ def gen_zonal_stats(
                         feature_stats = remap_categories(category_map, feature_stats)
                 else:
                     feature_stats = {}
+                
+            
+                if preserve_properties and 'properties' in feat:
+                    feature_stats['properties'] = feat['properties']
+                    if preserve_ids:
+                        if 'id' in feat['properties']:
+                            feature_stats['id'] = feat['properties']['id']
 
                 if 'min' in stats:
                     feature_stats['min'] = float(masked.min())
diff --git a/tests/test_zonal.py b/tests/test_zonal.py
index 48babc9..6cd5f64 100644
--- a/tests/test_zonal.py
+++ b/tests/test_zonal.py
@@ -29,6 +29,22 @@ def test_main():
     assert round(stats[0]['mean'], 2) == 14.66
 
 
+def test_polygon_id():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    features = read_features(polygons)
+
+    stats = zonal_stats(polygons, raster, preserve_ids=True)
+    for x, feature in enumerate(features):
+        assert stats[x]['id'] == feature['properties']['id']
+
+def test_polygon_properties():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    features = read_features(polygons)
+    stats = zonal_stats(polygons, raster, preserve_properties=True)
+    for x, feature in enumerate(features):
+        assert stats[x]['properties'] == feature['properties']
+
+
 # remove after band_num alias is removed
 def test_band_alias():
     polygons = os.path.join(DATA, 'polygons.shp')