The Google Maps / Bing Maps Spherical Mercator Projection

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:

image

Onto the flat, square map used as the base map layer in Google Maps / Bing Maps, such as this:

image

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:

image

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};
}
About these ads
This entry was posted in Bing Maps and tagged , , , . Bookmark the permalink.

57 Responses to The Google Maps / Bing Maps Spherical Mercator Projection

  1. Pingback: Tweets that mention The Google Maps / Bing Maps Spherical Mercator Projection | Alastair Aitchison -- Topsy.com

  2. Morten says:

    “The problem is that there is no standardised way to represent this “dual”-ellipsoid system”
    Ths makes no sense to me. You don’t have a dual ellipsoid. What you have is an input and an output spatial reference. Each SR has each their ellipsoid. In the above case you have WGS84 (epsg:4326) and you have the Mercator projection with a spherical ellipsoid (epsg:3857). You use one for input, and another for output. You will end up with two points (input and output), each representing the same point in the world, but in different “units”.

    The above WKT is also incorrect for 3857. Because its using an ellipsoid that’s a sphere, the 298.257223563 value is wrong, because this makes the WebMercator projection use a datum that has a flattening. Instead enter “0″ here to make it a sphere. The “semi_minor” parameter for Mercator is not really needed, and I doubt most projection engines would know what to do with that (they should pull the value from the ellipsoid instead). The ones that use it I think is using a hack for WebMercator, instead of properly implementing the datum-transform. In this case the datum transform is really easy: you don’t really have to do anything, since it shares same scale, center and orientation as WGS84, so the TOWGS84 parameter would contain all zeros. However in many cases the TOWGS84 parameter is important, since this tells the projection engine how to perform the datum shift (often they will just skip this step without it, kinda like what happens with proj.4 without the nadgrids ‘hack’)

    • alastaira says:

      Thanks for the comment, Morten – so much for me clearing the confusion, eh? :) The point I’m trying to make with the “dual-ellipsoid” comment is that the Bing Maps tiles are projected based on a different ellipsoid from that in which the coordinates are supplied (specifically, the “output” ellipsoid is a perfect sphere) – perhaps your phrasing is better than mine.
      However, I disagree with your comment about EPSG:3857 – the “input” geographic spatial reference system used by Bing / Google *is* WGS84, and therefore needs to be stated with the appropriate inverse-flattening ratio (as it is in the latest EPSG database). However, as the comment in the EPSG database states, EPSG:3857 uses “spherical development of ellipsoidal coordinates.” – in other words whilst the “input” is based on the WGS84 ellipsoid, the “output” spatial reference system is based on the sphere.

      Since projection libraries assume the ellipsoid used in the “input” spatial reference is the same as that in the “output” spatial reference, you need a hack to override this:- +nadgrids is one way to do this, semi_minor is another.

      Finally, as empirical support for my argument, my definition above for EPSG:3857 (including the “semi_minor” parameter) was the only way I could get Proj.NET to correctly project from WGS84 to EPSG:3857. If you can show me another way of achieving the correct result whilst stating 0 flattening for the input datum, I’d be very happy to be proved wrong :)

  3. Pingback: Loading Open Street Map POIs with SQL Server / Bing Maps | Alastair Aitchison

  4. Pingback: SQL Server : Type geographique - [Partie 2] Introduction précise du type géographique , Blog technique de Nicolas Boonaert

  5. Daniel harman says:

    I tried using the suggested projection for proj.net above and it gave me results orders of magnitude out. However when I switched to this wkt:

    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.0174532925199433, AUTHORITY["EPSG", "9122"]], AUTHORITY["EPSG", "4326"]]

    Which I took from

    http://www.synectics-tc.com/resources/downloads/coordinate-reprojection.html

    I got the same results as on that site, which are very close to:

    http://www.nearby.org.uk/coord.cgi?p=KT1+1AA&f=conv

    My test data was:

    “KT1 1AA”
    Original:518032,169094
    Proj.net:-0.304339270784791,51.4085960537841

    Whereas with the suggested WGS84 I got:

    “KT1 1AA”
    Original:518032,169094
    TDPG.GeoCoordConversion.GridReference
    Proj.net:-33878.8926521592,6693890.38230938

    Would love to know what the issue is, but failing that this info may help someone else trying to convert from the uk codepoint postcode database to google maps.

    • alastaira says:

      I don’t see any issue here… you’ve started with a set of coordinates from the OSGB National Grid system, which is EPSG:27700. You’ve then converted it into WGS84 geographic coordinates (EPSG:4326) in your first example, and Spherical Mercator projection in your second (EPSG:3857).
      (518032, 169094) (EPSG:27700) =
      (-0.304339270784791, 51.4085960537841) (EPSG:4326) =
      (-33878.8926521592, 6693890.38230938) (EPSG:3857)

      These are the same results I get using either CS2CS or Proj.NET with the WKT listed in this post.

  6. Daniel harman says:

    Ah thanks for clearing that up I realise I was confused thinking google maps used WGS84.

  7. Charles says:

    The published EPSG:3857 WKT text works for me too. I took a point in WGS84 long/lat and entered it into Bing Maps. I then converted it using both the published WKT and your revised version (with the semi-minor axis addition). I plotted both as single points on the map (in code), and the former appeared in the same place as the long/lat. The latter appeared some distnace due north.

  8. alfamale says:

    thanks ! you are the man

  9. GoldMelodyVN says:

    Great !. But how to convert lat=-90 to mecartor. ProjNet throw AgumentException: “Transformation cannot be computed at the poles.”
    In WGS84toGoogleBing function if lat=-90 –> log(0) = NA ?

    • alastaira says:

      Yep – you can’t plot a point at the poles in the Mercator projection. In fact,, Google/Bing Maps are both clipped so that the latitude only ranges between -85.05112878 to 85.05112878.

  10. Dennis Geasan says:

    ESRI has a projection named “Web Mercator Auxiliary Sphere” with WKID 102100. They are moving to this as their coordinate system in maps published in ArcGIS Online. 102100 is not listed in the EPSG database. 3857 is not included in the ESRI library of coordinate systems. Are these the same coordinate system?

    • alastaira says:

      Yeah – just to add even more confusion on the issue, ESRI appear to have created their own codes that mirror those of the EPSG. Comment on the issue from Melita Kennedy, Senior Product Engineer at ESRI:

      “The WKID was changed at version 10 to 3857. This is because this version of “web Mercator” closely aligns with the EPSG dataset’s entry (http://www.epsg.org). Doc’s out-of-date, I’ll ask that it be updated.
      The multiple codes are due to release timings:
      A. 102113 (because that was the only way we could emulate the web services’ definition)
      B. 102100 added next, using a differently built, but *completely* equivalent, definition
      C. EPSG then added 3857 (and deprecated their earlier 3785 entry which was closer to how 102113 was defined)
      D. We changed the ESRI 102100 code to EPSG 3857.”

      Taken from the ESRI support forums at http://forums.arcgis.com/threads/8762-WGS-1984-Web-Mercator-(Auxiliary-Sphere)-WKID-102100

      But yes, as far as I’m aware EPSG:3857 is completely equivalent to ESRI:102100

  11. Brian says:

    Great info. So where could I find the latest PROJ4 txt for EPSG 3857?

  12. Aline says:

    Hello,
    I am accessing the table sys.spatial_reference_systems, but there is not SRID 3857 … I use SQL Server 2008 Enterprise. Can you help me?

    Thank you.

    • alastaira says:

      @Aline – there won’t be – SQL Server’s sys.spatial_reference_systems table stores details of geographic coordinate systems only. To convert to a projected coordinate reference system such as 3857 you’ll have to use an external tool. OGR2OGR or Safe FME might help.

  13. Aline says:

    Hello,
    I need to convert the SRID 4326 to 3857. But I can not find the SRID 3857 sys.spatial_reference_systems in the reference table. Could you help me?

    Thank you

  14. Hello,
    i am one of the so-called “casual spatial developers”. i am trying to get the most functionality without the pain of understanding the details. I understand the concepts of datum and projection, but the least i need to know, the more i can do on what interests me (canoe maps).

    To me the 2 pieces of code that convert from Google to WGS84 have a great value, more than all newsgroups i have read about all this. But in addition to that, i have to convert between NAD27, NAD83, GWS84. Do you please have similar recipes ? Otherwise my only option is to check the corrections written on paper maps, and store, map-by-map, these corrections.

    So that would be of great help.

    Thx again !

    Charles

  15. bbirbo says:

    Saved me.
    I was using Proj.NET

  16. Pingback: Conversion of British National Grid (WKID:27700) to WGS84(WKID:4326) and then to WebMercator (WKID:102100)- Draft « osedok

  17. Daniel Góis de Oliveira Gomes says:

    Many thanks for this article Alastair!!
    Now I can go ahead;)

  18. Tom C says:

    Sorry, but you have actually added to the confusion – presumably references in the article to EPSG:3587 are typo errors. It really doesn’t help when you’re explaining issues related to the similarity of 3785 and 3857!

  19. Noel Zinn says:

    There is no possible “decision to operate on a perfect sphere”. Real-world coordinates are ellipsoidal. Spherical computations are a fiction. This is the heart of the difference between EPSG 3785 and 3857, viz. the recognition that the input geographicals are ellipsoidal and that the Web Mercator, though it uses spherical Mercator equations, is a different projection when applied to ellipsoidal geographicals, a non-conformal projection that (IMO) should not be called “Mercator”. -Noel

    • alastaira says:

      Noel – I’m not quite sure I agree with the statement that “spherical computations are a fiction“. And I don’t think the folks at MathWorld would agree either! http://mathworld.wolfram.com/SphericalCoordinates.html

      What I think you’re trying to argue is that the decision to use spherical Mercator equations makes Google Maps / Bing Maps et. al unsuitable for serious scientific purposes, but I don’t think anyone would disagree with that :)

      • Noel Zinn says:

        I agree, Alastair, that, given a sphere, the Wolfram equations are valid, no question. But where is that sphere? It’s a fiction. If you are gathering real-world geographicals with GPS, they are referenced to an ellipsoid, not a sphere. In this case, using the Wolfram equations would induce error. Web Mercator projection of ellipsoidal geographicals is valid, but, unlike the spherical and ellipsoidal Mercators, it’s non-conformal. In that sense Web Mercator wouldn’t be used for surveying, geodetic or (I guess) scientific purposes.

  20. alastaira says:

    But *all* spatial data is only an approximation defined in reference to an ellipsoid model. The “real-world” information you’re gathering using a GPS is an approximation based on the WGS84 ellipsoid. Using a spherical model is a simpler approximation, certainly, but it’s just as valid (after all, a sphere is just a special case of an ellipsoid in which all three axes are equal). Neither one truly reflects the shape of the geoid. So the “fictional” sphere is just as real as any other ellipsoid.

    WGS84 is not the single truth when it comes to geodesy (certainly not here in the UK, where the Ordnance Survey would refuse to let you taint their OSGB36 data by transforming it into “inferior” WGS84).

  21. Noel Zinn says:

    I concede that all real-world measurements have error and, therefore, “*all* spatial data is only an approximation”. Surveyors recognize 3 kinds of error: (1) unavoidable random error that can be minimized as technology improves but which cannot be eliminated, (1) blunders (like a numerical transposition or multipath) that can be eliminated (in theory) with good technique, and (3) bias or systematic error due to an un-calibrated instrument (always reading long, for example) or an incorrect mathematical model (using a sphere instead of an ellipsoid, for example). We get rid of the errors we can (blunders and bias) and use statistics to evaluate the errors we can’t (random). You are confusing random error with mathematical bias, Alastair. The geocenter is 20km closer to the pole than the equator in the WGS84 ellipsoid. That’s a lot of mathematical bias to accept just because “a sphere is just a special case of an ellipsoid” or because your GPS receiver has a few decimeters of random error. BTW, the only reason OSGB36 may fit the geoid better in the UK than WGS84 is that it was locally adjusted (rather than worldwide). The switched-on Brits I know have moved to ETRF89. Time to put your parochialisms aside. BTW … great blog. Thanks.

  22. In ESRI software, it’s easy enough to transform data from WGS84 to any spherical projection – but not the reverse. So the transformation to Web Mercator for Bing/Google Maps is a one way ticket in that environment, for now. Many thanks – best discussion of this topic on the Web.

    • Noel Zinn says:

      This surprises me. Web Mercator is completely reversible LLXY. The whole point of using Web Mercator with WGS84 rather than the ellipsoidal version of Mercator is that the reverse is closed (not iterative) and, thus, faster, though non-conformalities are introduced.

      • Web Mercator is completely reversible between LL and xy, but that’s just the reprojection part of the transformation. Datum shifts are separate, and happen only in geographic coordinates. There is no systematic way to transform data cast on a sphere to WGS84 – no equivalent of the nadcon software for this transformation. But ESRI will allow you to specify this transformation anyway; their Project module will reassign input data from sphere to WGS84 and complete the reprojection. This works just fine as long as the data were originally in WGS84 or another datum – which was the case for a number of datasets I had to transform from sphere to NAD83. I “Asked a Cartographer” at ESRI about this last December: http://mappingcenter.esri.com/index.cfm?fa=ask.answers&q=2233. Still getting my head around Web Mercator, thanks for the comment.

      • Noel Zinn says:

        For John Hutchinson below (no reply button): This quote from the Matrix is helpful in understanding Web Mercator.

        Prodigy: “Do not try to think outside the box. That’s impossible. Simply try to realize the truth.” Neo: “Which is?” Prodigy: “There is no box.”

        The truth is that there is no sphere in Web Mercator. Web Mercator is just a projection that maps from to ellipsoidal LL to XY and back to ellipsoidal LL. Yes, the Web Mercator uses spherical Mercator equations, but THERE IS NO SPHERE. You can still used NADCON to transform from NAD83 (or WGS84) ellipsoidal datum to NAD27 ellipsoidal datum and project on either end using Web Mercator. Web Mercator is just a non-conformal projection and there are many non-conformal projections. The sphere is an illusion. Get your head around that!

  23. Charles says:

    I have copied the following WKT from http://spatialreference.orgref/epsg/29902, but when I try to load it with Proj.net CreateFromWkt(…) I get an exception:

    PROJCS["TM65 / Irish Grid",
    GEOGCS["TM65",
    DATUM["TM65",
    SPHEROID["Airy Modified 1849",6377340.189,299.3249646,
    AUTHORITY["EPSG","7002"]],
    AUTHORITY["EPSG","6299"]],
    PRIMEM["Greenwich",0,
    AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
    AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4299"]],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",53.5],
    PARAMETER["central_meridian",-8],
    PARAMETER["scale_factor",1.000035],
    PARAMETER["false_easting",200000],
    PARAMETER["false_northing",250000],
    AUTHORITY["EPSG","29902"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

    The exception is “Expecting (‘PROJECTION’) but got a ‘UNIT’ at line 12 column 61.”

    Can anyone show me the correct version of the WKT or point me to a good source? I would have thought the site would be correct, or is it Proj.net at fault?

    TIA

  24. Tao zhang says:

    to Noel Zinn:
    impressive comments about “no sphere”, I agree with that

  25. Wes Keller says:

    The figure of the earth is complex. Even the geoid is only an approximation (think local surface roughness, the cell size of computed geoids is very coarse, there are multiple geoids in use). If your mapping the whole world on a mercator projection, a sphere is just as useful as an ellipsoid. At that scale most untrained people using a digital map won’t even notice the difference.

    • Noel Zinn says:

      The irregularity of the geoid is small (+/- 100m) compared with the difference between the sphere and the ellipsoid (about 20km at the poles).

  26. Claire Wordley says:

    Hi – I am a typical GoogleMaps user taht you describe and not a GIS professional – I’m a PhD student just learning to use ArcGIS. I have a problem working between the two – I don’t know how to apply the 3857 code you discussed.

    I’m digitizing features from Google Earth to use in ArcGIS (I work in India, where the mapping in GIS is at best scanty). I drew along a river in Google Earth, saved the line as a .kml file and then used the ‘from KML’ feature in ArcMaps10 to convert to a layer file. It worked, except my line is about 270m north of where it is meant to be. I tried downloading some python scripts from the ESRI website to convert .kml to .shp – they load fine into ArcToolbox but then won’t run my .kml files, citing bad magic numbers. I also tried using an online converter which still left the line north of where it is meant to be.

    All my inherited raster and vector files for the area should be in GCS_WGS_1984 as they were either in this format or I’ve converted them. I’ve also added XY data I took using a GPS and converted from DD.MM.MMMM to DD.DDDD which is all correctly aligned with the raster and vector files I inherited. The .kml file I converted to a layer file says it has the same coordinate system and datum (D_WGS_1984). I’m banging my head against a brick wall trying to see where the problem is! I’ve tried using bing maps as a base layer but the resolution for my area of rural India isn’t high enough, and I tried putting jpg files from google earth in but to get enough detail I’ll have to input millions of tiny files.

    • Noel Zinn says:

      Hi Claire! Since a “real” GIS professional hasn’t stepped in yet, I can offer a couple comments. First, the 3857 projection is native to Google Maps and not Google Earth (which presents as something closer to an orthographic). If you’re pulling your river .kml from Google Earth, you’re in WGS84 geographicals. But, have you set up your GIS projection as 3857? Second, the difference in Northing between Mercator and Web Mercator is way more than 270 meters north of Kanyakumari (that is, anywhere in India). You’d have to be a lot closer to the Equator for that “small” of a discrepancy. So, I think you need to look for a different solution to your problem. Enough to drive you batty, heh? Vanakkam.

  27. Guillermo Biot says:

    Hello,

    I’m having some precision issues between bing and google maps. I’m finding an 17m distance position difference between them for the same lat/lon coordinates. For instance, the center of the round about in this bing map http://binged.it/15fx4s3 gives the coordinates 40.232744, 0.272125. If I copy/paste the same numbers in google maps, the marker is mispositioned by 18m approx.
    Is it normal? I’m a bit confused.
    Thank you in advance.

    Regards,

    Guillermo

    • alastaira says:

      18 metres is not a huge discrepancy in map imagery of the entire world :) And neither Google nor Bing Maps are designed for precise positioning purposes. However, when such discrepancies occur, I often find it helpful to refer to a third source. In this case, Open Street Map shows that location to be more consistent with Bing Maps – http://www.openstreetmap.org/?mlat=40.232732&mlon=0.272176&zoom=17 – located in the middle of the roundabout. It doesn’t prove that Google is “wrong” and Bing Maps is “right”, but OSM and Bing appear to be more consistent at least.

  28. Yusuf Irzan says:

    Hi, I am trying to convert from WGS84 UTM Zone 50S to OpenStreetMap (EPSG 3857) using Proj.Net.
    My input are :
    X=285970.67279052734 and Y=9100209.1314086914.

    I am using WKT from :
    http://spatialreference.org/ref/epsg/32750/prettywkt/ (for UTM Zone 50S)
    http://spatialreference.org/ref/sr-org/7483/prettywkt/ (for EPSG 3857 / OpenStreetMap)

    Projected X is fine, but projected Y are always NaN.
    Can anyone find solution for kind problem?

    Thanks

  29. Guillermo Biot says:

    Hi alastaira, thank you very much for your answer. I agree that 18m is not a big deal, maybe it’s a matter of precision, but as you comment, a third source matches exactly one of them. Indeed the GIS software I use also matches the bing maps location, so I would dare to say that google maps is wrong by 18m, I wouldn’t expect that “error” from google maps. I tried different positions around that area and all of them have the same amount of offset in the same direction (always East, North is correct), it’s not just random. I found out that this missmatching doesn’t take place in other cities that I checked and maybe it will be corrected on the next render/process of the google maps tiles. Do you think it would be worthy to notify the people working on google maps?

    The big deal was that the missmatching was driving me crazy because I was using the lat/lon coordinates from google maps in that area and they didn’t match my maps on the GIS software. Obviously I thought the problem was at my side, on a datum translation.

    Thank you again for your time.

    Regards,

    Guillermo

  30. Pingback: The Google Maps / Bing Maps Spherical Mercator Projection | :nine

  31. Pingback: “Web mercator” coordinate system: explained! | geo.jot

  32. Pingback: OCAD Adds Co-ordinate System for Nearmap | MapSport CartographicMapSport Cartographic

  33. clayton says:

    I own a mapping SIRGAS 2000 and would like to know what to do for my bing map background is in the same projection.

  34. Cefn Hoile says:

    I did (what I thought was) a simple port of your code to python to use it within kartograph.py, but I found that using it to create a new GoogleBing projection within the kartograph proj module, cloning the Mercator projection at… https://github.com/kartograph/kartograph.py/blob/master/kartograph/proj/cylindrical.py
    …but the resulting coordinates created an upside-down projection.

    I imagine that it’s simply due to a difference in the logic of SVG x/y coordinates compared to .Net. Multiplying the y coordinate by -1 superficially solves the problem.

    The projection added to my local kartograph now looks like…

    # added by CH following http://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/
    class GoogleBing(Cylindrical):
    def __init__(self, lat0=0.0, lon0=0.0, flip=0):
    Cylindrical.__init__(self, lon0=lon0, flip=flip)
    self.minLat = -85
    self.maxLat = 85
    def project(self, lon, lat):
    lon, lat = self.ll(lon, lat)
    x = lon * 20037508.34 / 180
    y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180)
    y = y * 20037508.34 / 180
    y = y * -1 # added because appeared upside down
    return (x, y)

  35. Lars Bilke says:

    Hi alastaira,

    thanks for this post! I am trying your fixed projection with the ogr2ogr-tool. I am converting from EPSG:31468 but the tool says that there is “No translation for Mercator to PROJ.4 format is known”. Please see this page for the full output: https://gist.github.com/bilke/7942425

    It would be very nice if you can help me!
    Thanks
    Lars

    • Lars Bilke says:

      I just looked a little bit at your other posts in this blog and saw that you are working with Unity and terrains as well. To explain my motivation for the comment a little more: I am working in environmental research center and I am trying to combine Unity terrains created with WorldComposer / TerrainComposer (http://www.terraincomposer.com, which rely on Bing map data) with data from GIS. Therefore I try to convert the shape-date into the Bing coordinate system / projection. Do you have any experience in combining Unity with GIS data?

      • alastaira says:

        Yes – I’ve used GIS data in Unity a little – specifically GB Ordnance Survey data. As this is already projected in Eastings/Northings expressed in metres, it was pretty straightforward to import into Unity (which, although not explicitly documented, is generally assumed to have a scale in which 1 unit = 1 metre). I added an offset to all coordinate values so that (0,0,0) in Unity worldspace corresponded to my desired origin (Norwich, which is about 623000,308000 in OSGB) and kept scale the same.

        I’ve not used TerrainComposer at all, but so long as you can export data with a similar flat projection it should be quite straightforward.

  36. Pingback: Displaying WRF data using OpenLayers | Computing.io

  37. Roy Jackson says:

    Thank you for the formula – WGS84 to Google

  38. Pingback: Nautilytics - Custom Tile Server for iOS and Web

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s