Split tiles on-the-fly to improve over scale vector handling

[HarelM]

Launch Checklist

• Related to Proposal: Better support for overzooming via very large tile extents #2507
• Fixes Safari crashes when zoomed in and there are a lot of labels #6524

This is an experimental feature allowing to over-scale vector tiles by splitting them instead of over scaling them.
It is only used by default for Safari because it changes a bit the results of query render features geometry, since the tiles are no longer capped at max zoom, so the geometry is sliced as well in high zoom tiles, thus creating a "sliced" geometries, while this might be expected up to max-zoom, this is might not be expected in higher zoom levels.
It also changes the rendering of labels in the middle of areas for the same reasons.
The performance of this seems to work well and it solves the crash on iOS as far as my tests went.

What is in this PR:

Once an over-scaled tile request is sent, the code slices the original tile and returns it.
This is very similar to how geojson works, but geojson works only on a single layer, while this is working on multiple layers and doesn't create an index, just a small cache.
Some of the code here was copied from geojson-vt package.
The experimentalOverzoomingWithGeojsonVt flag was added to allow controlling this behavior and avoid backward compatibility issues.

• We should consider removing it as part of v6 breaking changes features #6427, maybe

It is applied by default to Safari browser to avoid the crash that is linked at the beginning of this post.

CC: @wayofthefuture @ChrisLoer

• Confirm your changes do not include backports from Mapbox projects (unless with compliant license) - if you are not sure about this, please ask!
• Briefly describe the changes in this PR.
• Link to related issues.
• Write tests for all new functionality.
• Add an entry to CHANGELOG.md under the ## main section.

[HarelM]

All the failing render test are failing due to placement of labels, which is expected as part of this PR.
I've looks at all the labels diff and I can't say which one is more "correct" between the old and the new, so I tend to say that this is OK.
The query tests are failing too, from what I looked on a single test the query geometry is returned as a sliced geometry instead of the full one, which I think is also expected as part of this change, but can be considered a breaking change for anyone that expects to get the geometry from the tile at max zoom and not a further down sliced geometry.

[HarelM]

I've enabled debug in geojson-vt and I see that it creates indexes from z0 to z14 where most of the data is duplicated, and there's a lot of points there (although nothing it shouldn't handle I believe).
I'll try and see if there's a way to maybe avoid the toGeoJSON part of the code and back to vector tiles and maybe keep everything in tile coordinates somehow...

[HarelM]

So I wrote the following code to "simply" split a tile and take the relevant geometry out of it.
It seems to prevent the crash on iOS.
I copied some of the code from geojson-vt related to slicing/clipping, so @wayofthefuture if you plan to update the code there it might be worth considering exposing some of the clipping methods, maybe.
For some reason it's not perfect and I'm not sure why, I see both contour labels in the wrong place and some line are not scaled correctly (see image below bottom left corner - orange lines), I'm not sure what I did wrong and I'm not sure it's worth the debugging, so I'm leaving it for now.

import {type VectorTile, VectorTileFeature, VectorTileLayer} from '@mapbox/vector-tile';
import Protobuf from 'pbf';
import Point from '@mapbox/point-geometry';
import {fromVectorTileJs} from '@maplibre/vt-pbf';
import {type LoadVectorTileResult} from './vector_tile_worker_source';
import { CanonicalTileID } from './tile_id';

class OverzoomedFeatureWrapper extends VectorTileFeature {
    myGeometry: Point[][];

    constructor(type: 0 | 1 | 2 | 3, geometry: Point[][], properties: any, id: number, extent: number) {
        super(new Protobuf(), 0, extent, [], []);
        this.type = type;
        this.properties = properties ? properties : {};
        this.extent = extent;
        this.myGeometry = geometry;
        this.id = id;
    }

    loadGeometry() {
        return structuredClone(this.myGeometry);
    }
}

class OverzoomedTileLayer extends VectorTileLayer {
    private _myFeatures: OverzoomedFeatureWrapper[];
    name: string;
    extent: number;
    version: number = 2;
    length: number;

    constructor(features: OverzoomedFeatureWrapper[], layerName: string, extent: number) {
        super(new Protobuf());
        this._myFeatures = features;
        this.name = layerName;
        this.length = features.length;
        this.extent = extent;
    }

    feature(i: number): VectorTileFeature {
        return this._myFeatures[i];
    }
}

export class OverzoomedVectorTile implements VectorTile {
    layers: Record<string, VectorTileLayer> = {};

    addLayer(layer: OverzoomedTileLayer) {
        this.layers[layer.name] = layer;
    }
}

/**
 * Encodes the virtual tile into binary vector tile form. 
 * This is a convenience that allows `FeatureIndex` to operate the same way across `VectorTileSource` and `GeoJSONSource` data.
 * @param virtualVectorTile - a VectorTile created from GeoJSON data using geojson-vt
 * @returns 
 */
export function toVirtualVectorTile(virtualVectorTile: VectorTile): LoadVectorTileResult {
    
    let pbf: Uint8Array = fromVectorTileJs(virtualVectorTile);
    if (pbf.byteOffset !== 0 || pbf.byteLength !== pbf.buffer.byteLength) {
        pbf = new Uint8Array(pbf);  // Compatibility with node Buffer (https://github.com/mapbox/pbf/issues/35)
    }
    
    return {
        vectorTile: virtualVectorTile,
        rawData: pbf.buffer
    };
}

export function sliceTileLayer(sourceLayer: VectorTileLayer, maxZoomTileID: CanonicalTileID, targetTileID: CanonicalTileID): OverzoomedTileLayer {
    const { extent } = sourceLayer;
    if (extent !== 4096) {
        console.log(`Warning: overzooming only tested with extent=4096, got extent=${extent}`);
    }
    const dz = targetTileID.z - maxZoomTileID.z;
    const scale = Math.pow(2, dz);
    const startXRelative = targetTileID.x / scale - maxZoomTileID.x;
    const startYRelative = targetTileID.y / scale - maxZoomTileID.y;

    const clipMinX = Math.round(startXRelative * extent);
    const clipMinY = Math.round(startYRelative * extent);
    const clipMaxX = Math.round((startXRelative + (1 / scale)) * extent);
    const clipMaxY = Math.round((startYRelative + (1 / scale)) * extent);

    // The clipping box coordinates within the source tile (0-extent)
    const clippingBox = {
        minX: clipMinX - 128,
        minY: clipMinY - 128,
        maxX: clipMaxX + 128,
        maxY: clipMaxY + 128
    };

    const featureWrappers: OverzoomedFeatureWrapper[] = [];
    for (let index = 0; index < sourceLayer.length; index++) {
        const feature: VectorTileFeature = sourceLayer.feature(index);
        let geometry = feature.loadGeometry();
        geometry = clip(geometry, feature.type, clippingBox.minX, clippingBox.maxX, AxisType.X);
        geometry = clip(geometry, feature.type, clippingBox.minY, clippingBox.maxY, AxisType.Y);
        // scale coordinates to the new tile extent
        for (const ring of geometry) {
            for (const point of ring) {
                point.x = (point.x - clipMinX) * scale;
                point.y = (point.y - clipMinY) * scale;
            }
        }
        if (geometry.length === 0) {
            continue;
        }
        featureWrappers.push(new OverzoomedFeatureWrapper(
            feature.type,
            geometry,
            feature.properties,
            feature.id,
            extent
        ));
    }
    return new OverzoomedTileLayer(featureWrappers, sourceLayer.name, extent);
}

const enum AxisType {
    X = 0,
    Y = 1
}

/* clip features between two vertical or horizontal axis-parallel lines:
 *     |        |
 *  ___|___     |     /
 * /   |   \____|____/
 *     |        |
 *
 * k1 and k2 are the line coordinates
 * axis: 0 for x, 1 for y
 * minAll and maxAll: minimum and maximum coordinate value for all features
 */
export default function clip(geometry: Point[][], type: number, start: number, end: number, axis: AxisType): Point[][] {
    //const min = axis === AxisType.X ? geometry.map(ring => ring.reduce((min, p) => Math.min(min, p.x), Infinity)).reduce((a, b) => Math.min(a, b), Infinity) : geometry.map(ring => ring.reduce((min, p) => Math.min(min, p.y), Infinity)).reduce((a, b) => Math.min(a, b), Infinity);
    //const max = axis === AxisType.Y ? geometry.map(ring => ring.reduce((max, p) => Math.max(max, p.y), -Infinity)).reduce((a, b) => Math.max(a, b), -Infinity) : geometry.map(ring => ring.reduce((max, p) => Math.max(max, p.x), -Infinity)).reduce((a, b) => Math.max(a, b), -Infinity);

    //if (max < start || min >= end) { // trivial reject
    //    return [];
    //}
    switch (type) {
        case 1: // POINT
            return clipPoints(geometry, start, end, axis);
        case 2: // LINESTRING
            return clipLines(geometry, start, end, axis, false);
        case 3: // POLYGON
            return clipLines(geometry, start, end, axis, true);
    }
    
    return null;
}

function clipPoints(geometry: Point[][], start: number, end: number, axis: AxisType): Point[][] {
    const newGeometry: Point[][] = [];
    for (const ring of geometry) {
        for (const point of ring) {
            const a = axis === AxisType.X ? point.x : point.y;
            if (a >= start && a <= end) {
                newGeometry.push([point]);
            }
        }
    }
    return newGeometry;
}

function clipLine(line: Point[], start: number, end: number, axis: AxisType, isPolygon: boolean): Point[][] {
    const intersect = axis === AxisType.X ? intersectX : intersectY;

    let slice: Point[] = [];
    const newLine: Point[][] = [];
    for (let i = 0; i < line.length - 1; i++) {
        const p1 = line[i];
        const p2 = line[i + 1];
        const a = axis === AxisType.X ? p1.x : p1.y;
        const b = axis === AxisType.X ? p2.x : p2.y;
        let exited = false;

        if (a < start) {
            // ---|-->  | (line enters the clip region from the left)
            if (b > start) {
                slice.push(intersect(p1, p2, start));
            }
        } else if (a > end) {
            // |  <--|--- (line enters the clip region from the right)
            if (b < end) {
                slice.push(intersect(p1, p2, end));
            }
        } else {
            slice.push(p1);
        }
        if (b < start && a >= start) {
            // <--|---  | or <--|-----|--- (line exits the clip region on the left)
            slice.push(intersect(p1, p2, start));
            exited = true;
        }
        if (b > end && a <= end) {
            // |  ---|--> or ---|-----|--> (line exits the clip region on the right)
            slice.push(intersect(p1, p2, end));
            exited = true;
        }

        if (!isPolygon && exited) {
            newLine.push(slice);
            slice = [];
        }
    }

    // add the last point
    let last = line.length - 1;
    const a = axis === AxisType.X ? line[last].x : line[last].y;
    if (a >= start && a <= end) {
        slice.push(line[last]);
    }

    // close the polygon if its endpoints are not the same after clipping
    if (isPolygon && slice.length > 0 && !slice[0].equals(slice[slice.length - 1])) {
        slice.push(new Point(slice[0].x, slice[0].y));
    }
    if (slice.length > 0) {
        newLine.push(slice);
    }
    return newLine;
}


function clipLines(geometry: Point[][], start: number, end: number, axis: AxisType, isPolygon: boolean): Point[][] {
    const newGeometry: Point[][] = [];
    for (const line of geometry) {
        const clippedLines = clipLine(line, start, end, axis, isPolygon);
        if (clippedLines.length > 0) {
            newGeometry.push(...clippedLines);
        }
    }
    return newGeometry;
}

function intersectX(p1: Point, p2: Point, x: number): Point {
    const t = (x - p1.x) / (p2.x - p1.x);
    return new Point(x, p1.y + (p2.y - p1.y) * t);
}

function intersectY(p1: Point, p2: Point, y: number): Point {
    const t = (y - p1.y) / (p2.y - p1.y);
    return new Point(p1.x + (p2.x - p1.x) * t, y);
}

[wayofthefuture]

vector_tile_worker_source.ts

overzoomed-geojson-vector-tile.ts

[HarelM]

Nice! Thanks for this! So this was about when to do the coordinates conversion basically?
Are you sure about the need to cache the calculated tiles? I'm considering removing the cache all together...?

[wayofthefuture]

Yes... I didn't get a chance to calculate performance but processing could be heavy and it is a nice o(1) lru cache that we can also use in source cache eventually. It's simple and I don't see how it wouldn't improve performance.

[HarelM]

I think the problem here is that iOS might kill things that use too much memory, this is what I think is causing the crash with the geojson index, although it might be a magnitude bigger...

[wayofthefuture]

Many magnitudes larger especially with tolerance: 0. I was able to monitor memory with the worker, and the previous version was growing well over 10GB in memory and quickly. The current version with 10,000 (notice I added a 0!) in cache didn't exceed 2GB with heavy automated movement. With 0 cache the page was about 1GB in total. You had 1 pyramid per layer per max zoom tile, that's like n^n^n or something! The current cache is simply flat vector tiles so I think it will work nicely.

[HarelM]

I found this when searching for something related to clipping in the code base as I assumed there should already be something:
https://github.com/maplibre/maplibre-gl-js/blob/main/src/symbol/clip_line.ts

It's not exactly the same, but looks very similar, I wonder if I can use it or adopt it to our needs...

[HarelM]

I moved the clipping code to the same file as the other method.
They are very similar, but they differ a bit with the rounding of the intersection points, so I'll keep them side by side at the moment.
I also changed the name of the experimental variable to match the current implementation.

[andrewharvey]

Thanks for your work on this @HarelM.

It is only used by default for Safari because it changes a bit the results of query render features geometry, since the tiles are no longer capped at max zoom, so the geometry is sliced as well in high zoom tiles, thus creating a "spliced" geometries, while this might be expected up to max-zoom, this is might not be expected in higher zoom levels.
The performance of this seems to work well and it solves the crash on iOS as far as my tests went.

I'm not a fan of having the library function differently on different browsers like this. We should set a clear expectation to users about what queryRenderedFeatures returns. We already say in the API docs

Because features come from tiled vector data or GeoJSON data that is converted to tiles internally, feature geometries may be split or duplicated across tile boundaries and, as a result, features may appear multiple times in query results. For example, suppose there is a highway running through the bounding rectangle of a query. The results of the query will be those parts of the highway that lie within the map tiles covering the bounding rectangle, even if the highway extends into other tiles, and the portion of the highway within each map tile will be returned as a separate feature. Similarly, a point feature near a tile boundary may appear in multiple tiles due to tile buffering.

I think we should just accept that geometries returned from queryRenderedFeatures at an overzoom will have further clipping applied beyond the sources maxZoom. That should still represent the geometries visible on the map and therefore works the same as any other zoom level under the maxZoom of the source.

[HarelM]

I'm not a fan of having the library function differently on different browsers like this

Neither am I but currently Safari crashes while Chrome doesn't, so the difference is already there to some extent.

The main difference I saw is for areas that has a label placed in the middle, this label is now repeated in every sliced tile, it becomes apparent in high zoom levels (20+) since this is most of what you see in the screen, since there's not a lot of data in one screen.

I suggest to try the branch out with a style your are familiar with and see if you can spot the changes.

[HarelM]

Ok, So I wanted to polish this a bit more and I think I came up with an idea to set the number of zoom level to overscale instead of doing it for all zoom levels.
This would mean that you can specify how many zoom levels to overscale and the rest would be using the split technique.
I'm currently testing a default of 4 - meaning, the last 4 zoom levels will be overscaled and the rest will be sliced.
This gives more control on when to overscale and when to slice, which I think is a good middle ground.

[wayofthefuture]

So if you have map tile boundaries enabled would you see it switch from 1 tile to 256 subdivided tiles at maxzoom + 4 -> maxzoom + 5? Not sure how such a switch would look and it could have drastic style changes maybe?

[HarelM]

I'll push it to a different branch so you could play with it if you want.
Let's use number to illustrate it better:

Source max zoom is 10.
Map max zoom is 20.
We are over-scaling for only 3 zoom levels.

In that case the slicing will happen for zoom 11 to 17, and after that, a level 17 "sliced tile" will be used as the base for over scaling, similar to what happens now with the tiles at max zoom, only in this implementation we first slice the tile and later on, if needed scale it.

It gives more control on when to stop slicing and start scaling.
A private case would be the behavior before this PR which is to always scale.
It also removes the "bypass covering tiles" part by using the maxzoom parameter.

[HarelM]

See the following PR:

• Define how much to slice and how much to scale #6617

[wayofthefuture]

Got it I was thinking of it backwards...