diff --git a/README.md b/README.md index ee22978..366278b 100644 --- a/README.md +++ b/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 [, , ] tiles, return only those tiles that are on the edge. - union Returns the unioned shape of a stream of [, , ] 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 [, , ] 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 [, , ] 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. diff --git a/supermorecado/__init__.py b/supermorecado/__init__.py index 9b9b006..dc9fe52 100644 --- a/supermorecado/__init__.py +++ b/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 diff --git a/supermorecado/densitytiles.py b/supermorecado/densitytiles.py new file mode 100644 index 0000000..58e33ae --- /dev/null +++ b/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, + ) + ] diff --git a/supermorecado/edge_finder.py b/supermorecado/edge_finder.py index ee2645d..001f49b 100644 --- a/supermorecado/edge_finder.py +++ b/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) diff --git a/supermorecado/scripts/cli.py b/supermorecado/scripts/cli.py index 707ba12..b843e5e 100644 --- a/supermorecado/scripts/cli.py +++ b/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)) diff --git a/supermorecado/super_utils.py b/supermorecado/super_utils.py index 24e29a1..089cdb5 100644 --- a/supermorecado/super_utils.py +++ b/supermorecado/super_utils.py @@ -27,7 +27,7 @@ """ 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 @@ -35,7 +35,7 @@ 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)) diff --git a/supermorecado/uniontiles.py b/supermorecado/uniontiles.py index 3a49426..a50d2e7 100644 --- a/supermorecado/uniontiles.py +++ b/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)) diff --git a/tests/expected/heatmap.txt b/tests/expected/heatmap.txt new file mode 100644 index 0000000..9c0c428 --- /dev/null +++ b/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"} diff --git a/tests/fixtures/heatmap.txt b/tests/fixtures/heatmap.txt new file mode 100644 index 0000000..b903abd --- /dev/null +++ b/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] diff --git a/tests/test_cli.py b/tests/test_cli.py index f35df85..df6cba8 100644 --- a/tests/test_cli.py +++ b/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)