Skip to content

MultiVectorLayer in single SQL query with PostGIS #43

@nitrag

Description

@nitrag

I wrote this up for our purposes but perhaps you can incorporate it in a cleaner way than I could:

@dataclass
class VectorLayerQueryset:
    layer_name: str
    queryset: QuerySet
    queryset_limit = None
    fields: tuple = None
    geom_name: str = 'geom'


class MultiLayerVectorTile:
    """
    Base Mixin to handle vector tile generation
    """

    vector_tile_content_type = "application/x-protobuf"
    vector_layers: list[VectorLayerQueryset] = None

    vector_tile_extent = 4096  # define tile extent
    vector_tile_buffer = (
        256  # define buffer around tiles (intersected polygon display without borders)
    )
    vector_tile_clip_geom = (
        True  # define if feature geometries should be clipped in tile
    )

    @classmethod
    def get_bounds(cls, x, y, z):
        return mercantile.xy_bounds(x, y, z)

    def get_tile(self, x, y, z):
        # get tile coordinates from x, y and z
        xmin, ymin, xmax, ymax = self.get_bounds(x, y, z)

        queries = []
        params = []
        for layer in self.vector_layers:
            queryset = layer.queryset

            # keep features intersecting tile
            filters = {
                # GeoFuncMixin implicitly transforms to SRID of geom
                f"{layer.geom_name}__intersects": MakeEnvelope(
                    xmin, ymin, xmax, ymax, 3857
                )
            }
            queryset = queryset.filter(**filters)
            # annotate prepared geometry for MVT
            queryset = queryset.annotate(
                geom_prepared=AsMVTGeom(
                    Transform(layer.geom_name, 3857),
                    MakeEnvelope(xmin, ymin, xmax, ymax, 3857),
                    self.vector_tile_extent,
                    self.vector_tile_buffer,
                    self.vector_tile_clip_geom,
                )
            )
            fields = (
                layer.fields + ("geom_prepared",)
                if layer.fields
                else ("geom_prepared",)
            )
            # limit feature number if limit provided
            limit = layer.queryset_limit
            if limit:
                queryset = queryset[:limit]
            # keep values to include in tile (extra included_fields + geometry)
            queryset = queryset.values(*fields)
            # generate SQL/Params for layer
            sql, p = queryset.query.sql_with_params()
            queries.append(
                f"SELECT ST_ASMVT(subquery.*, '{layer.layer_name}', {self.vector_tile_extent}, 'geom_prepared') as mvt "
                f"FROM ({sql}) as subquery"
            )
            params = params + list(p)

        with connection.cursor() as cursor:
            cursor.execute(
                "SELECT string_agg(mvt, '') as mvt FROM ({}) sub".format(' UNION '.join(queries)),
                params=[*params]
            )
            row = cursor.fetchone()[0]
            # psycopg2 returns memoryview, psycopg returns bytes
            return row.tobytes() if type(row) == memoryview else row or b""

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or request

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions