static image rendering moves bbox (wrong zoom level by calcZForBBox)

[cholin]

I'm trying to use the static render export of tileserver-gl through bounding boxes with custom widths+heights. It seems the resulting rendered image does not correspond to my provided bbox but somehow moved a bit. Let's look at the following example:

I have a bounding box which is defined through the following coords [13.383408, 52.483568, 13.404351, 52.492578]. It looks like the following in leaflet:

I render the image through tileserver-gl (v3.0.0)+osm-liberty style with the following link: http://localhost:8080/styles/osm-liberty/static/13.383408,52.483568,13.404351,52.492578/1754x1240.png. For easier comparison, I put the grid on each image. The ratio of the size 1240/1754 corresponds to the the one of the bounding box. This results in the following rendered image:

If you have a deeper look at both images you can see that the rendered image include parts (look at the street on the top right corner) which where not part in the original bounding box. It seems the whole bbox moved a bit to the north. I rendered the same selection as well with mapnik. Here it works as expected (easier to see if you take the mapnik image and the tileserver-gl one and look it in an extra image viewer):

I tried to look for reasons why tileserver-gl behaves this way. Seems you translate the bounding box into a centered point of the box itself, width+height and a zoom level. The latter one is a float and seems calculated in calcZForBBox:

https://github.com/maptiler/tileserver-gl/blob/d5a079d8f4e5bf697286033ee66d3814cc3c0cb0/src/serve_rendered.js#L176

At least for me it is not clear how you come up with the calculation (like why you divide by Math.LN2 for instance). Somehow it looks a bit like magic numbers but maybe I'm just missing something here. On the other hand if you take the distance per pixel equation by osm, transform the equation into z to be able to calculate the zoom level it works. In the following you find a snippet how to do so (you need an additional haversine function to calculate the distance between two points in WSG84):

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
  //
  // 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) )

  // 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 distancePerTile = distancePerPixel * 256;
  const circumferenceEarth = 2 * Math.PI * 6378137;
  const z = Math.log2(Math.abs(circumferenceEarth * Math.cos(latRadian) / distancePerTile));

  // bounds enforcement [0,25]
  return Math.min(Math.abs(z), 25)
};

This produces the following fixed static image:

Unfortunately it seems the current implementation is buggy. If you are interested in the snippet, I could create as well a PR. Otherwise I would be happy if you can tell me a way how to circumvent this behavior to get the bbox selection I want for rendered static images. :)

Edit: You may think you can just use the Mercator coordinates to calculate the distance. This would be easier as they are already in Meters. I tried this first but because of the deviation when you are far away from the equator, you end up with worse results. Here is a mercator snippet:

const latCenter = bbox[1]+(bbox[3]-bbox[1])/2;
const latRadian = latCenter*Math.PI/180;
const left = mercator.forward([bbox[0], latCenter]);
const right = mercator.forward([bbox[2], latCenter]);

const distancePerPixel = Math.abs(left[0]-right[0]) / w;
const distancePerTile = distancePerPixel * 256;
const circumferenceEarth = 2 * Math.PI * 6378137;
const z = Math.log2(Math.abs(circumferenceEarth * Math.cos(latRadian) / distancePerTile));

For my bbox it calculates the following distances in Meter:

• mercator: 2331.3640956836753
• haversine: 1418.0409892023902

You see there is already a big difference for points of small distances in Berlin. If you look at the following rendered image with Mercator distance, you see it is even worse than the current implementation:

In my opinion calculating the zoom level with haversine as distance function seems to be the fix.

[cholin]

I looked again at the current implementation of calcZForBBox. I hope I got it right:

It seems you calculate the distance in Meters of the bbox in Mercantile coordinates of either width or height. Then you derive the distance per pixel by dividing again through width or height. Then you apply log2 through a change of base to the natural logarithm (a bit complicated to read).

So I think the problems are:

• Mercantile distance measurements are vague for points above or below the equator. They even get worse the more you are away.
• No consideration of latitude and the implication of it (distortion of tiles)