add heatmap function

README.md

@@ -58,9 +58,10 @@ Options:
   --help  Show this message and exit.
 
 Commands:
-  burn   Burn a stream of GeoJSONs into a output stream of the tiles they intersect for a given zoom.
-  edges  For a stream of [<x>, <y>, <z>] tiles, return only those tiles that are on the edge.
-  union  Returns the unioned shape of a stream of [<x>, <y>, <z>] tiles in GeoJSON.
+  burn     Burn a stream of GeoJSONs into a output stream of the tiles they intersect for a given zoom.
+  edges    For a stream of [<x>, <y>, <z>] tiles, return only those tiles that are on the edge.
+  heatmap  Creates a vector `heatmap` of tile densities.
+  union    Returns the unioned shape of a stream of [<x>, <y>, <z>] tiles in GeoJSON.
 ```
 
 ### `supermorecado burn`
@@ -122,6 +123,35 @@ cat tests/fixtures/france.geojson | supermorecado burn 6 --identifier WGS1984Qua
 
 ![](https://user-images.githubusercontent.com/10407788/236115946-2dfabe05-e51d-473b-9957-26bbbe9e61bb.jpg)
 
+
+### `supermorecado heatmap`
+
+```
+<[x, y, z] stream> | supermorecado heatmap --identifier {tms Identifier} | <{geojson} stream>
+```
+Outputs a stream of heatmap GeoJSON from an input stream of `[x, y, z]`s.
+
+Using default TMS (`WebMercatorQuad`)
+```
+cat tests/fixtures/heatmap.txt| supermorecado heatmap | fio collect | geojsonio
+```
+
+![](https://github.com/developmentseed/supermorecado/assets/10407788/fec78304-12e8-47de-bca3-af58f04e8824)
+
+Using other TMS (e.g `WGS1984Quad`)
+
+```
+# create a list of tiles
+cat tests/fixtures/france.geojson | supermorecado burn 6 --identifier WGS1984Quad > france_wgs84_z6.txt
+# randomly append more tiles
+for run in {1..10}; do cat france_wgs84_z6.txt | sort -R | head -n 2 >> france_wgs84_z6.txt; done
+
+cat france_wgs84_z6.txt |  supermorecado heatmap --identifier WGS1984Quad | fio collect | geojsonio
+```
+
+![](https://github.com/developmentseed/supermorecado/assets/10407788/dec3dbdb-6030-4dd8-9efc-5c8d4b7c83d3)
+
+
 ## API migration
 
 `supermorecado` is really similar to `supermercado` (it reuse most of the code) but with the addition of multiple TMS support from morecantile.

supermorecado/__init__.py

@@ -6,5 +6,6 @@
 __version__ = "0.1.0"
 
 from .burntiles import burnTiles  # noqa
+from .densitytiles import densityTiles  # noqa
 from .edge_finder import findedges  # noqa
 from .uniontiles import unionTiles  # noqa

supermorecado/densitytiles.py

@@ -0,0 +1,74 @@
+"""Tiles heatmap but for other TMS.
+
+adapted from https://github.com/mapbox/supermercado/pull/53
+"""
+
+from typing import Dict, List
+
+import attr
+import morecantile
+import numpy
+import numpy.typing as npt
+from affine import Affine
+from rasterio import features
+
+from supermorecado import super_utils as sutils
+
+WEB_MERCATOR_TMS = morecantile.tms.get("WebMercatorQuad")
+
+
+def create_density(
+    tiles: npt.NDArray, xmin: int, xmax: int, ymin: int, ymax: int, pad: int = 1
+) -> numpy.ndarray:
+    """Given an ndarray or tiles and range, create a density raster of tile frequency."""
+    burn = numpy.zeros(
+        (xmax - xmin + (pad * 2 + 1), ymax - ymin + (pad * 2 + 1)), dtype=numpy.uint32
+    )
+    tiles, counts = numpy.unique(tiles, return_counts=True, axis=0)
+    burn[(tiles[:, 0] - xmin + pad, tiles[:, 1] - ymin + pad)] = counts
+    return burn
+
+
+@attr.s
+class densityTiles:
+    """heatmap."""
+
+    tms: morecantile.TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
+
+    def heatmap(self, tiles: npt.NDArray) -> List[Dict]:
+        """Creates a vector "heatmap" of tile densities"""
+        xmin, xmax, ymin, ymax = sutils.get_range(tiles)
+
+        zoom = sutils.get_zoom(tiles)
+
+        burn = create_density(tiles, xmin, xmax, ymin, ymax, 0)
+        density = numpy.asarray(
+            numpy.flipud(numpy.rot90(burn)).astype(numpy.uint8),
+            order="C",
+        )
+
+        nw = self.tms.xy(*self.tms.ul(xmin, ymin, zoom))
+        se = self.tms.xy(*self.tms.ul(xmax + 1, ymax + 1, zoom))
+        afftrans = 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(self.tms)
+
+        return [
+            {
+                "geometry": unprojecter.unproject(feature),
+                "properties": {"n": count},
+                "type": "Feature",
+            }
+            for feature, count in features.shapes(
+                density,
+                mask=density > 0,
+                transform=afftrans,
+            )
+        ]

supermorecado/edge_finder.py

@@ -27,11 +27,12 @@
 """
 
 import numpy
+import numpy.typing as npt
 
 from supermorecado import super_utils as sutils
 
 
-def findedges(tiles: numpy.array):
+def findedges(tiles: npt.NDArray):
     """Find edges."""
     xmin, xmax, ymin, ymax = sutils.get_range(tiles)
 

supermorecado/scripts/cli.py

@@ -6,7 +6,7 @@
 import cligj
 import morecantile
 
-from supermorecado import burnTiles, findedges
+from supermorecado import burnTiles, densityTiles, findedges
 from supermorecado import super_utils as sutils
 from supermorecado import unionTiles
 
@@ -105,3 +105,38 @@ def burn(features, sequence, zoom, identifier, tms):
     burntiles = burnTiles(tms=tilematrixset)
     for t in burntiles.burn(features, zoom):
         click.echo(t.tolist())
+
+
+@cli.command(short_help="Creates a vector `heatmap` of tile densities.")
+@click.argument("inputtiles", default="-", required=False)
+@click.option("--parsenames", is_flag=True)
+@click.option(
+    "--identifier",
+    type=click.Choice(morecantile.tms.list()),
+    default="WebMercatorQuad",
+    help="TileMatrixSet identifier.",
+)
+@click.option(
+    "--tms",
+    help="Path to TileMatrixSet JSON file.",
+    type=click.Path(),
+)
+def heatmap(inputtiles, parsenames, identifier, tms):
+    """
+    Creates a vector `heatmap` of tile densities.
+    """
+    try:
+        inputtiles = click.open_file(inputtiles).readlines()
+    except IOError:
+        inputtiles = [inputtiles]
+
+    tiles = sutils.tile_parser(inputtiles, parsenames)
+
+    tilematrixset = morecantile.tms.get(identifier)
+    if tms:
+        with open(tms, "r") as f:
+            tilematrixset = morecantile.TileMatrixSet(**json.load(f))
+
+    tiledensity = densityTiles(tms=tilematrixset)
+    for u in tiledensity.heatmap(tiles):
+        click.echo(json.dumps(u))

supermorecado/super_utils.py

@@ -27,15 +27,15 @@
 """
 import json
 import re
-from typing import Any, Dict, Generator, Sequence, Tuple
+from typing import Any, Dict, Generator, List, Sequence, Tuple
 
 import attr
 import morecantile
 import numpy
 import numpy.typing as npt
 
 
-def parseString(tilestring, matcher):
+def parseString(tilestring, matcher) -> List[int]:
     """Parse Tile."""
     tile = [int(r) for r in matcher.match(tilestring).group().split("-")]
     tile.append(tile.pop(0))

supermorecado/uniontiles.py

@@ -26,6 +26,8 @@
 
 """
 
+from typing import Dict, List
+
 import attr
 import morecantile
 import numpy
@@ -44,13 +46,12 @@ class unionTiles:
 
     tms: morecantile.TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
 
-    def union(self, tiles: npt.NDArray):
+    def union(self, tiles: npt.NDArray) -> List[Dict]:
         """Union of tiles."""
         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 = self.tms.xy(*self.tms.ul(xmin, ymin, zoom))

tests/expected/heatmap.txt

@@ -0,0 +1,8 @@
+{"geometry": {"type": "Polygon", "coordinates": [[[-5.625000000000113, 52.482780222078254], [-5.625000000000113, 48.9224992637583], [-1.2131008896728105e-13, 48.9224992637583], [-1.2131008896728105e-13, 52.482780222078254], [-5.625000000000113, 52.482780222078254]]]}, "properties": {"n": 6.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[-1.2131008896728105e-13, 52.482780222078254], [-1.2131008896728105e-13, 48.9224992637583], [5.62499999999987, 48.9224992637583], [5.62499999999987, 52.482780222078254], [-1.2131008896728105e-13, 52.482780222078254]]]}, "properties": {"n": 4.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[5.62499999999987, 52.482780222078254], [5.62499999999987, 48.9224992637583], [11.24999999999986, 48.9224992637583], [11.24999999999986, 52.482780222078254], [5.62499999999987, 52.482780222078254]]]}, "properties": {"n": 6.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[-5.625000000000113, 48.9224992637583], [-5.625000000000113, 45.08903556483109], [-1.2131008896728105e-13, 45.08903556483109], [-1.2131008896728105e-13, 48.9224992637583], [-5.625000000000113, 48.9224992637583]]]}, "properties": {"n": 2.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[-1.2131008896728105e-13, 48.9224992637583], [-1.2131008896728105e-13, 45.08903556483109], [11.24999999999986, 45.08903556483109], [11.24999999999986, 48.9224992637583], [-1.2131008896728105e-13, 48.9224992637583]]]}, "properties": {"n": 1.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[-5.625000000000113, 45.08903556483109], [-5.625000000000113, 40.97989806962022], [-1.2131008896728105e-13, 40.97989806962022], [-1.2131008896728105e-13, 45.08903556483109], [-5.625000000000113, 45.08903556483109]]]}, "properties": {"n": 4.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[-1.2131008896728105e-13, 45.08903556483109], [-1.2131008896728105e-13, 40.97989806962022], [5.62499999999987, 40.97989806962022], [5.62499999999987, 45.08903556483109], [-1.2131008896728105e-13, 45.08903556483109]]]}, "properties": {"n": 3.0}, "type": "Feature"}
+{"geometry": {"type": "Polygon", "coordinates": [[[5.62499999999987, 45.08903556483109], [5.62499999999987, 40.97989806962022], [11.24999999999986, 40.97989806962022], [11.24999999999986, 45.08903556483109], [5.62499999999987, 45.08903556483109]]]}, "properties": {"n": 2.0}, "type": "Feature"}

tests/fixtures/heatmap.txt

@@ -0,0 +1,29 @@
+[31, 21, 6]
+[32, 21, 6]
+[33, 21, 6]
+[31, 22, 6]
+[32, 22, 6]
+[33, 22, 6]
+[31, 23, 6]
+[32, 23, 6]
+[33, 23, 6]
+[31, 22, 6]
+[31, 21, 6]
+[31, 21, 6]
+[31, 21, 6]
+[33, 21, 6]
+[31, 23, 6]
+[31, 21, 6]
+[31, 21, 6]
+[32, 21, 6]
+[32, 23, 6]
+[33, 21, 6]
+[33, 21, 6]
+[31, 23, 6]
+[31, 23, 6]
+[33, 21, 6]
+[33, 21, 6]
+[32, 21, 6]
+[32, 21, 6]
+[33, 23, 6]
+[32, 23, 6]

tests/test_cli.py

@@ -74,6 +74,22 @@ def test_burn_cli():
     assert result.output == expected
 
 
+def test_heatmap_cli():
+    filename = os.path.join(os.path.dirname(__file__), "fixtures/heatmap.txt")
+    expectedFilename = os.path.join(os.path.dirname(__file__), "expected/heatmap.txt")
+
+    runner = CliRunner()
+    result = runner.invoke(cli, ["heatmap", filename])
+    assert result.exit_code == 0
+
+    with open(expectedFilename) as ofile:
+        expected = ofile.readlines()
+        expected_counts = [json.loads(f)["properties"]["n"] for f in expected]
+
+    count = [json.loads(f)["properties"]["n"] for f in result.output.splitlines()]
+    assert count == expected_counts
+
+
 def test_burn_tile_center_point_roundtrip():
     tile = [83885, 202615, 19]
     w, s, e, n = mercantile.bounds(*tile)