I see a lot of people get confused or report problems when attempting to project spatial data to create raster images that line up with the base map layers used in Bing Maps / Google Maps. It’s not helped by the fact that there is a lot of misinformation on the net on the subject. I hope that this post goes someway to helping explain the confusion, and not adding to it…
Note that I’m not describing the process of cutting up tiles or anything to do with quadkey numbering systems – these topics are already described well in this MSDN article. Instead, I’m just talking about how you project geographic coordinate data defined on a 3d model of the world, such as this:
Onto the flat, square map used as the base map layer in Google Maps / Bing Maps, such as this:
Once you’ve done this projection onto a flat image, then you can get on with the business of chopping it up into smaller tiles as appropriate.
The Problem(s)
This seems like it should be a relatively straightforward projection, and there are plenty of spatial libraries out there that can perform transformation and projection between geographic and projected coordinate systems, so why does this particular conversion seem to present so many problems?
I think there are three main reasons:
1. The sort of people that visualise data on Google Maps / Bing Maps are (typically) not GIS professionals. This comment is not in anyway meant to belittle the quality of their data or output. However, the fact is that most GIS professionals still use the commercial (expensive) ESRI ArcGIS product set, and free tools such as Google / Bing are always more likely to attract casual spatial developers. As such, these developers are less likely to know (or want to know) about things like datums, coordinate systems, projections – they just want a map that works.
2. The EPSG code assigned to the projection has changed several times and much of the information displayed on the internet, as well as used in spatial applications themselves, is out-of-date. You’ll see references to EPSG 900913, 3857, 3785, 3587, and various different definitions of each of those systems.
3. The actual projection used by Microsoft / Google itself is (a bit of) a hack, which can lead to misalignment issues when your code actually does the projection the “correct” way. For example, it is common to find issues raised across a number of tools that conversion between WGS84 and Bing / Google Maps leads to Y values that are “out” by about 20km – see the Proj.NET issues here, here, here, here, and here, the ArcGIS API for Silverlight issue here, or the issue for the Proj.4 library reported here, for example.
To provide a practical demonstration of this problem, consider the following location, expressed in WGS84 (latitude, longitude) coordinates:
(52.4827802220782, -5.625 )
This point, incidentally, is the WGS84 coordinate of the point that lies at the bottom left hand corner of Bing Maps’ quadkey 031311. If you project this location into the “official” Google /Bing Maps projection (EPSG:3857) using the Proj.NET reprojection library, you’ll get the following result:
(-626172.13571216376, 6853955.508199729)
You can confirm this by using the online tool at http://www.synectics-tc.com/resources/downloads/coordinate-reprojection.html which is based on the Proj.NET library. (Note, however, that this still uses the outdated 3785 code to refer to the projection).
However, if you transform this point using, say the latest version (that’s important) of the CS2CS tool, as follows:
echo -5.625 52.4827802220782 | cs2cs -f “%.10f” +init=epsg:4326 +to +init=epsg:3857
then you’ll get the following result instead:
(-626172.1357121646, 6887893.4928337997)
Notice that, whilst the X coordinate remains the same, the Y coordinate differs by something like 35km between the two conversions. Projecting a simplified polygon of the United Kingdom using the same projection parameters in these two tools results in differences as shown in the following images:
This same discrepancy appears across a whole variety of conversion tools – some returning the first result, and some returning the second (which lines up correctly with Bing / Google Maps). So, which is right?
The Explanation
Issue 1.) really needs no explanation – as spatial data becomes more mainstream, it is inevitable that new developers will want to start doing things with spatial data, and displaying it on a Google Map / Bing Map is typically one of the very first things you might attempt. These tools both offer great, intuitive UIs to navigate and analyse spatial data, and seeing as WGS84 is by far the most prolific source of data, it seems reasonable to expect that conversion from WGS84 <-> Bing / Google should not be too hard to do.
As for the other issues…
Confusion over the EPSG Code
Coordinate systems are defined using a set of parameters, including the type of projection used (Mercator, Albers etc.), coordinate offset (False Easting, False Northing), the unit of measurement (Metre, Foot etc.), the underlying datum (WGS84, NAD27 etc.) and many more.
Because it’s quite a pain to have to keep quoting this complete set of parameters, every system is commonly referred to by a unique integer code (also called it’s SRID) instead. The set of SRIDs assigned to spatial reference systems in mainstream usage was originally set up by the European Petroleum Survey Group, and they are generally referred to as EPSG codes.
“EPSG:” 900913
Developers started projecting their own spatial data to overlay on Google and Bing Maps ever since their APIs allowed you to do so, which was back in about 2005. Unfortunately, however, EPSG were rather dismissive of the system used by Microsoft and Google (more details in a bit) and refused to assign it an official EPSG code. An EPSG statement (as reported by Morten Nielsen ) read:
“We have reviewed the coordinate reference system used by Microsoft, Google, etc. and believe that it is technically flawed. We will not devalue the EPSG dataset by including such inappropriate geodesy and cartography.”
As a result, during the first few years, developers needed some other way to refer to the set of projection parameters used on the Google Map, and the code 900913 (i.e. GOOGLE) came into common usage. It’s worth noting that this code was never developed, nor supported, by the EPSG (hence the quotation marks on the heading to this section).
The problem that arose was, with no official body to oversee the standard, many different variations appeared – all generally using the same parameters, but with different naming conventions etc., and these had to be manually added into any databases that were otherwise generally populated with official recognised systems from the EPSG.
EPSG: 3785
EPSG took some time to warm to the system used by Microsoft and Google, but eventually, in mid-2008 they decided to add the system to their code registry, assigning it the code EPSG: 3785. There was a small fanfare among the spatial community, and people started recoding the old 900913-hack code to use the new “official” code instead.
The set of projection parameters corresponding to EPSG:3785 was as follows:
PROJCS["Popular Visualisation CRS / Mercator",
GEOGCS["Popular Visualisation CRS",
DATUM["Popular Visualisation Datum",
SPHEROID["Popular Visualisation Sphere", 6378137, 0, AUTHORITY["EPSG",7059]],
TOWGS84[0, 0, 0, 0, 0, 0, 0],
AUTHORITY["EPSG",6055]],
PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]],
AXIS["E", EAST],
AXIS["N", NORTH],
AUTHORITY["EPSG",4055]],
PROJECTION["Mercator"],
PARAMETER["False_Easting", 0],
PARAMETER["False_Northing", 0],
PARAMETER["Central_Meridian", 0],
PARAMETER["Latitude_of_origin", 0],
UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],
AXIS["East", EAST],
AXIS["North", NORTH],
AUTHORITY["EPSG",3785]]
(As an aside, the projection was listed on the EPSG registry with an accompanying note reading “It is not a recognised geodetic system”).
EPSG: 3857
However, less than a year after the official EPSG: 3785 code was finally introduced, EPSG issued a change request (EPSG::2008.114, issued in Dec 2008), deprecating the 3785 code and introducing a new code, 3857, with corresponding parameters as follows:
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich", 0, AUTHORITY["EPSG","8901"]],
UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre", 1, AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
AUTHORITY["EPSG","3857"]]
It’s not unusual for EPSG codes in the registry to be updated every now and again, but this change was particularly odd – firstly because it happened so soon after the original code had been introduced, and secondly because the “replacement” code 3857 was such a close permutation of the previous code 3785 (coincidence?).
Whereas I can easily remember many of the common EPSG codes I use (4326, 27700, 4269 etc.), I always get muddled up as to which one of 3857 and 3785 is the new and which is the old, and many sites on the internet incorrectly list the parameters of one under the code of the other. Also, despite the fact that EPSG registry normally retains details of deprecated codes, 3785 has mysteriously vanished from the database completely, with the only note in the database being “OGP revised its approach to description of Popular Visualisation CRS.”….
The Projection
Once we’ve solved the mystery over which EPSG code is which, we then come on to the problem of understanding the discrepancies between the different results of projection which the various codes above represent.
The projection parameters defined above for EPSG:3785 describes a Mercator projection of geographic coordinates defined on a spherical model of the earth, of radius 6,378,137m, as specified on the following line:
SPHEROID["Popular Visualisation Sphere", 6378137, 0, AUTHORITY["EPSG","7059"]]
The EPSG:3857 projection, in contrast, suggests a Mercator projection of coordinates defined on the WGS84 ellipsoid, with an inverse-flattening ratio of 298.257223563:
SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG","7030"]]
The problem is that neither of these really describes the Google / Bing Maps projection, in which the underlying geographic coordinates are defined using WGS84 (as in 3857), but projected as if they were defined on a sphere (as in 3785).
The problem is that there is no standardised way to represent this “dual”-ellipsoid system, in which a different ellipsoid is used to interpret the geographic coordinates than is used to project them. Different software libraries handle this anomaly in different ways, which explains the differences in the coordinates obtained from using this projection in various tools.
Systems that treat the supplied coordinates as measured on a sphere will be incorrect, but so will systems that project the WGS84 ellipsoid directly using Mercator projection. In fact, this second class of errors explains the apparent “20km out” in the y coordinate only experienced in many of the posts linked to earlier. Since the WGS84 ellipsoid is oblate (like a sphere that’s been sat on) and less high than the sphere, features on a Mercator map projected directly from the WGS84 ellipsoid will appear to be closer to the equator than those projected from a sphere (i.e. placed too far south in the Northern Hemisphere, or too far North in the Southern Hemisphere).
The Solution
Firstly, you should start by making sure that you’re supplying the correct EPSG code, 3857, and that it corresponds to the projection parameters listed above.
However, unless you’re using a conversion library that has been specifically designed to handle the special case of the Bing Maps / Google Maps projection, you’re also going to have to provide a little bit more information to explain that the projection should occur from the sphere rather than the WGS84 ellipsoid specified in the geographic coordinates system. That additional information can be supplied in different ways:
In the PROJ.4 library, you can specify a null gridshift using the +nadgrids=@null parameter as explained here. In fact, the reason why, earlier this post, I specified that you need to use the latest version of CS2CS is that there was a change made in version 4.7.0 of the PROJ.4 library. This change means that that the +nadgrids=@null parameter is automatically added to the internal definition of the 3857 projection, which means that you no longer need to manually include it every time -you’ll get the correct result using the EPSG 3857 projection as-is.
If you’re using geoserver, you need to edit the projection parameters for 3857, replacing the projection from “Mercator” to “Mercator_1SP_Google”.
If you’re using any other system, you need to add an additional “semi_minor” parameter to the projection parameters in the WKT. This parameter allows you to manually specify the semi-minor axis of the ellipsoid used for the projection, which by default is otherwise assumed to be the same ellipsoid specified in the GEOGCS.
To create a projection based on the sphere that has the same radius as the WGS84 ellipsoid, as Google/Microsoft use, you should therefore state a semi-minor axis of 6,378,137. (Note – you can also provide a semi-major parameter if required, but in this case we want a sphere based on the same semi-major axis as the WGS84, so it is not necessary). This will override the WGS84 semi-minor axis that would have been used for the projection otherwise – 6,356,752.314245m.
So, the following WKT works, for example in Proj.NET, to specify the correct projection to Bing / Google Maps.
PROJCS["Popular Visualisation CRS / Mercator",
GEOGCS["Popular Visualisation CRS",
DATUM["WGS84",
SPHEROID["WGS84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7059"]],
AUTHORITY["EPSG","6055"]],
PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]],
AXIS["E", EAST], AXIS["N", NORTH], AUTHORITY["EPSG","4055"]],
PROJECTION["Mercator"],
PARAMETER["semi_minor",6378137],
PARAMETER["False_Easting", 0],
PARAMETER["False_Northing", 0],
PARAMETER["Central_Meridian", 0],
PARAMETER["Latitude_of_origin", 0],
UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],
AXIS["East", EAST], AXIS["North", NORTH],
AUTHORITY["EPSG","3785"]]
The final option, of course, is to forego the general projection routines and write your own specific mathematical transformation designed to convert only between Bing / Google Maps and WGS84. Thankfully, the decision to operate on a perfect sphere makes these functions very easy to write, using only a bit of spherical trigonometry:
double[] WGS84toGoogleBing(double lon, double lat) {
double x = lon * 20037508.34 / 180;
double y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
y = y * 20037508.34 / 180;
return new double[] {x, y};
}
double[] GoogleBingtoWGS84Mercator (double x, double y) {
double lon = (x / 20037508.34) * 180;
double lat = (y / 20037508.34) * 180;
lat = 180/Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);
return new double[] {lon, lat};
}
