Web Mercator Tile Visibility
This notebook demonstrates a method to work out what Web Mercator tiles (also known as Slippy tiles) are visible on an arbitrary map projection P. Since slippy tiles are arranged in a quadtree, it's possible to traverse tiles recursively to determine visibility. This enables fast visibility detection, even at high zoom levels (where there may be thousands of candidate tiles). Inspired by Jason Davies' work in the same domain. Click and drag the globe!
demo = { let demo = DOM.element('div'); demo.appendChild(globeView.canvas); demo.appendChild(mercatorView.canvas); globeView.canvas.style.display = 'inline-block'; globeView.canvas.style.width = '45%'; globeView.canvas.style.margin = '5px'; mercatorView.canvas.style.display = 'inline-block'; mercatorView.canvas.style.width = '45%'; mercatorView.canvas.style.margin = '5px'; return demo; }
viewof scale = { const input = DOM.element('input'); input.type = 'range'; input.min = 1; input.max = 10; input.step = 'any'; input.value = 1; input.addEventListener('input', () => { mutable globe = mutable globe.scale(300 * input.value); }); return input; }
viewof zoom = { const scaleToZoom = scale => Math.round(2 * Math.log2(scale + 1)); let form; if (this) { form = this; if (form.auto.checked) { form.zoom.value = scaleToZoom(scale); } form.value = form.zoom.valueAsNumber; } else { form = html`<form> <input type='range' name='zoom' min=0 max=10 value=4></input> <input type='checkbox' name='auto'></input> <label>auto</label> </form>`; form.zoom.oninput = () => { form.value = form.zoom.valueAsNumber; }; form.value = form.zoom.valueAsNumber; } form.auto.oninput = () => { form.zoom.disabled = form.auto.checked; form.zoom.value = scaleToZoom(scale); form.value = form.zoom.valueAsNumber; }; return form; }
In the demo above, the yellow rectangle on the globe is the area of interest; slippy tiles that intersect this area are shown in red. Try dragging the globe to rotate the Earth, adjusting the scale slider (which will zoom in on the globe), or adjusting the zoom slider (which will change the size of the mercator tiles being considered). When "auto" is checked, an appropriate zoom level is chosen based on the current scale.
The heavy lifting here is done in getVisibleTiles (below), which is commented to help explain the process.
getVisibleTiles = function(projection, zoom) { function getVisibleChildren(x, y, z) { console.log(x, y, z); let visible = false; // Define a projection stream based on the input projection. This stream first does all of // the operations (rotating, clipping, and resampling) that the input projection would do. // Then, any geometries that survive that process (i.e. are visible in the output projection) // get passed on to our custom stream, which just sets `visible = true` if it sees a point. // We'll use this stream to determine if a tile is visible when projected through `projection`. let stream = projection.stream({ point: () => { visible = true; }, lineStart: () => {}, lineEnd: () => {}, polygonStart: () => {}, polygonEnd: () => {} }); // calculate screen-space coordinates of tile (x, y, z) in the `mercator` projection below const extent = mercator.clipExtent(); // we assume that extent is a square, and thus all tiles are also squares const size = (extent[1][0] - extent[0][0]) / (2 ** z); const epsilon = 0.001; const x0 = x * size + epsilon, y0 = y * size + epsilon, x1 = (x + 1) * size - epsilon, y1 = (y + 1) * size - epsilon; const step = Math.min(size, 1); const tile = { type: "Polygon", coordinates: [ [] // Interpolate along each of the four edges of the Mercator polygon. This is a rather // crude way to do interpolation (with a constant step size) but it's easier than the // adaptive resampling that d3 does natively, which is not available to us here since // we're not using a projection. .concat(d3.range(x0, x1 + step / 2, +step).map(x => [x, y0])) .concat(d3.range(y0, y1 + step / 2, +step).map(y => [x1, y])) .concat(d3.range(x1, x0 - step / 2, -step).map(x => [x, y1])) .concat(d3.range(y1, y0 - step / 2, -step).map(y => [x0, y])) // convert this mercator polygon to a spherical polygon on the Earth .map(point => mercator.invert(point)) ], }; // Project the spherical polygon representing the tile by feeding it through the stream. d3.geoStream(tile, stream); // If any point in the tile was visible in the output projection, then our stream will have // received that point via its `point` method and set `visible = true`. if (visible) { if (z < zoom) { // if we're not yet at the desired zoom level, recurse on each of the four sub-tiles let visibleChildren = [] .concat(getVisibleChildren(2 * x, 2 * y, z + 1)) // upper left .concat(getVisibleChildren(2 * x + 1, 2 * y, z + 1)) // upper right .concat(getVisibleChildren(2 * x, 2 * y + 1, z + 1)) // lower left .concat(getVisibleChildren(2 * x + 1, 2 * y + 1, z + 1)) // lower right return visibleChildren; } else { return [tile]; // this tile is visible and is of the desired zoom level } } else { return []; // this tile is not visible } } return getVisibleChildren(0, 0, 0); }
globeView = { let context; if (this) { context = this; //this.getContext("2d"); context.clearRect(0, 0, size, size); } else { context = DOM.context2d(size, size); d3.select(context.canvas).call(drag); } const path = d3.geoPath(globe, context); context.strokeStyle = 'black'; context.fillStyle = 'black'; context.lineWidth = 1; context.beginPath(), path(sphere), context.fillStyle = "#fff", context.fill(); context.beginPath(), path(land), context.fillStyle = "#000", context.fill(); context.beginPath(), path(graticule), context.strokeStyle = 'rgba(100, 100, 100, 0.25)', context.stroke(); context.beginPath(), path(quads), context.strokeStyle = 'red', context.stroke(); context.beginPath(), path(quads), context.fillStyle = 'rgba(255, 0, 0, 0.25)', context.fill(); context.beginPath(), path(sphere), context.strokeStyle = 'black', context.stroke(); context.strokeStyle = 'yellow'; context.strokeRect(rect.x, rect.y, rect.width, rect.height); // context.fillStyle = 'rgba(0, 0, 0, 0.25)'; // context.fillRect(0, 0, size, rect.y); // context.fillRect(0, rect.y, rect.x, rect.height); // context.fillRect(rect.x + rect.width, rect.y, size, rect.height); // context.fillRect(0, rect.y + rect.height, size, size); return context; }
mercatorView = { let context; if (this) { context = this; //this.getContext("2d"); context.clearRect(0, 0, size, size); } else { context = DOM.context2d(size, size); } const path = d3.geoPath(mercator, context); context.fillStyle = 'black'; context.strokeStyle = 'black'; context.beginPath(), path(sphere), context.fillStyle = "#fff", context.fill(); context.beginPath(), path(land), context.fillStyle = "#000", context.fill(); context.beginPath(), path(graticule), context.strokeStyle = 'rgba(100, 100, 100, 0.25)', context.stroke(); for (var quad of quads.features) { context.beginPath(), path(quad), context.strokeStyle = 'red', context.stroke(); context.beginPath(), path(quad), context.fillStyle = 'rgba(255, 0, 0, 0.25)', context.fill(); } context.beginPath(), path(sphere), context.strokeStyle = 'black', context.stroke(); return context; }
quads = { globe.clipExtent([[rect.x, rect.y], [rect.x + rect.width, rect.y + rect.height]]); let features = getVisibleTiles(globe, zoom) .map(geom => turf.feature(geom)); globe.clipExtent([[0, 0], [size, size]]); return turf.featureCollection(features); }
drag = { let v0, q0, r0; function dragstarted() { v0 = versor.cartesian(mutable globe.invert(d3.mouse(this))); q0 = versor(r0 = mutable globe.rotate()); } function dragged() { const v1 = versor.cartesian(mutable globe.rotate(r0).invert(d3.mouse(this))); const q1 = versor.multiply(q0, versor.delta(v0, v1)); mutable globe = mutable globe.rotate(versor.rotation(q1)); } return d3.drag() .on("start", dragstarted) .on("drag", dragged); }
rect = Object({ x: 1/6 * size, y: 2/6 * size, width: 4/6 * size, height: 2/6 * size, })
size = 600
mutable globe = d3.geoOrthographic() .translate([size / 2, size / 2]) .clipExtent([[0, 0], [size, size]]) .rotate([-17, -13, 30]) .precision(0.1)
mercator = d3.geoMercator() .fitSize([size, size], sphere) .clipExtent([[0, 0], [size, size]])
graticule = d3.geoGraticule10()
sphere = Object({ type: 'Sphere' })
land = topojson.feature(world, world.objects.land)
world = d3.json("https://unpkg.com/world-atlas@1/world/110m.json")
versor = { var acos = Math.acos, asin = Math.asin, atan2 = Math.atan2, cos = Math.cos, max = Math.max, min = Math.min, PI = Math.PI, sin = Math.sin, sqrt = Math.sqrt, radians = PI / 180, degrees = 180 / PI; // Returns the unit quaternion for the given Euler rotation angles [λ, φ, γ]. function versor(e) { var l = e[0] / 2 * radians, sl = sin(l), cl = cos(l), // λ / 2 p = e[1] / 2 * radians, sp = sin(p), cp = cos(p), // φ / 2 g = e[2] / 2 * radians, sg = sin(g), cg = cos(g); // γ / 2 return [ cl * cp * cg + sl * sp * sg, sl * cp * cg - cl * sp * sg, cl * sp * cg + sl * cp * sg, cl * cp * sg - sl * sp * cg ]; } // Returns Cartesian coordinates [x, y, z] given spherical coordinates [λ, φ]. versor.cartesian = function(e) { var l = e[0] * radians, p = e[1] * radians, cp = cos(p); return [cp * cos(l), cp * sin(l), sin(p)]; }; // Returns the Euler rotation angles [λ, φ, γ] for the given quaternion. versor.rotation = function(q) { return [ atan2(2 * (q[0] * q[1] + q[2] * q[3]), 1 - 2 * (q[1] * q[1] + q[2] * q[2])) * degrees, asin(max(-1, min(1, 2 * (q[0] * q[2] - q[3] * q[1])))) * degrees, atan2(2 * (q[0] * q[3] + q[1] * q[2]), 1 - 2 * (q[2] * q[2] + q[3] * q[3])) * degrees ]; }; // Returns the quaternion to rotate between two cartesian points on the sphere. versor.delta = function(v0, v1) { var w = cross(v0, v1), l = sqrt(dot(w, w)); if (!l) return [1, 0, 0, 0]; var t = acos(max(-1, min(1, dot(v0, v1)))) / 2, s = sin(t); // t = θ / 2 return [cos(t), w[2] / l * s, -w[1] / l * s, w[0] / l * s]; }; // Returns the quaternion that represents q0 * q1. versor.multiply = function(q0, q1) { return [ q0[0] * q1[0] - q0[1] * q1[1] - q0[2] * q1[2] - q0[3] * q1[3], q0[0] * q1[1] + q0[1] * q1[0] + q0[2] * q1[3] - q0[3] * q1[2], q0[0] * q1[2] - q0[1] * q1[3] + q0[2] * q1[0] + q0[3] * q1[1], q0[0] * q1[3] + q0[1] * q1[2] - q0[2] * q1[1] + q0[3] * q1[0] ]; }; function cross(v0, v1) { return [ v0[1] * v1[2] - v0[2] * v1[1], v0[2] * v1[0] - v0[0] * v1[2], v0[0] * v1[1] - v0[1] * v1[0] ]; } function dot(v0, v1) { return v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2]; } return versor; }
topojson = require("topojson-client@3")
turf = require("@turf/turf@5")
d3 = require("d3-drag", "d3-array", "d3-tile", "d3-fetch@0.1", "d3-geo@1", "d3-geo-projection", "d3-selection", "d3-color")
pi = Math.PI
tau = 2 * Math.PI