Common Mercator Projection formulas for Google Maps not working correctly

Question

I am building a Tile Overlay server for Google maps in C#, and have found a few different code examples for calculating Y from Latitude. After getting them to work in general, I started to notice certain cases where the overlays were not lining up properly. To test this, I made a test harness to compare Google Map's Mercator LatToY conversion against the formulas I found online. As you can see below, they do not match in certain cases.

Case #1

Zoomed Out: The problem is most evident when zoomed out. Up close, the problem is barely visible.

Case #2

Point Proximity to Top & Bottom of viewing bounds: The problem is worse in the middle of the viewing bounds, and gets better towards the edges. This behavior can negate the behavior of Case #1

The Test:

I created a google maps page to display red lines using the Google Map API's built in Mercator conversion, and overlay this with an image using the reference code for doing Mercator conversion. These conversions are represented as black lines. Compare the difference.

The Results: Equator http://www.kayak411.com/Mercator/MercatorComparison%20-%20Equator.png North Zoomed Out http://www.kayak411.com/Mercator/MercatorComparison%20-%20North%20Zoomed%20Out.png

Check out the top-most and bottom-most lines: North Top & Bottom Example http://www.kayak411.com/Mercator/MercatorComparison%20-%20North%20Zoomed%20Out%20-%20TopAndBottom.png

The problem gets visually larger but numerically smaller as you zoom in: alt text http://www.kayak411.com/Mercator/MercatorComparison%20-%20North%20Zoomed%20Midway.png

And it all but disappears at closer zoom levels, regardless of screen orientation. alt text http://www.kayak411.com/Mercator/MercatorComparison%20-%20North%20Zoomed%20In.png

The Code:

Google Maps Client Side Code:

            var lat = 0;
        for (lat = -80; lat <= 80; lat += 5) {
            map.addOverlay(new GPolyline([new GLatLng(lat, -180), new GLatLng(lat, 0)], "#FF0033", 2));
            map.addOverlay(new GPolyline([new GLatLng(lat, 0), new GLatLng(lat, 180)], "#FF0033", 2));
        }

Server Side Code:

Tile Cutter : http://mapki.com/wiki/Tile_Cutter

OpenStreetMap Wiki : http://wiki.openstreetmap.org/wiki/Mercator

 protected override void ImageOverlay_ComposeImage(ref Bitmap ZipCodeBitMap)
        {
            Graphics LinesGraphic = Graphics.FromImage(ZipCodeBitMap);

            Int32 MapWidth = Convert.ToInt32(Math.Pow(2, zoom) * 255);

            Point Offset =
                Cartographer.Mercator2.toZoomedPixelCoords(North, West, zoom);

            TrimPoint(ref Offset, MapWidth);

            for (Double lat = -80; lat <= 80; lat += 5)
            {
                Point StartPoint = Cartographer.Mercator2.toZoomedPixelCoords(lat, -179, zoom);
                Point EndPoint = Cartographer.Mercator2.toZoomedPixelCoords(lat, -1, zoom);

                TrimPoint(ref StartPoint, MapWidth);
                TrimPoint(ref EndPoint, MapWidth);

                StartPoint.X = StartPoint.X - Offset.X;
                EndPoint.X = EndPoint.X - Offset.X;

                StartPoint.Y = StartPoint.Y - Offset.Y;
                EndPoint.Y = EndPoint.Y - Offset.Y;


                LinesGraphic.DrawLine(new Pen(Color.Black, 2),
                    StartPoint.X,
                    StartPoint.Y,
                    EndPoint.X,
                    EndPoint.Y);

                LinesGraphic.DrawString(
                    lat.ToString(),
                    new Font("Verdana", 10),
                    new SolidBrush(Color.Black),
                    new Point(
                        Convert.ToInt32((width / 3.0) * 2.0),
                        StartPoint.Y));
            }
        }

        protected void TrimPoint(ref Point point, Int32 MapWidth)
        {
            point.X = Math.Max(point.X, 0);
            point.X = Math.Min(point.X, MapWidth - 1);

            point.Y = Math.Max(point.Y, 0);
            point.Y = Math.Min(point.Y, MapWidth - 1);
        }

So, Anyone ever experienced this? Dare I ask, resolved this? Or simply have a better C# implementation of Mercator Project coordinate conversion?

Thanks!

Gwenngwenneth

you should check out sharpmap and proj.net, they are on codeplex – Objection

@Muad'Dib Thanks for the suggestion. I grabbed the Proj.net stuff, but their projection formula doesn't take zoom factor into account, and I can't figure out how to apply it. (projnet.codeplex.com/Thread/View.aspx?ThreadId=77458) It also looks like SharpMap uses Proj.net as their projection engine, so that probably won't be of much help. – Gwenngwenneth

I also found Microsoft's own implementation of Mercator that they use on Bing Maps, and it is also displaying the same problem on Google Maps. I'm almost wondering if it's because I'm using the north west viewing corner as the offset, maybe I need to use the center as my offset... msdn.microsoft.com/en-us/library/bb259689.aspx – Gwenngwenneth

There's been a lot of discussion of this topic. I'm a GIS n00b, so I can't say I understand what's really going on, but try googling "ESPG 900913" for some discussion. – Monumentalize

Answer 1

Thank you all for your suggestions & assistance.

What I eventually found out is that it's not a formula or technical problem, I believe it's a methodology problem.

You can't define the viewing area in Lat/Lng format, and expect to populate it with the appropriate Mercator projections. That's where the distortion happens. Instead, you have to define the correct viewing box in Mercator, and project Mercator.

Doing that I was able to correctly match up with Google maps.

Gwenngwenneth

Answer 2

You may have to create several points along the longtitude in order for the points to be projected corretly along the latitude. In your examples you are only really projecting two points at the start and the end of the line and connecting the two.

The problem will be more apparent at the equator due to the more significant curvature of the earth. It will be less when zoomed in for the same reason.

Have a look at http://code.google.com/apis/maps/documentation/overlays.html#Great_Circles

Try creating your Google polylines with the geodsic parameter to see if this makes a difference. I think this adds points along the line and projects them automatically:

var lat = 0;
var polyOptions = {geodesic:true};
for (lat = -80; lat <= 80; lat += 5) {
    map.addOverlay(new GPolyline([new GLatLng(lat, -180), new GLatLng(lat, 0)], "#FF0033", 2, polyOptions));
    map.addOverlay(new GPolyline([new GLatLng(lat, 0), new GLatLng(lat, 180)], "#FF0033", 2, polyOptions));
}

I had to read into this as all my distance measurements were wrong in OpenLayers for similar reasons: http://geographika.co.uk/watch-out-for-openlayer-distances (more links/explanations)

Slayton