How to transform EPSG:3857 to tile pixel coordinates at zoom factor 0?

Question

I'm creating a Leaflet map with an WebGL overlay (using L.CanvasOverlay) which renders a street network on the GPU.

Let's assume I have a point (1497143.62, 6896952.18) in EPSG:3857 and want to transform it to a pixel coordinate on a 256x256 raster, similar to a google maps tile at zoom level 0 (compare: this and that).

I can convert the same point (52.51904, 13.44907) to EPSG:4326 and use the following JavaScript code to do the reprojection to the 256x256 tile raster:

/**
 * Converts Lat/Lon to tile pixel X/Y at zoom level 0 
 * for 256x256 tile size and inverts y coordinates.
 *
 * @param {float} lat Latitude
 * @param {float} lon Longitued
 * @return Point with x and y corrdinates.
 */
function latLonToPixels(lat, lon) {
    var sinLat = Math.sin(lat * Math.PI / 180.0);
    var pixelX = ((lon + 180) / 360) * 256;
    var pixelY = (0.5 - Math.log((1 + sinLat) / (1 - sinLat)) / (Math.PI * 4)) * 256;
    return { x: pixelX, y: pixelY };
}

... which leaves me with the desired point (137.56378, 83.94214) on the raster.

So, using leaflet I could transform my point from 3857 to 4326 and afterwards calculate the pixel coordinates using the formula above. But as I understand the forumlas, this would require two spherical reprojections which is jus too much computing for large geometries.

I suspect there is a much more simple and efficient approach for this transformation. But I couldnt figure it out yet. Any idea?

In addition, searching the internet I can't figure out what the EPSG:3857 coordinates actually denote. This and that answer try to explain it a bit but it is still unclear how to programmatically solve this.

score 7 · Accepted Answer · edited May 22 '19 at 11:56

This is the answer. Sometimes you have to go through all the process of asking a question to understand the solution. JavaScript:

/**
 * Converts spherical web mercator to tile pixel X/Y at zoom level 0 
 * for 256x256 tile size and inverts y coordinates.
 *
 * @param {L.point} p Leaflet point in EPSG:3857
 * @return {L.point} Leaflet point with tile pixel x and y corrdinates.
 */
function mercatorToPixels(p)  {
  var equator = 40075016.68557849;
  var pixelX = (p.x + (equator / 2.0)) / (equator / 256.0);
  var pixelY = ((p.y - (equator / 2.0)) / (equator / -256.0));
  return L.point(pixelX, pixelY);
}

Explanation

The tile at zoom level 0 has a size of 40,075,016m * 40,075,016m. The tile size is 256px * 256px. The center of the web mercator is (0, 0).

Translate the X-coordinate by a positiv half equator.
Divide the X-result by the positive pixel size in meter.
Translate the Y-coordinate by a negative half equator.
Divide the Y-result by the negative pixel size in meter.

That's it.

Please mark it as answer, even though it is yourself (everyone does that, duck typing :) ). — bugmenot123, Jul 09 '15 at 12:39
@bugmenot123 i can not mark my own answer as accepted within 48h after asking the question, but thanks for the heads up. — q9f, Jul 10 '15 at 09:53

score 1 · Answer 2 · answered Jul 09 '15 at 11:36

I would look into using the proj4.js library to transform your coordinates programatically using javascript. From the Github Readme:

Basic usage: proj4(fromProjection[, toProjection2, coordinates])

eg:

```

var firstProjection = 'PROJCS["NAD83 / Massachusetts Mainland",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",42.68333333333333],PARAMETER["standard_parallel_2",41.71666666666667],PARAMETER["latitude_of_origin",41],PARAMETER["central_meridian",-71.5],PARAMETER["false_easting",200000],PARAMETER["false_northing",750000],AUTHORITY["EPSG","26986"],AXIS["X",EAST],AXIS["Y",NORTH]]';

var secondProjection = "+proj=gnom +lat_0=90 +lon_0=0 +x_0=6300000 +y_0=6300000 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";

//I'm not going to redefine those two in latter examples.
proj4(firstProjection,secondProjection,[2,5]);
// [-2690666.2977344505, 3662659.885459918]

It's either not answering my question or so badly explained that I don't understand what you are suggesting, sorry. — q9f, Jul 10 '15 at 09:55
wow @donSchoe thanks for down voting my answer when I was attempting to be helpful. I realize now that you weren't looking for a library, mis-read your question. I linked to the documentation so don't see how that's "poorly explained". No need to be a jerk. — clhenrick, Jul 10 '15 at 14:56

How to transform EPSG:3857 to tile pixel coordinates at zoom factor 0?

2 Answers2