diff --git a/supermercado/__init__.py b/supermercado/__init__.py
index c30af33..316a4ec 100644
--- a/supermercado/__init__.py
+++ b/supermercado/__init__.py
@@ -2,3 +2,4 @@
 from . import edge_finder as edgetiles
 from . import burntiles
 from . import uniontiles
+from . import tile_density
diff --git a/supermercado/scripts/cli.py b/supermercado/scripts/cli.py
index 9f6c12c..4e452d4 100644
--- a/supermercado/scripts/cli.py
+++ b/supermercado/scripts/cli.py
@@ -1,6 +1,6 @@
 import click, json
 import cligj
-from supermercado import edge_finder, uniontiles, burntiles, super_utils
+from supermercado import edge_finder, uniontiles, burntiles, super_utils, tile_density
 
 
 @click.group("supermercado")
@@ -49,6 +49,26 @@ def union(inputtiles, parsenames):
 cli.add_command(union)
 
 
+@click.command("heatmap")
+@click.argument("inputtiles", default="-", required=False)
+@click.option("--parsenames", is_flag=True)
+def heatmap(inputtiles, parsenames):
+    """
+    Returns the unioned shape of a stream of [<x>, <y>, <z>] tiles in GeoJSON.
+    """
+    try:
+        inputtiles = click.open_file(inputtiles).readlines()
+    except IOError:
+        inputtiles = [inputtiles]
+    density_feautures = tile_density.heatmap(inputtiles, parsenames)
+
+    for u in density_feautures:
+        click.echo(json.dumps(u))
+
+
+cli.add_command(heatmap)
+
+
 @click.command("burn")
 @cligj.features_in_arg
 @cligj.sequence_opt
diff --git a/supermercado/tile_density.py b/supermercado/tile_density.py
new file mode 100644
index 0000000..869b86b
--- /dev/null
+++ b/supermercado/tile_density.py
@@ -0,0 +1,79 @@
+import numpy as np
+
+from supermercado import super_utils as sutils
+from rasterio import features, Affine
+import mercantile
+
+
+def create_density(tiles, xmin, xmax, ymin, ymax, pad=1):
+    """Given an ndarray or tiles and range,
+    create a density raster of tile frequency
+    """
+    burn = np.zeros(
+        (xmax - xmin + (pad * 2 + 1), ymax - ymin + (pad * 2 + 1)), dtype=np.uint32
+    )
+
+    tiles, counts = np.unique(tiles, return_counts=True, axis=0)
+
+    burn[(tiles[:, 0] - xmin + pad, tiles[:, 1] - ymin + pad)] = counts
+
+    return burn
+
+
+def heatmap(inputtiles, parsenames):
+    """Creates a vector "heatmap" of tile densities
+
+    Parameters
+    ----------
+    inputtiles: sequence of [x, y, z] or strings
+        tiles to create a density heatmap from
+    parsenames: boolean
+        parse tile names from strings
+
+    Returns
+    -------
+    density_features: generator
+        sequence of features
+    """
+    tiles = sutils.tile_parser(inputtiles, parsenames)
+
+    xmin, xmax, ymin, ymax = sutils.get_range(tiles)
+
+    zoom = sutils.get_zoom(tiles)
+
+    # make an array of shape (xrange + 3, yrange + 3)
+    burn = sutils.burnXYZs(tiles, xmin, xmax, ymin, ymax, 0)
+
+    nw = mercantile.xy(*mercantile.ul(xmin, ymin, zoom))
+
+    se = mercantile.xy(*mercantile.ul(xmax + 1, ymax + 1, zoom))
+
+    aff = Affine(
+        ((se[0] - nw[0]) / float(xmax - xmin + 1)),
+        0.0,
+        nw[0],
+        0.0,
+        -((nw[1] - se[1]) / float(ymax - ymin + 1)),
+        nw[1],
+    )
+
+    unprojecter = sutils.Unprojecter()
+
+    density = create_density(tiles, xmin, xmax, ymin, ymax)
+
+    density = np.asarray(np.flipud(np.rot90(density)).astype(np.uint8), order="C")
+
+    density_features = (
+        {
+            "geometry": unprojecter.unproject(feature),
+            "properties": {"n": count},
+            "type": "Feature",
+        }
+        for feature, count in features.shapes(
+            density,
+            mask=density > 0,
+            transform=aff,
+        )
+    )
+
+    return density_features