feat(geo): vector tiles of features for maps

baseapp/geo/README.md

@@ -80,6 +80,79 @@ Malformed filter values (wrong arity, out-of-range coordinates, radius over the
 cap, unknown target IDs) raise field-scoped validation errors instead of
 silently returning unfiltered results.
 
+### Vector tiles (maps)
+
+`geoFeatures` returns at most 100 features per page, which is too few to draw a
+dense map. For maps, the plugin also serves features as
+[Mapbox Vector Tiles](https://github.com/mapbox/vector-tile-spec) (MVT), built by
+PostGIS (`ST_AsMVT`), through the project's `v1` REST URLs:
+
+```
+GET /v1/geo/tiles/{z}/{x}/{y}.mvt
+```
+
+- **From zoom 13 up**, the tile's `features` layer holds every feature in it,
+  points as points and polygons as polygons, with `id` (the feature's public ID,
+  the one the GraphQL API uses) and `feature_type` attributes.
+- **Below zoom 13**, the `clusters` layer aggregates features into a 16×16 grid
+  per tile: one point per non-empty cell, at the centroid of its features, with
+  a `point_count`. Each feature counts in exactly one tile, so totals add up.
+- **Filters**: the same query parameters as `geoFeatures` (`feature_type`,
+  `target_object_id`, `bbox`, `near`), e.g. `?feature_type=store`. Malformed
+  values return `400` with field errors.
+- **Access**: requires `view_geojsonfeature`, through `has_perm` (public with the
+  default backend).
+- **Responses**: empty tiles return `204`. Tiles carry
+  `Cache-Control: public, max-age=300`, so a CDN can absorb map pans.
+- **Cost of zooms 0–3**: those 85 tiles skip the spatial index (on tiles that
+  large, a geography envelope's great-circle edges diverge from the map's), so
+  each one scans every feature the filters allow. Fine for tens of thousands of
+  features; on much larger tables, cache those tiles in front of the API or set
+  the map's `minzoom` to 4.
+
+With MapLibre GL:
+
+```js
+map.addSource("features", {
+  type: "vector",
+  tiles: [`${API_URL}/v1/geo/tiles/{z}/{x}/{y}.mvt?feature_type=store`],
+});
+map.addLayer({
+  id: "clusters",
+  type: "circle",
+  source: "features",
+  "source-layer": "clusters",
+  paint: { "circle-radius": ["interpolate", ["linear"], ["get", "point_count"], 1, 4, 100, 18] },
+});
+map.addLayer({ id: "features", type: "circle", source: "features", "source-layer": "features" });
+```
+
+To add attributes or filters, subclass the view and route it in your project:
+
+```python
+# myproject/geo/views.py
+from django.db.models import F, Q
+from baseapp.geo.rest_framework.views import GeoJSONFeatureTileView
+
+
+class StoreTileView(GeoJSONFeatureTileView):
+    cluster_max_zoom = 12  # individual features from zoom 12
+    cluster_grid_size = 8  # fewer, larger clusters
+
+    def filter_queryset(self, queryset):
+        queryset = super().filter_queryset(queryset)
+        if self.request.query_params.get("open_now") == "true":
+            queryset = queryset.filter(Q(name__icontains="24h"))
+        return queryset
+
+    def get_feature_properties(self):
+        # Values must be strings, numbers or booleans.
+        return {**super().get_feature_properties(), "name": F("name")}
+```
+
+`baseapp.geo.tiles.build_features_tile` / `build_clusters_tile` build tiles from
+any `GeoJSONFeature` queryset, for projects that need their own endpoint.
+
 ### Geometry input
 
 Mutations accept geometry through the `Geometry` scalar as any of:

baseapp/geo/plugin.py

@@ -25,7 +25,13 @@ def get_settings(self) -> PackageSettings:
             graphql_mutations=[
                 "baseapp.geo.graphql.mutations.GeoMutations",
             ],
+            # REST: vector tiles for maps
+            v1_urlpatterns=self.v1_urlpatterns,
             # Deps
             required_packages=[],
             optional_packages=[],
         )
+
+    @staticmethod
+    def v1_urlpatterns(include, path, re_path) -> list:
+        return [re_path(r"", include("baseapp.geo.rest_framework.urls"))]

baseapp/geo/rest_framework/permissions.py

@@ -0,0 +1,14 @@
+import swapper
+from rest_framework.permissions import BasePermission
+
+
+class CanViewGeoJSONFeatures(BasePermission):
+    """Read access to features, as decided by the `view_geojsonfeature` permission.
+
+    Goes through `has_perm`, so a project that overrides `GeoPermissionsBackend` (e.g. to hide
+    features from anonymous users) gets the same rule on the REST endpoints as on GraphQL.
+    """
+
+    def has_permission(self, request, view) -> bool:
+        GeoJSONFeature = swapper.load_model("baseapp_geo", "GeoJSONFeature")
+        return request.user.has_perm(f"{GeoJSONFeature._meta.app_label}.view_geojsonfeature")

baseapp/geo/rest_framework/urls.py

@@ -0,0 +1,11 @@
+from django.urls import path
+
+from .views import GeoJSONFeatureTileView
+
+urlpatterns = [
+    path(
+        "geo/tiles/<int:z>/<int:x>/<int:y>.mvt",
+        GeoJSONFeatureTileView.as_view(),
+        name="geo-feature-tiles",
+    ),
+]

baseapp/geo/rest_framework/views.py

@@ -0,0 +1,96 @@
+import swapper
+from django.core.exceptions import ValidationError
+from django.db.models import CharField, Expression, F, QuerySet
+from django.db.models.functions import Cast
+from django.http import Http404, HttpResponse
+from django.utils.translation import gettext_lazy as _
+from rest_framework import status
+from rest_framework.request import Request
+from rest_framework.response import Response
+from rest_framework.views import APIView
+
+from baseapp.geo.graphql.filters import GeoJSONFeatureFilter
+from baseapp.geo.tiles import Tile, build_clusters_tile, build_features_tile
+
+from .permissions import CanViewGeoJSONFeatures
+
+MVT_CONTENT_TYPE = "application/vnd.mapbox-vector-tile"
+
+
+class GeoJSONFeatureTileView(APIView):
+    """Mapbox Vector Tile of GeoJSON features: `GET .../geo/tiles/{z}/{x}/{y}.mvt`.
+
+    From `cluster_max_zoom` up, the tile's `features` layer holds every feature, with the
+    attributes of `get_feature_properties()` (by default its public `id` and `feature_type`). Below it, the `clusters`
+    layer aggregates them into grid cells with a `point_count`. Query parameters are the same
+    filters as the `geoFeatures` GraphQL connection (e.g. `feature_type`, `target_object_id`).
+
+    Extension points for projects: `get_queryset`, `filter_queryset` and
+    `get_feature_properties` (e.g. to add attributes of the target object).
+    """
+
+    permission_classes = [CanViewGeoJSONFeatures]
+    filterset_class = GeoJSONFeatureFilter
+
+    features_layer = "features"
+    clusters_layer = "clusters"
+    cluster_max_zoom = 13
+    cluster_grid_size = 16
+    # Short enough for new features to show up soon, long enough for a CDN to absorb map pans.
+    cache_max_age = 300
+
+    def get_queryset(self) -> QuerySet:
+        return swapper.load_model("baseapp_geo", "GeoJSONFeature").objects.all()
+
+    def filter_queryset(self, queryset: QuerySet) -> QuerySet:
+        """Apply `filterset_class` to the query parameters; raises ValidationError on bad input."""
+        filterset = self.filterset_class(self.request.query_params, queryset=queryset)
+        if not filterset.is_valid():
+            raise ValidationError(filterset.errors)
+        return filterset.qs
+
+    def get_feature_properties(self) -> dict[str, Expression]:
+        """MVT attributes of each feature: name → ORM expression (string, number or boolean).
+
+        `id` is the feature's public ID (its DocumentId's UUID, as text): the same ID the
+        GraphQL API uses, and it doesn't reveal how many features exist. `document` joins at
+        most one DocumentId per feature (unique per content type and object).
+        """
+        return {
+            "id": Cast("document__public_id", output_field=CharField()),
+            "feature_type": F("feature_type"),
+        }
+
+    def get(self, request: Request, z: int, x: int, y: int) -> HttpResponse:
+        tile = Tile(z, x, y)
+        if not tile.is_valid:
+            raise Http404(_("Tile out of range."))
+
+        try:
+            queryset = self.filter_queryset(self.get_queryset())
+            prefilter = tile.prefilter_polygon()
+            if prefilter is not None:
+                queryset = queryset.filter(geometry__intersects=prefilter)
+
+            if z < self.cluster_max_zoom:
+                content = build_clusters_tile(
+                    tile, queryset, grid_size=self.cluster_grid_size, layer=self.clusters_layer
+                )
+            else:
+                content = build_features_tile(
+                    tile,
+                    queryset,
+                    properties=self.get_feature_properties(),
+                    layer=self.features_layer,
+                )
+        except ValidationError as exc:
+            detail = exc.message_dict if hasattr(exc, "error_dict") else {"detail": exc.messages}
+            return Response(detail, status=status.HTTP_400_BAD_REQUEST)
+
+        response = HttpResponse(
+            content,
+            content_type=MVT_CONTENT_TYPE,
+            status=status.HTTP_200_OK if content else status.HTTP_204_NO_CONTENT,
+        )
+        response["Cache-Control"] = f"public, max-age={self.cache_max_age}"
+        return response

baseapp/geo/tests/mvt.py

@@ -0,0 +1,135 @@
+"""Minimal Mapbox Vector Tile decoder for tests.
+
+Reads the protobuf wire format of the MVT 2.1 spec directly, so tests need no protobuf-based
+dependency (mapbox-vector-tile caps protobuf below the version baseapp locks). Decodes layers,
+their extent, and each feature's type, properties and geometry in tile coordinates.
+https://github.com/mapbox/vector-tile-spec/tree/master/2.1
+"""
+
+import struct
+from collections.abc import Iterator
+
+GEOMETRY_TYPES = {1: "Point", 2: "LineString", 3: "Polygon"}
+MOVE_TO, LINE_TO, CLOSE_PATH = 1, 2, 7
+
+
+def _varint(data: bytes, pos: int) -> tuple[int, int]:
+    result = shift = 0
+    while True:
+        byte = data[pos]
+        pos += 1
+        result |= (byte & 0x7F) << shift
+        if not byte & 0x80:
+            return result, pos
+        shift += 7
+
+
+def _fields(data: bytes) -> Iterator[tuple[int, int, int | bytes]]:
+    """Yield `(field_number, wire_type, value)`; length-delimited values are bytes."""
+    pos = 0
+    while pos < len(data):
+        key, pos = _varint(data, pos)
+        field, wire_type = key >> 3, key & 7
+        if wire_type == 0:
+            value, pos = _varint(data, pos)
+        elif wire_type == 1:
+            value, pos = data[pos : pos + 8], pos + 8
+        elif wire_type == 2:
+            length, pos = _varint(data, pos)
+            value, pos = data[pos : pos + length], pos + length
+        elif wire_type == 5:
+            value, pos = data[pos : pos + 4], pos + 4
+        else:
+            raise ValueError(f"Unsupported wire type {wire_type}")
+        yield field, wire_type, value
+
+
+def _packed(data: bytes) -> list[int]:
+    values, pos = [], 0
+    while pos < len(data):
+        value, pos = _varint(data, pos)
+        values.append(value)
+    return values
+
+
+def _zigzag(value: int) -> int:
+    return (value >> 1) ^ -(value & 1)
+
+
+def _value(data: bytes) -> str | float | int | bool | None:
+    for field, _wire_type, raw in _fields(data):
+        if field == 1:
+            return raw.decode()
+        if field == 2:
+            return struct.unpack("<f", raw)[0]
+        if field == 3:
+            return struct.unpack("<d", raw)[0]
+        if field == 4:  # int64: two's complement in the varint
+            return raw if raw < 2**63 else raw - 2**64
+        if field == 5:  # uint64
+            return raw
+        if field == 6:
+            return _zigzag(raw)
+        if field == 7:
+            return bool(raw)
+    return None
+
+
+def _geometry(commands: list[int]) -> list[list[tuple[int, int]]]:
+    """Decode drawing commands into parts (points, or rings), in tile coordinates."""
+    parts, current, x, y, pos = [], [], 0, 0, 0
+    while pos < len(commands):
+        command, count = commands[pos] & 7, commands[pos] >> 3
+        pos += 1
+        if command == CLOSE_PATH:
+            continue
+        for _ in range(count):
+            x += _zigzag(commands[pos])
+            y += _zigzag(commands[pos + 1])
+            pos += 2
+            if command == MOVE_TO and current:
+                parts.append(current)
+                current = []
+            current.append((x, y))
+    if current:
+        parts.append(current)
+    return parts
+
+
+def _feature(data: bytes, keys: list[str], values: list) -> dict:
+    feature = {"type": None, "properties": {}, "geometry": []}
+    for field, _wire_type, raw in _fields(data):
+        if field == 2:
+            tags = _packed(raw)
+            for key_index, value_index in zip(tags[::2], tags[1::2]):
+                feature["properties"][keys[key_index]] = values[value_index]
+        elif field == 3:
+            feature["type"] = GEOMETRY_TYPES.get(raw)
+        elif field == 4:
+            feature["geometry"] = _geometry(_packed(raw))
+    return feature
+
+
+def decode(tile: bytes) -> dict[str, dict]:
+    """`{layer name: {"extent": int, "features": [feature, ...]}}`."""
+    layers = {}
+    for field, _wire_type, layer_data in _fields(tile):
+        if field != 3:
+            continue
+        name, extent, keys, values, raw_features = None, 4096, [], [], []
+        for layer_field, _layer_wire_type, raw in _fields(layer_data):
+            if layer_field == 1:
+                name = raw.decode()
+            elif layer_field == 2:
+                raw_features.append(raw)
+            elif layer_field == 3:
+                keys.append(raw.decode())
+            elif layer_field == 4:
+                values.append(_value(raw))
+            elif layer_field == 5:
+                extent = raw
+        layers[name] = {
+            "extent": extent,
+            "features": [_feature(raw, keys, values) for raw in raw_features],
+        }
+    return layers