From af2b9f5c0186260f2969450c4303f0605d29d3f0 Mon Sep 17 00:00:00 2001 From: Kohl Listi Date: Thu, 4 Mar 2021 00:28:12 +0100 Subject: [PATCH 1/2] Fix calcZForBBox Static image rendering moves bbox (wrong zoom level by calcZForBBox). This changesets calculates the correct zoom level for a given bounding box. See maptiler/tileserver-gl#490 --- src/serve_rendered.js | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/serve_rendered.js b/src/serve_rendered.js index 11b6bfb..2c84762 100644 --- a/src/serve_rendered.js +++ b/src/serve_rendered.js @@ -174,24 +174,33 @@ const renderOverlay = (z, x, y, bearing, pitch, w, h, scale, }; const calcZForBBox = (bbox, w, h, query) => { - let z = 25; + // Use equation for horizontal distance per tile given zoom and latitude + // see https://wiki.openstreetmap.org/wiki/Zoom_levels#Distance_per_pixel_math + // + // S_pixel * T = S_tile = C * cos(l) / (2^z) = + // S_tile := horizontal distance per tile + // S_pixel := horizontal distance per pixel + // T := tile size + // C := Equatorial circumference of Earth + // l := centered latitude of bbox + // z := zoom level + // + // => z = log2( C * cos(l) / (S_pixel*T) ) - const padding = query.padding !== undefined ? - parseFloat(query.padding) : 0.1; + // points (use centered latitude of bbox) + const latCenter = bbox[1]+(bbox[3]-bbox[1])/2; + const latRadian = latCenter*Math.PI/180; + const east = [bbox[0], latCenter]; + const west = [bbox[2], latCenter]; - const minCorner = mercator.px([bbox[0], bbox[3]], z), - maxCorner = mercator.px([bbox[2], bbox[1]], z); - const w_ = w / (1 + 2 * padding); - const h_ = h / (1 + 2 * padding); + // calculate zoom + const distancePerPixel = haversine(east, west) / w; + const distancePerTile = distancePerPixel * 256; + const circumferenceEarth = 2 * Math.PI * 6378137; + const z = Math.log2(Math.abs(circumferenceEarth * Math.cos(latRadian) / distancePerTile)); - z -= Math.max( - Math.log((maxCorner[0] - minCorner[0]) / w_), - Math.log((maxCorner[1] - minCorner[1]) / h_) - ) / Math.LN2; - - z = Math.max(Math.log(Math.max(w, h) / 256) / Math.LN2, Math.min(25, z)); - - return z; + // bounds enforcement [0,25] + return Math.min(Math.abs(z), 25) }; const existingFonts = {}; From 29b2187af5596f4cd7eb7e1c709228a3fd26f71f Mon Sep 17 00:00:00 2001 From: iom Date: Thu, 4 Mar 2021 12:52:56 +0100 Subject: [PATCH 2/2] calcZForBBox - add missing distance function Add distance function based on haversine formula (see https://en.wikipedia.org/wiki/Haversine_formula ) --- src/serve_rendered.js | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/serve_rendered.js b/src/serve_rendered.js index 2c84762..1a0664a 100644 --- a/src/serve_rendered.js +++ b/src/serve_rendered.js @@ -27,6 +27,8 @@ const utils = require('./utils'); const FLOAT_PATTERN = '[+-]?(?:\\d+|\\d+\.?\\d+)'; const httpTester = /^(http(s)?:)?\/\//; +const RADIUS_EARTH_EQUATORIAL_IN_M = 6378137; + const getScale = scale => (scale || '@1x').slice(1, 2) | 0; mbgl.on('message', e => { @@ -173,6 +175,19 @@ const renderOverlay = (z, x, y, bearing, pitch, w, h, scale, return canvas.toBuffer(); }; +const toRadian = (deg) => deg*Math.PI/180; + +const getDistanceForPoints = (a, b) => { + // calculates great-circle distance between two points on a sphere in + // meters (see https://en.wikipedia.org/wiki/Haversine_formula) + const phi = [toRadian(a[1]), toRadian(b[1])]; // latitudes + const lambda = [toRadian(a[0]), toRadian(b[0])]; // longitudes + const h = (Math.sin((phi[1]-phi[0])/2))**2 + + Math.cos(phi[0]) * Math.cos(phi[1]) * + (Math.sin((lambda[1]-lambda[0])/2))**2; + return 2 * RADIUS_EARTH_EQUATORIAL_IN_M * Math.asin(Math.sqrt(h)); +}; + const calcZForBBox = (bbox, w, h, query) => { // Use equation for horizontal distance per tile given zoom and latitude // see https://wiki.openstreetmap.org/wiki/Zoom_levels#Distance_per_pixel_math @@ -189,15 +204,14 @@ const calcZForBBox = (bbox, w, h, query) => { // points (use centered latitude of bbox) const latCenter = bbox[1]+(bbox[3]-bbox[1])/2; - const latRadian = latCenter*Math.PI/180; const east = [bbox[0], latCenter]; const west = [bbox[2], latCenter]; // calculate zoom - const distancePerPixel = haversine(east, west) / w; + const distancePerPixel = getDistanceForPoints(east, west) / w; const distancePerTile = distancePerPixel * 256; - const circumferenceEarth = 2 * Math.PI * 6378137; - const z = Math.log2(Math.abs(circumferenceEarth * Math.cos(latRadian) / distancePerTile)); + const circumferenceEarth = 2 * Math.PI * RADIUS_EARTH_EQUATORIAL_IN_M; + const z = Math.log2(Math.abs(circumferenceEarth * Math.cos(toRadian(latCenter)) / distancePerTile)); // bounds enforcement [0,25] return Math.min(Math.abs(z), 25)