Calculating Mercator d3 scaling to match Leaflet zoom level

Question

Given a Leaflet zoom level (0..22), how would I calculate D3 scaling value for geoMercator projection? At zoom=0, the whole world fits within a single tile (256x256). In tiles, the world size is 2^zoom x 2^zoom tiles.

Answer 1

For comparison with d3 scales for Mercator maps:

A mercator map uses this formula:

var point =  [x, Math.log(Math.tan(Math.PI / 4 + y / 2))];  

The scale factor in a d3 map essentially applies the following transform:

   point[0] =  point[0] * k;
   point[1] =  -point[1] * k;

The default map scale for d3 mercator projections is 961/2π, or 961 pixels by 360 degrees.

One of the quirks of the Mercator formula is that the a square view of the "whole" world will actually cap out at ±85.05113 degrees North/South. Web Mercators can push this as they are not conformal projections. Leaflet pushes the extent of the "whole" world to ±89.15 degrees or so North/South.

So, d3 using a proper Mercator and Leaflet using a web mercator means that scale values might not mesh nicely. This answer on GIS.stackexchange will provide more info on discrepencies between the two.

But, you can still align the two (most of the time).

---

One method, as noted by Yurik in the comments below is to use the zoom level and tile size to get a scale factor:

Web map services use a tiled grid where increasing a zoom level quadruples the number of tiles shown. At zoom level one, the world fits in one tile, at zoom level two, the world fits in four square tiles. The general formula for the number of tiles is then:

Number of Tiles = 4^zoomLevel

Most importantly, each time we zoom in, the number of tiles that it takes to go across the map (one row) doubles.

With this as a starting point, we can figure out how to match the two.

A d3 geoMercator uses 961/(2*Math.PI) as the default scale - this spreads the equator's 360 degrees (or 2 radians) over 961 pixels. To set this to work for a tile based layer, we only need to know the tile size and the zoom level.

We need to know how many pixels the equator is spread over, to get that we use:

tileSize * Math.pow(2,zoomLevel)

This gives us the width of all the tiles that circle the equator. Then we can divide that by 2Pi and get our d3 scale:

projection.scale(tileSize * Math.pow(2,zoomLevel) / (2 * Math.PI))

With the differences between a d3 Mercator and a web mercator, there may be issues of distortion depending on where you are located and how far in you zoom, but this should offer good alignment in most situations.

---

Alternatively, we can use the actual corners of the leaflet map to determine an appropriate scale:

D3 features a method that allows a projection to fit a geographic extent to a svg extent: projection.fitExtent(). This method takes geojson or geometries. In this answer I'm using geojson.

.getBounds() will return the leaftlet map extent, so you can create a geojson bounding box fairly easily:

var bbox = {
    "type": "Polygon",
    "coordinates": [
        [
            [mymap.getBounds()._northEast.lng, mymap.getBounds()._northEast.lat], 
            [mymap.getBounds()._northEast.lng, mymap.getBounds()._southWest.lat], 
            [mymap.getBounds()._southWest.lng, mymap.getBounds()._southWest.lat], 
            [mymap.getBounds()._southWest.lng, mymap.getBounds()._northEast.lat],
            [mymap.getBounds()._northEast.lng, mymap.getBounds()._northEast.lat]
        ]
    ]
}

Note that winding order is actually important in d3 - outer rings are counter clockwise

Then all you have to do is set your projection to fit that extent:

var projection = d3.geoMercator()
   .fitSize([600, 400], bbox);

Using a slightly modified leaflet example (changed centering point and zoom), the whole thing looks like this:

	var features = {
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -0.17565250396728513,
              51.510118018740904
            ],
            [
              -0.17586708068847653,
              51.509744084227556
            ],
            [
              -0.17584562301635742,
              51.50871574849058
            ],
            [
              -0.17359256744384766,
              51.50613812964363
            ],
            [
              -0.17204761505126953,
              51.50552375337273
            ],
            [
              -0.16966581344604492,
              51.50481587478995
            ],
            [
              -0.16599655151367188,
              51.50454874793857
            ],
            [
              -0.1624774932861328,
              51.504001132997686
            ],
            [
              -0.16058921813964844,
              51.5039744199054
            ],
            [
              -0.16033172607421875,
              51.50426826305929
            ],
            [
              -0.16013860702514648,
              51.5043884710761
            ],
            [
              -0.16016006469726562,
              51.50465559886706
            ],
            [
              -0.15996694564819336,
              51.50510971251776
            ],
            [
              -0.16282081604003906,
              51.505737450406535
            ],
            [
              -0.16466617584228516,
              51.5058710105437
            ],
            [
              -0.16835689544677734,
              51.50588436653591
            ],
            [
              -0.1705455780029297,
              51.506098061878475
            ],
            [
              -0.17273426055908203,
              51.506672363145654
            ],
            [
              -0.17282009124755857,
              51.50681927626061
            ],
            [
              -0.17468690872192383,
              51.508729103648925
            ],
            [
              -0.17511606216430664,
              51.50999782583918
            ],
            [
              -0.17526626586914062,
              51.510144728231545
            ],
            [
              -0.17565250396728513,
              51.510118018740904
            ]
          ]
        ]
      }
    }
  ]
};



	var mymap = L.map('mapid').setView([51.5, -0.171], 14);

	L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpejY4NXVycTA2emYycXBndHRqcmZ3N3gifQ.rJcFIG214AriISLbB6B5aw', {
		maxZoom: 18,
		attribution: 'Map data &copy; <a href="http://openstreetmap.org">OpenStreetMap</a> contributors, ' +
			'<a href="http://creativecommons.org/licenses/by-sa/2.0/">CC-BY-SA</a>, ' +
			'Imagery © <a href="http://mapbox.com">Mapbox</a>',
		id: 'mapbox.streets'
	}).addTo(mymap);
	

	var svg = d3.select('#mapid')
	  .append('svg')
	  .attr('width',600)
	  .attr('height',400);

	
	
	
	// Create a geojson bounding box:
	var bbox = {
		"type": "Polygon",
		"coordinates": [
			[
				[mymap.getBounds()._northEast.lng, mymap.getBounds()._northEast.lat], 
				[mymap.getBounds()._northEast.lng, mymap.getBounds()._southWest.lat], 
				[mymap.getBounds()._southWest.lng, mymap.getBounds()._southWest.lat], 
				[mymap.getBounds()._southWest.lng, mymap.getBounds()._northEast.lat],
			    [mymap.getBounds()._northEast.lng, mymap.getBounds()._northEast.lat]
			]
		]
    }
	
	
	var projection = d3.geoMercator()
	   .fitSize([600, 400], bbox);
	
	var path = d3.geoPath().projection(projection);
	
	svg.append("path")
	  .datum(features)
	  .attr('d',path);
	

		svg {
			z-index: 10000;
			position: relative;
		}
	

<div id="mapid" style="width: 600px; height: 400px;"></div>

	
	<link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />

    <link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css" integrity="sha512-07I2e+7D8p6he1SIM+1twR5TIrhUQn9+I6yjqD53JQjFiMf8EtC93ty0/5vJTZGF8aAocvHYNEDJajGdNx1IsQ==" crossorigin=""/>
    <script src="https://unpkg.com/[email protected]/dist/leaflet.js" integrity="sha512-A7vV8IFfih/D732iSSKi20u/ooOfj/AGehOKq0f4vLT1Zr2Y+RX7C+w8A1gaSasGtRUZpF/NZgzSAu4/Gc41Lg==" crossorigin=""></script>
	
	<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.9.1/d3.js"></script>

I have not added the code to update the projection when the map changes, but this is simply recalculating the bounding box and re-applying the fitSize method.

Answer 2

I am trying to do this from first principles.

This information has given me a great start. Thank you.

Based on your advice, I can see at level 18 there are 0.00137329 degrees per tile.

if I take a known point: x: 234472, y:157536 and it's gps coordinates -34.125714, 141.999165

234472*0.00137329 = 321.9983 or 141.9983 south. EXCELLENT!

however, if I do the same with Y

157536*0.00137329 = 216.3428 or -36.3428, which is 2.22 degrees away from where it should be.

What am I doing wrong? Similar errors are found for other locations

(PS, if I use the extents of the map, X is correct, but Y is wrong, but much closer - I'd rather not use the extents if I didn't have to)

Please help.

KudaPucat