raster encoder [do not merge]

rio_vectortiles/__init__.py

@@ -2,17 +2,16 @@
 import mercantile
 import rasterio
 
-from vtzero.tile import Tile, Layer, Polygon
 import gzip
 from io import BytesIO
-from affine import Affine
+
 import numpy as np
 from rasterio.warp import reproject
 from rasterio.io import MemoryFile
-from rasterio.features import shapes
 from rasterio.transform import from_bounds
 from rasterio.enums import Resampling
 import warnings
+from PIL import Image
 
 warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning)
 
@@ -70,37 +69,25 @@ def read_transform_tile(
 
     with rasterio.open(src_path) as src:
         src_band = rasterio.band(src, bidx=1)
-        vtile = Tile()
-
-        layer = Layer(vtile, layer_name.encode(), extent=extent)
 
         with MemoryFile() as mem:
             with mem.open(**dst_kwargs) as dst:
                 dst_band = rasterio.band(dst, bidx=1)
                 reproject(src_band, dst_band, resampling=Resampling.mode)
-                dst.transform = Affine.identity()
-                if interval is None and not len(filters):
-                    vectorizer = shapes(dst_band)
-                else:
-                    data = dst.read(1)
-                    for f in filters:
-                        data = f(data)
-                    data = (data // interval * interval).astype(np.int32)
-                    vectorizer = shapes(data)
-
-                for s, v in vectorizer:
-                    feature = Polygon(layer)
-                    for ring in s["coordinates"]:
-
-                        feature.add_ring(len(ring))
-                        for x, y in ring:
-                            feature.set_point(x, y)
-                        feature.close_ring()
-                    feature.add_property(b"v", f"{int(v)}".encode())
-                    feature.commit()
+                data = dst.read(1)
+                A = data / 256
+                B = data // 256
+                C = B / 256
+                D = B // 256
+
+                blue = ((A - B) * 256).astype(np.uint8)
+                green = ((C - D) * 256).astype(np.uint8)
+                red = (((D / 256) - (D // 256)) * 256).astype(np.uint8)
+
+                img = Image.fromarray(np.dstack([red, green, blue]))
 
     with BytesIO() as dst:
-        with gzip.open(dst, mode="wb") as gz:
-            gz.write(vtile.serialize())
+        img.save(dst, format="png")
+
         dst.seek(0)
         return dst.read(), tile

rio_vectortiles/scripts/cli.py

@@ -7,8 +7,18 @@
 import rasterio
 from rasterio.warp import transform_bounds
 from functools import partial
+
 from concurrent.futures import ThreadPoolExecutor
-from rio_vectortiles import read_transform_tile, decompress_tile
+from rio_vectortiles import read_transform_tile
+from pathlib import Path
+
+
+def make_tile_iterable(w, s, e, n, zoom, buffer=0):
+    UL = mercantile.tile(w, n, zoom)
+    LR = mercantile.tile(e, s, zoom)
+    for x in range(max(UL.x - buffer, 0), min(LR.x + buffer, 2**zoom)):
+        for y in range(max(UL.y - buffer, 0), min(LR.y + buffer, 2**zoom)):
+            yield mercantile.Tile(x, y, zoom)
 
 
 @click.group("vectortiles")
@@ -42,8 +52,14 @@ def dump_tiles(mbtiles, output_directory):
         "SELECT tile_column, tile_row, zoom_level, tile_data from tiles"
     ):
         tiley = int(2**z) - y - 1
-        with open(f"{output_directory}/{z}-{x}-{tiley}.mvt", "wb") as dst:
-            dst.write(decompress_tile(b))
+        dst_path = Path(f"{output_directory}/{z}/{x}/{tiley}.png")
+        for i, p in enumerate(dst_path.parts[1:-1], 1):
+            dst_folder = Path(*dst_path.parts[: i + 1])
+            if not os.path.isdir(dst_folder):
+                os.mkdir(dst_folder)
+
+        with open(dst_path, "wb") as dst:
+            dst.write(b)
 
 
 @cli.command("clump", short_help="Clump and index tiles from an mbtiles")
@@ -100,8 +116,11 @@ def vectortiles(
         maxzoom = get_maxzoom(*wm_bounds, *src.shape, max_extent) + zoom_adjust
 
         wgs_bounds = transform_bounds(src.crs, "EPSG:4326", *src.bounds)
-
-        tiles = list(mercantile.tiles(*wgs_bounds, range(maxzoom + 1)))
+        zoom_groups = (
+            make_tile_iterable(*wgs_bounds, z, 1) for z in range(maxzoom + 1)
+        )
+        tiles = [t for z in zoom_groups for t in z]
+        # tiles = list(mercantile.tiles(*wgs_bounds, range(maxzoom + 1)))
         click.echo(f"Generating {len(tiles):,} tiles from zooms 0-{maxzoom}", err=True)
         dst_profile = dict(
             driver="GTiff", count=1, dtype=src.meta["dtype"], crs="EPSG:3857"