Vectorize shapes (instead of pixel)

[vincentsarago]

Just wanted to see if it was possible. Answer: Not really but at least it was fun.

from rasterio._features import _shapes
from rasterio.rio.helpers import coords
from rasterio.transform import IDENTITY


cpdef rgb_to_hex(rgba):
    return '#{:02x}{:02x}{:02x}'.format(*rgba)


cpdef bytes shapes_encoder(
    data,
    mask,
    str layer_name = "layer",
    colormap = {},
):
    cdef float x, y
    cdef tuple rgba

    mvt = Tile()
    mvt_layer = Layer(mvt, layer_name.encode())

    for (p, v) in _shapes(data, mask, connectivity=4, transform=IDENTITY):
        feature = Polygon(mvt_layer)
        xs, ys = zip(*coords(p))
        feature.add_ring(len(xs))
        for (x, y) in zip(xs, ys):
            feature.set_point(int(x), int(y))

        if colormap:
            rgba = colormap.get(v, (0, 0, 0, 0))
            feature.add_property("rgba".encode(), rgb_to_hex(rgba).encode())
        else:
            feature.add_property("value".encode(), str(v).encode())
        feature.commit()

    return mvt.serialize()

problem:

• polygon returned by rasterio.features._shapes are quite complex and do need better handling (interior ring, multi polygon)
• Polygon are interesting the tile bounds and should not be polygon but LineString 🤷‍♂️

[jqtrde]

[vincentsarago]

update:

from rasterio._features import _shapes
from rasterio.rio.helpers import coords
from rasterio.transform import IDENTITY
from shapely import geometry

cpdef rgb_to_hex(rgba):
    return '#{:02x}{:02x}{:02x}'.format(*rgba)


cpdef bytes shapes_encoder(
    data,
    mask,
    str layer_name = "layer",
    colormap = {},
):
    cdef float x, y
    cdef tuple rgba

    mvt = Tile()
    mvt_layer = Layer(mvt, layer_name.encode())

    cdef int sc = 4096 // data.shape[1]

    for (p, v) in _shapes(data, mask, connectivity=4, transform=IDENTITY.scale(sc)):
        feature = Polygon(mvt_layer)
        feature.set_id(v)

        polygon = geometry.shape(p)

        feature.add_ring(len(polygon.exterior.coords))
        for x,y in polygon.exterior.coords:
            feature.set_point(int(x), int(y))

        for polygon_interior in list(polygon.interiors):
            feature.add_ring(len(polygon_interior.coords))
            for x,y in polygon_interior.coords:
                feature.set_point(int(x), int(y))

        if colormap:
            rgba = colormap.get(v, (0, 0, 0, 0))
            feature.add_property("color".encode(), rgb_to_hex(rgba).encode())

        feature.add_property("value".encode(), str(v).encode())
        feature.commit()

    return mvt.serialize()

By using shapely we can split interior/exterior rings. This is still pretty fast but maybe using GEOS directly will also be faster 🤷‍♂️

To Do:

• fix polygon crossing tile boundary
• Rasterio/GDAL produce self intersecting polygon which then don't get rendered properly
  

some reading on polygon validation:

• Polygon validity and performance tilezen/mapbox-vector-tile#100
• https://github.com/lycantropos/wagyu

[vincentsarago]

merge #2 in master but I'm getting shapely warnings Self-intersection at or near point! I guess I should really try to use wagyu as mentioned in #2 (comment)

[vincentsarago]

maybe we should just use pygeos instead of shapely 🤷‍♂️

[vincentsarago]

pygeos 0.10 has some nice improvement and may be worth to try https://github.com/pygeos/pygeos/releases/tag/0.10